new;
/* proc to generate VAR(1) y(t)=A*y(t-1)+v(t) from A,*/

proc(1) = genvar(A, y0, V);
	local lam, B, Bstar;

	{lam,B} = eigv(A);
	if scalerr(lam[1]); "eig fails in genvar "; stop; endif;

	Bstar = inv(B)';
	
	retp( real(recserar( V*Bstar, y0'Bstar, lam')*B') );
endp;

n = 100;
lam = {0.9, 0.8, 0.7}; 
p=rows(lam); 
D = diagrv(eye(p),lam);D;
B = {1 1 1, 0 1 1, 1 0 1}; 
A = B*D*inv(B); A;
y0 = rndn(p,1); 
V = rndn(n,p);
Y = genvar(A, y0, V);
meanc(Y)~stdc(Y);