/* PROC FINDV2 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; thjei version uses recursive update of B matrix. usage: {B,S} = findv2(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) = findv2(nl,Al,Cl); clearg V,j,B,S,n,A,C,CC,An,t1,t2,d,dinv,chk; n = nl; a=al; c=cl; if rows(n) /= cols(A); "findv2: n and A have inconsistent dimensions"; endif; T1 = A*c*c'*a'; B = inv(n'*(c*c'+t1)*n)*(n'*(c*c'+t1)*a'*n); "B" B; chk= c*c'+ A*c*c'*a'; "check" inv(n'chk*n)*(n'*chk*a'*n); wait;; Dinv = inv(n'*t1*n)*(n'A*c*c'*a'*n); j = 1; do while j <= 1; t2 = a*t1*a'; "t2" t2 ;"check" a*a*c*c'*a'*a'; wait;; dinv = inv(n't2*n)*(n'*t1*n)*dinv + eye(cols(n)); d = inv(dinv); "d" d; wait;; t1 = t2; B = (eye(cols(n))-D)*B + D*inv(n'T1*n)*(n'*t1*a'*n); format /RZ 15,6;"new B" ; B; chk= c*c'*a*chk*a'; "check" inv(n'chk*n)*(n'*chk*a'*n); wait;; j = j+1; endo; "new B" ; B; 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); "old B" ;B; S = (n'A-B'n')*V*(A'n-n*B) + n'C*C'*n; retp(B,S); endp; /* test call*/ {B,S} = findv2(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); */