/* PROC FINDVL2 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; THis version lets you have lagged variables in the VAR This version accomodates unit roots better. usage: {B,S} = findvl2(n,A,C,lags); inputs: n : y = n'x (multiple series ok, e.t. ny~nc) A,C : x(t) = Ax(t-1) + Ce(t); E(ee') = I lags: number of lagged y on rhs of desired var ouptuts: B: y = B'[y(-1) y(-2) ... y(-lag)]' + 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. also, is inelegant: calculates Exx' three times! --------------------------------------------------------------------------- */ proc(2) = findvl2(nl,Al,Cl,lagsl); clearg n,A,C,lags,oneva,CC,V,An,j,Ezz,i,newcol,apj,aj,newrow,Ezz1,Bbig,B,Exx,S, oneigv; "locals cleared!"; @local n,A,C,lags,onveaCC,V,An,j,Ezz,i,newcol,apj,aj, newrow,Ezz1,Bbig,B,Exx,S oneigv;@ n = nl; /* this lets you see inputs while debugging*/ a=al; c=cl; lags=lagsl; if rows(n) /= cols(A); "findvlag: n and A have inconsistent dimensions"; endif; {va,ve}=eigv(A); oneva = ( abs(real(va)-1) .< 1E-10 ); /* vector of unit eigenvectors */ if maxc( oneva ); "findvl2: found unit eigenvalues in A"; if maxc(va) > 1; " Biggest eig. of A is > 1"; endif; /* limit of A^j*c is eig. value corresp to 1; */ oneigv = selif(ve',oneva)'; V = oneigv*oneigv'; /* note v is not EEx but lim T->inf eex/T */ /* times some arbitrary constant */ else; "findvl2: no unit roots in A"; /* find E(xx') */ 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; endif; v = (v+v')/2; /* split the differentce keep it symmetric */ /* construct E(zz') where z = [z z(-1) z(-2) .. ] */ Ezz = n'V*n; i = 2; do while i <= lags; newcol=n'v*n; aj = A; j = 2; do while j <= i; newcol = (n'Aj*v*n)|newcol; Aj = A*Aj; j = j+1; endo; newrow=n'*v*A'*n; apj = A'; j = 3 ; do while j <= i; Apj = A'*Apj; newrow = (n'*v*apj*n)~(newrow); j = j+1; endo; Ezz = (Ezz|newrow)~(newcol); i = i+1; endo; ezz= (ezz+ezz')/2; /* construct first row of E(zz(-1)') where z = [z z(-1) z(-2) .. ] */ aj = a; Ezz1 = n'a*V*n; i = 2; do while i <= lags; aj = a*aj; Ezz1 = Ezz1~(n'aj*V*n); i = i+1; endo; B = inv(Ezz)*(Ezz1'); /* find variance covariance matrix */ /* construct E( [x(t) x(t-1)] ... [ x(t) x(t-1)...]' ) */ Exx = V; i = 2; do while i <= lags; newcol=v; aj = A; j = 2; do while j <= i; newcol = (Aj*v)|newcol; Aj = A*Aj; j = j+1; endo; newrow=v*A'; apj = A'; j = 3 ; do while j <= i; Apj = A'*Apj; newrow = (v*apj)~(newrow); j = j+1; endo; Exx = (Exx|newrow)~(newcol); i = i+1; endo; if lags > 1; Bbig = B~ ( eye(cols(n)*(lags-1)) | zeros(cols(n),cols(n)*(lags-1)) ); else; Bbig=B; endif; S = eye(lags).*.(n'A) - Bbig'(eye(lags).*.(n')); S = S*Exx*S'; S = S[1:cols(n),1:cols(n)]; S = S+ (n'CC*n); retp(B,S); endp; /* test call*/ output file=findvl2.out reset; format /RD 10,2; lags = 2; {B,S} = findvl2(nc~ny,Am,Cm,lags); "B'" B'; "S " S; @ Bbig = B~ ( eye(cols(n)*(lags-1)) | zeros(cols(n),cols(n)*(lags-1)) );@ Bbig = B~ ( eye(2*(lags-1)) | zeros(2,2*(lags-1)) ); "eigenvals of B" real(eig(bbig)'); output off;