DSGEHISTORY.RPF shows how to compute the historical decomposition for a DSGE model that's been fit to data. It is based upon the CAGAN.RPF example.


The Cagan model has two shocks: a money growth shock and a money demand shock. The historical decomposition computes the cumulated effects during the estimation period of the estimated shocks. The start of the program (through the estimation of the parametes) is identical to CAGAN.RPF. The following Kalman filters through the sample using the estimated model, including the WHAT option on the DLM to compute the filtered shocks.


dlm(start=%(EvalModel(),sw=%diag(||sig_eps^2,sig_eta^2||)),$

  a=adlm,f=fdlm,y=||x,mu||,c=cdlm,sw=sw,presample=ergodic,$

  type=filter,method=solve,what=what)

 

Because the model has (one) unit root, the WHAT vectors aren't defined for the first period, so this sets the range over which the shocks exist:

 

compute gstart=%regstart()+1,gend=%regend()


The Cagan model has the general form:


(1) \({{\bf{X}}_t} = {{\bf{A}}}{{\bf{X}}_{t - 1}} + {{\bf{F}}}{{\bf{W}}_t}\)


(2) \({{\bf{Y}}_t} = {{\bf{C}}}^\prime {{\bf{X}}_t}\)


where all the shocks are included in (1); the observation equation (2) is an identity. Because of linearity, the \({\bf{X}}_t \) can be decomposed into a sum of the evolution of the states from \({\bf{X}}_0 \) (with zero \({\bf{W}}_t \)), those with \({\bf{X}}_0  = 0\) and only the 1st components of \({\bf{W}}_t \), and those with \({\bf{X}}_0  = 0\) and only the 2nd components of \({\bf{W}}_t \); from that, \({\bf{Y}}_t \) can decomposed similarly into three parts by premultiplying those three by \({{\bf{C}}}^\prime\). The first of these (the effects of initial conditions) isn't as interesting, so we just compute how the two sets of shocks impact the evolution of \({\bf{Y}}_t \).


YHIST is what we will generate. This is a 2-vector on the outside (one per shock) of a SERIES[VECT] where each VECTOR is a 2-vector representing the dependent variables. The model is computed over the sample with a zeroed out initial X vector, and the WHAT(TIME) zeroed out except for the component being analyzed (by an elementwise product (.* operation) with a unit vector.

 

dec vect[series[vect]] yhist(2)

gset yhist(1) = %zeros(2,1)

gset yhist(2) = %zeros(2,1)

do i=1,2

   compute x0=%zeros(%rows(adlm),1)

   do time=gstart,gend

      compute x0=adlm*x0+fdlm*(what(time).*%unitv(%size(what(time)),i))

      *

      * Extract the part of the observables due to the state

      *

      compute yhist(i)(time)=tr(cdlm)*x0

   end do time

end do i


This does a graph of the cumulative effects of the shocks on the two observables:


spgraph(footer="Cumulative Effects of Shocks",vfields=2,$

    key=below,style=lines,klabels=||"Money Growth Shock","Money Demand Shock"||)

  set x_to_eps = yhist(1)(t)(1)

  set x_to_eta = yhist(2)(t)(1)

  graph(header="Inflation") 2

  # x_to_eps

  # x_to_eta

  set mu_to_eps = yhist(1)(t)(2)

  set mu_to_eta = yhist(2)(t)(2)

  graph(header="Money Growth") 2

  # mu_to_eps

  # mu_to_eta

spgraph(done)


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)
*************************************************************************
*
* Up to this point, this is the same as the CAGAN.RPF example. For this point
* on is the calculation of the historical decomposition.
*
* Compute the Kalman filter w shocks (into WHAT).
*
dlm(start=%(EvalModel(),sw=%diag(||sig_eps^2,sig_eta^2||)),$
  a=adlm,f=fdlm,y=||x,mu||,c=cdlm,sw=sw,presample=ergodic,$
  type=filter,method=solve,what=what)
*
* Because the model has (one) unit root, <<what>> is not defined for period 1
*
compute gstart=%regstart()+1,gend=%regend()
*
* Compute the incremental effect of the shocks, given the dynamics of the model.
* This starts with a zero presample value for the states and extrapolates that
* using the state equation, adding in only one set of historical shocks in each
* pass.
*
dec vect[series[vect]] yhist(2)
gset yhist(1) = %zeros(2,1)
gset yhist(2) = %zeros(2,1)
do i=1,2
   compute x0=%zeros(%rows(adlm),1)
   do time=gstart,gend
      compute x0=adlm*x0+fdlm*(what(time).*%unitv(%size(what(time)),i))
      *
      * Extract the part of the observables due to the state
      *
      compute yhist(i)(time)=tr(cdlm)*x0
   end do time
end do i
*
spgraph(footer="Cumulative Effects of Shocks",vfields=2,$
    key=below,style=lines,klabels=||"Money Growth Shock","Money Demand Shock"||)
  set x_to_eps = yhist(1)(t)(1)
  set x_to_eta = yhist(2)(t)(1)
  graph(header="Inflation") 2
  # x_to_eps
  # x_to_eta
  set mu_to_eps = yhist(1)(t)(2)
  set mu_to_eta = yhist(2)(t)(2)
  graph(header="Money Growth") 2
  # mu_to_eps
  # mu_to_eta
spgraph(done)

Output


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.410835100  2.523460742     -1.74793  0.08047599

2.  LAMBDA                        0.660060303  0.054608059     12.08723  0.00000000

3.  SIG_ETA                       0.219441069  0.175680074      1.24909  0.21163040

4.  SIG_EPS                       0.073583790  0.010355457      7.10580  0.00000000


Graph