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;
proc(1) = genvars(A, y0, V);
	local T, i, Y;
	T = rows(V);
	Y = zeros(T, cols(V));
	Y[1,.] = y0';
	i=2;
	do until i > T;
		Y[i,.]= Y[i-1,.]*A' + V[i,.];
		i = i+1;
	endo;

	retp( Y );
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);

nreps=1000;
tim = hsec;
i = 1;
do until i > nreps;
	Y = genvar(A, y0, V);
	i = i+1;
endo;
tim = hsec - tim;

meanc(Y)~stdc(Y);

tim2 = hsec;
i = 1;
do until i > nreps;
	Y2 = genvars(A, y0, V);
	i = i+1;
endo;
tim2 = hsec - tim2;

"times " 0.01*(tim~tim2);
"ratio " tim2/tim;
"last rows" Y[n,.]|Y2[n,.];