/* PROC FINDVAR finds VAR representations of lower-dimensional vectors. given: 1) state vector evolves as x(t) = ax(t-1) + ce(t) E(ee') = I 2) series you want are expressed as n'x; usage: {B,S} = findvar(n,A,C); inputs: n : y = n'x A,C : x(t) = Ax(t-1) + Ce(t); E(ee') = I ouptuts: B: y = B'y(-1) + d S: E(dd') = S note: for now it ignores explosion in the E(xx') matrix under unit roots. get recursive formulas for finite functions of E(xx'), like B and S that don't essentially divide infinity by infinity. --------------------------------------------------------------------------- */ proc(2) = findvar(nl,Al,Cl); local V,j,B,S,n,A,C,CC,An; n = nl; a=al; c=cl; if rows(n) /= cols(A); "findvar: n and A have inconsistent dimensions"; endif; CC = C*C'; V = CC; An = A; j = 1; do while (maxc(maxc(V)) < 1E6) and (j <= 100); /* use doubling algorithm */ V = An*V*An'+ V; An = An*An; j = j+1; endo; B = inv(n'*V*n)*(n'V*A'n); S = (n'A-B'n')*V*(A'n-n*B) + n'C*C'*n; retp(B,S); endp; /* test call*/ /* {B,S} = findvar(ny~nc,Am,Cm); B; S; */ /* note: the following code produces correct B for nonsingular systems (e.g. y not y and c ) when there are unit roots. */ /* {va, ve} = eigv(A) ; /* ve*diagrv(0*A,va)*inv(ve) produces eigenvalue decomposition */ lbig = diagrv( zeros(rows(A),cols(A)), va .== 1); V = ve*lbig*inv(ve)*C*C'*inv(ve)'*lbig*ve'; B = inv(n'*V*n)*(n'V*A'n); */