/*--------------------------------------------------------------------------- PROC FINDVK finds Wold representation using kalman filter. Model: x(t) = Ax(t-1) + Ce(t) E(ee') = I z(t) = n'x(t). Method: P(z(t)|z(t-1),...) = n'Ax^(t-1) x^(t) = (I-Kn'A)x^(t-1) + Kz(t) K = (ASA' + CC')n[n'(ASA' + CC')n]^(-1) S = (I-Kn')(ASA' + CC')(I-Kn')'; Usage: {ir,sigvar} = findvk(n,A,C); Returns: ir = impulse response sigvar = variance covariance matrix of z wold innovations ----------------------------------------------------------------------------- */ proc(2) = findvk(nl,al,cl); local n,nsm,a,c,sig,term,j,signew,dettest, test,k,sigvar,q,ir,xhat,shock,z,t, maxsit,convtol,irlags,sigtol,sigvtol; maxsit = 1000; /* sigma iterations? */ convtol = 1E-9; /* sigma convergence tolerance */ sigtol = 1E-6; /* test for sigma = 0 */ irlags = 300; /* lags in impulse-response */ sigvtol = 1E-10; /* det of var vcv flags singular */ n = nl; a = al; c = cl; nsm = n; /* possibly smaller n to use if proj(x|z) is singular*/ /* try sigma = 0 solution first */ term = c*c'; nsm = n; do while det(nsm'*term*nsm) < 1E-8; "findvk: pruning n in zero solution"; nsm = nsm[.,1:cols(nsm)-1]; endo; test = maxc(maxc(abs(term-term*nsm*inv(nsm'*term*nsm)*nsm'*term))); if test < sigtol; "sigma=0 solution works"; sig = 0; else; /* loop on ricatti equation to find sigma */ sig = eye(rows(a)); nsm = n; j = 1; do while j <= maxsit; term = A*sig*A' + c*c'; dettest = det(nsm'*term*nsm); if dettest < 1E-8 ; "findvk: warning. Sigma loop trying to invert near-singular mx"; "pruning a column of n "; nsm = nsm[.,1:cols(nsm)-1]; sig = eye(rows(a)); j = 1; continue; endif; signew = term - term * nsm * inv(nsm'*term*nsm) * nsm' * term'; signew = (signew+signew')/2; test = maxc(maxc(abs(signew-sig))); if test < convtol; break; endif; if j >= maxsit-10; "findvk: Sig loop did not converge.last test" test; endif; sig=signew; j = j+1; endo; /* check if sigma = 0 : xhat = x */ if maxc(maxc(abs(sig))) < sigtol; "findvk: ERROR sigma looks like zero sigma matrix, but zero"; " solution did not work. setting to zero but check"; sig = 0; endif; endif; /* find VAR innovation var cov matrix */ sigvar = n'*term*n; /* check if sigma = 0 : xhat = x */ if sig == 0; term = c*c'; sigvar = n'*term*n; if det(sigvar) < sigvtol; Q = n'c; /* kludge. better diagonalization Q for singular case? */ if rows(q) > cols(q); q = q~zeros(rows(q),rows(q)-cols(q)); endif; if cols(q) > rows(q); if not maxc(maxc(dotfne(q[.,rows(q)+1:cols(q)],0))); q = q[.,1:rows(q)]; else; "findvk: can't prune q. find better diag. method"; endif; endif; endif; endif; /* find K */ k = term*nsm*inv(nsm'*term*nsm); /* choleski to make orthogonal shocks */ if det(sigvar) > sigvtol; Q = (chol(sigvar))'; /* gauss gives upper triangular. qq' = sigvar */ endif; if Q == 0; "findvk: failed to specify Q, will bomb"; endif; /* simulate impulse-responses */ j = 1; do while j <= cols(n); shock = zeros(cols(n),1); shock[j] = 1; z = (Q*shock)'; xhat = k*(z[1:cols(nsm)]'); t = 2; do while t <= irlags; z = z|((n'a*xhat)') ; xhat = a*xhat; /* why does the following not work? */ /* xhat = (eye(cols(nsm))-K*nsm')*A*xhat + K*(z[rows(z),1:cols(nsm)]'); /* only nsm are used to form xhat if singular */ z = z|((n'*A*xhat)'); */ t = t+1; endo; if j == 1; ir = z; else; ir = ir~z; endif; j = j+1; endo; retp(ir,sigvar); endp; /* test call*/ /*{ir,sigvar}= findvk(eye(2),(.9~0)|(0~.8),eye(2));*/ /* {ir,sigvar} = findvk(nc~ny,Am,Cm); ir[1:10,.]; */