* * DLMCYCLE.RPF * Estimation of a state-space model with a common growth cycle * * RATS User's Guide, Example from Section 10.7. * cal(q) 1947 open data coninc.rat data(format=rats) 1947:1 1989:3 gyd82 gc82 * * Convert to growth rates * set ygr = 100.0*log(gyd82/gyd82{1}) set cgr = 100.0*log(gc82/gc82{1}) * dec rect g(1,2) nonlin phic phi1 phi2 mu1 mu2 g sigc=1.0 sig11 sig22 * * The DLM instruction's "A" matrix (H in the text) will be diagonal with * the phi's down the diagonal. The SW matrix will similarly be diagonal, * with the three sigma components. The "C" matrix (A in the text) is a * vector of the g parameters in the first row (remember, it's input in * transposed form) and an identity matrix below that. * * These are all time-invariant, but depend upon the free parameters, so * we set up a function to use with the START option on DLM; that will * copy the current settings into fixed matrices for use on the * instruction. * function %%DLMStart compute a=%diag(||phic,phi1,phi2||) compute sw=%diag(||sigc,sig11,sig22||) compute c=g~~%identity(2) end %%DLMStart * * Estimate an AR(1) on the sum of the two series. Use that to get a * guess for the PHIC. * set combined = ygr+cgr linreg combined # constant combined{1} compute phic=%beta(2) * * Get a preliminary estimate of the cycle by taking the fitted values and * rescaling them to have a unit residual variance. * prj testc set testc = testc/sqrt(%seesq) * * Estimate AR(1) regressions on the two observables and pull out guess * values. * ar1 ygr # constant testc compute mu1=%beta(1),g(1,1)=%beta(2),phi1=%rho,sig11=%seesq ar1 cgr # constant testc compute mu2=%beta(1),g(1,2)=%beta(2),phi2=%rho,sig22=%seesq * * Estimate the parameters by maximum likelihood. Given the estimated * model, use Kalman smoothing to estimate the business cycle. * dlm(start=%%DLMStart(),type=smooth,y=||ygr,cgr||,sw=sw,c=c,a=a,\$ mu=||mu1,mu2||,presample=ergodic,pmethod=simplex,piters=10,\$ method=bfgs) 1947:2 1989:3 xstates * * The cycle series is the first component of the state vector. Note that * the model doesn't make the cycle process sufficiently "stiff". * Something with more persistence than an AR(1) would probably give * results that are more in line with our notion of a business cycle. * set cycle 1947:2 1989:3 = xstates(t)(1) graph(footer="Estimated State of Business Cycle") # cycle