*
* 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