This is an example of estimation of the parameters underlying a small DSGE model (Cagan's model of hyperinflation) It is from the User's Guide, Section 10.8.5.
The series MU and X are actual data. A1, A2, EPS and ETA are SERIES used in solving the model. This sets up the equations of the model:
declare series a1 a2 eps eta
declare real alpha lambda sig_eta sig_eps
*
frml(identity) f1 x = x{1}+a1-lambda*a1{1}
frml(identity) f2 mu = (1-lambda)*x{1}+lambda*mu{1}+a2-lambda*a2{1}
frml(identity) f3 a1 = 1.0/(lambda+(1-lambda)*alpha)*(eps-eta)
frml(identity) f4 a2 = (1.0/(lambda+(1-lambda)*alpha)*$
((1+alpha*(1-lambda))*eps-(1-lambda)*eta))
frml d1 = eps
frml d2 = eta
*
group cagan f1 f2 f3 f4 d1 d2
ALPHA, LAMBDA, SIG_ETA and SIG_EPS are the underlying parameters which will need to be estimated. While it's possible to eliminate A1 and A2 from this with some (tedious) algebra, there is no reason to do so. The guess values need to be set carefully because the F3 and F4 equations both divide by a function of the parameters, so we have to avoid zero values for that, and, in fact, need to stay on the correct (here, negative) side of it.
The following FUNCTION solves the model given values for the parameters, creating the state-space matrices ADLM and FDLM. (The model is entirely mean zero, so there is no need for a "Z").
function EvalModel
dsge(model=cagan,a=adlm,f=fdlm) x mu a1 a2 eps eta
end EvalModel
Because DSGE is a bit of a "black-box", you won't generally know how many states will be in the solved dynamic system. You do know where the endogenous variables themselves are, since you control that by the position on the DSGE instruction (thus X will be state 1, MU will be state 2, etc.). This does an evaluation of the model, and then creates a C matrix for the state-space model which pulls out just those first two states.
compute EvalModel()
compute cdlm=%identity(2)~~%zeros(%rows(adlm)-2,2)
Because we really have no idea what the scale is on the variances, we first estimate the model (using DLM) with the alpha and lambda fixed. This uses only a small number of simplex iterations to get the sigmas into the right zone. The START option first does a solution of the DSGE model (using the EvalModel function above) and sets the (diagonal) covariance matrix of the shocks—both of those depend only upon the parameters and are time-invariant. PRESAMPLE=ERGODIC is important here because this has a mix of unit (in F1) and non-unit roots.
nonlin sig_eta sig_eps
dlm(start=%(EvalModel(),sw=%diag(||sig_eps^2,sig_eta^2||)),$
a=adlm,f=fdlm,y=||x,mu||,c=cdlm,sw=sw,presample=ergodic,$
method=simplex,iters=5,noprint)
The final estimates are done with the full parameter set:
nonlin alpha lambda sig_eta sig_eps
dlm(start=%(EvalModel(),sw=%diag(||sig_eps^2,sig_eta^2||)),$
a=adlm,f=fdlm,y=||x,mu||,c=cdlm,sw=sw,presample=ergodic,$
pmethod=simplex,piters=5,method=bfgs)
Full Program
open data cagan_data.prn
data(format=prn,org=cols) 1 34 mu x
*
declare series a1 a2 eps eta
declare real alpha lambda sig_eta sig_eps
*
frml(identity) f1 x = x{1}+a1-lambda*a1{1}
frml(identity) f2 mu = (1-lambda)*x{1}+lambda*mu{1}+a2-lambda*a2{1}
frml(identity) f3 a1 = 1.0/(lambda+(1-lambda)*alpha)*(eps-eta)
frml(identity) f4 a2 = (1.0/(lambda+(1-lambda)*alpha)*$
((1+alpha*(1-lambda))*eps-(1-lambda)*eta))
frml d1 = eps
frml d2 = eta
*
group cagan f1 f2 f3 f4 d1 d2
*
* lambda+(1-lambda)*alpha needs to stay clear of its zero point. It
* appears that it needs to be negative, so alpha must be less than
* -lambda/(1-lambda) on the guess values.
*
compute alpha=-3.00,lambda=.7,sig_eta=.001,sig_eps=.001
*****
function EvalModel
dsge(model=cagan,a=adlm,f=fdlm) x mu a1 a2 eps eta
end EvalModel
*****
compute EvalModel()
compute cdlm=%identity(2)~~%zeros(%rows(adlm)-2,2)
*
* Because we really have no idea what the scale is on the variances, we
* first estimate the model with the alpha and lambda fixed. This uses
* only a small number of simplex iterations to get the sigmas into the
* right zone.
*
nonlin sig_eta sig_eps
dlm(start=%(EvalModel(),sw=%diag(||sig_eps^2,sig_eta^2||)),$
a=adlm,f=fdlm,y=||x,mu||,c=cdlm,sw=sw,presample=ergodic,$
method=simplex,iters=5,noprint)
*
* Now re-estimate freeing up the alpha and lambda
*
nonlin alpha lambda sig_eta sig_eps
dlm(start=%(EvalModel(),sw=%diag(||sig_eps^2,sig_eta^2||)),$
a=adlm,f=fdlm,y=||x,mu||,c=cdlm,sw=sw,presample=ergodic,$
pmethod=simplex,piters=5,method=bfgs)
Output
Note that the rank of the observables is 67 = 2*34-1. There are two observables per data point, but there is one unit root.
DLM - Estimation by BFGS
Convergence in 14 Iterations. Final criterion was 0.0000014 <= 0.0000100
Usable Observations 34
Rank of Observables 67
Log Likelihood 34.1493
Variable Coeff Std Error T-Stat Signif
************************************************************************************
1. ALPHA -4.410835089 2.523452045 -1.74794 0.08047494
2. LAMBDA 0.660060303 0.054607913 12.08726 0.00000000
3. SIG_ETA 0.219441069 0.175679627 1.24910 0.21162924
4. SIG_EPS 0.073583790 0.010355382 7.10585 0.00000000