*
* SS_9_3
* State-space and DSGE, 2nd edition
* Non-linear Kalman filter
*
* Based upon Matheson and Stavrev(2013), "The Great Recession and
* the inflation puzzle," Economics Letters, vol. 120, no 3, pp
* 468-472 with a reconstructed (and extended) data set.
*
* This does the extended Kalman filter but doesn't impose
* constraints on the time-varying coefficients.
*
******************
*
* Get styles to do thin and fat colored lines for the estimates and
* bounds.
*
open styles mscolors.txt
grparm(import=styles)
*
cal(q) 1960
*
* This is a reconstructed data set extended to 2016:4. (Note:
* PTRPCE is the long-term expected PCE inflation which is not used).
*
open data msextended.rat
data(format=rats) 1960:01 2016:04 ptrcpi ptrpce cpi unrate mdef
*
set u      = unrate
set pi     = 400.0*log(cpi/cpi{1})
set pim    = 400.0*log(mdef/mdef{1})-pi
set pibar  = ptrcpi
set pistar = .25*(pi+pi{1}+pi{2}+pi{3})
*
@nbercycles(down=nber)
*
* States are u-u* (ugap), u*, kappa, theta, gamma in that order.
* The observables are pi and u.
*
* These are the expansion points for the measurement equation for
* pi.
*
dec series kappa_0 ugap_0
dec series theta_0 gamma_0
*
* Initialize the linearization with an HP filtered unemployment for ustar
*
filter(type=hp) u / ustar_0
*
set ugap_0 = u-ustar_0
*
* Do rolling regressions
*
nonlin(parmset=rollparms) kappa_r theta_r gamma_r
frml rollpi pi = pistar{1}+theta_r*(pibar-pistar{1})-$
       kappa_r*ugap_0+gamma_r*pim
*
* Span for rolling regressions, first entry available for
* estimation (we lose four data points in computing pistar, then
* use a lag of pistar as well), first end point for which rolling
* regressions can be done and last entry available for estimation.
*
compute span   =40
compute sstart =1961:02
compute rstart =sstart+span-1
compute send   =2016:04
*
* This does the full sample regression
*
compute kappa_r=0.1,gamma_r=0.1,theta_r=0.5
nlls(parmset=rollparms,frml=rollpi) pi sstart send
compute sv=%sigmasq
*
* This does rolling sample regressions and saves the series of
* estimated coefficients.
*
clear(zeros) kappa_re theta_re gamma_re
do time=rstart,send
   compute kappa_r=0.1,gamma_r=0.1,theta_r=0.5
   nlls(parmset=rollparms,frml=rollpi,noprint) pi time-(span-1) time
   compute kappa_re(time)=kappa_r
   compute theta_re(time)=theta_r
   compute gamma_re(time)=gamma_r
end do time
*
set kappa_0 = %if(t<=rstart,kappa_re(rstart),kappa_re)
set theta_0 = %if(t<=rstart,theta_re(rstart),theta_re)
set gamma_0 = %if(t<=rstart,gamma_re(rstart),gamma_re)
*
* Compute variances of changes in the coefficients to use as upper
* bound on coefficient shock variance.
*
sstats(mean) rstart+1 send (kappa_0-kappa_0{1})^2>>sigsqup_kappa $
                           (theta_0-theta_0{1})^2>>sigsqup_theta $
                           (gamma_0-gamma_0{1})^2>>sigsqup_gamma
*
linreg ugap_0
# ugap_0{1}
compute rho=%beta(1)
compute swugap=%sigmasq
*
dec frml[vect] muf
dec frml[rect] cf
dec frml[symm] svf
dec frml[symm] swf
*
frml muf = ||+kappa_0*ugap_0+pistar{1},0.0||
*
frml cf  = ||-kappa_0   ,1.0|$
             0.0        ,1.0|$
             -ugap_0    ,0.0|$
       (pibar-pistar{1}),0.0|$
              pim       ,0.0||
*
dec vector swtvc(4)
*
frml svf = %diag(||sv,0.0||)
frml swf = %diag(swugap~~swtvc)
*
nonlin(parmset=base) rho sv swugap swtvc
*
* These are the two variance ratio pegs
*
nonlin(parmset=stable)   swugap=15.0*swtvc(1)
nonlin(parmset=flexible) swugap=5.0*swtvc(1)
*
* A is an identity matrix other than the 1,1 element. %%DLMSetup will
* reset that to the current test value of rho.
*
dec rect a(5,5)
compute a=%na~\%identity(4)
*
function %%DLMSetup
compute a(1,1)=rho
end
*
nonlin(parmset=limits) swtvc(1)>=0.0 rho<=.95 $
                       swtvc(2)>=0.0 swtvc(2)<=sigsqup_kappa $
                       swtvc(3)>=0.0 swtvc(3)<=sigsqup_theta $
                       swtvc(4)>=0.0 swtvc(4)<=sigsqup_gamma
compute swtvc(1)=.1
compute swtvc(2)=.5*sigsqup_kappa
compute swtvc(3)=.5*sigsqup_theta
compute swtvc(4)=.5*sigsqup_gamma
*
* Do one Kalman smooth pass to get expansion points for future
* non-linear approximations.
*
dlm(start=%%DLMSetup(),a=a,c=cf,y=||pi,u||,mu=muf,sv=svf,sw=swf,$
  presample=ergodic,type=smooth) sstart send xstates vstates
set ugap_0  = xstates(t)(1)
set kappa_0 = xstates(t)(3)
*
dec hash[parmset] pegs
compute pegs("stab")=stable
compute pegs("flex")=flexible
*
dec hash[hash[hash[series]]] h_ustar h_kappa h_theta h_gamma
dec hash[hash[series]] h_pie
dec hash[series[vect]] h_xstates
dec hash[series[symm]] h_vstates
*
dec string root ktype
dofor root = "stab" "flex"
   dlm(start=%%DLMSetup(),a=a,c=cf,y=||pi,u||,mu=muf,sv=svf,sw=swf,type=filter,$
      parmset=base+pegs(root)+limits,presample=ergodic,method=bfgs,condition=10) $
         sstart send h_xstates("kf") h_vstates("kf")
   *
   dlm(start=%%DLMSetup(),a=a,c=cf,y=||pi,u||,mu=muf,sv=svf,sw=swf,type=smooth,$
      presample=ergodic) $
         sstart send h_xstates("sm") h_vstates("sm")
   *
   dofor ktype = "kf" "sm"
      *
      * NAIRU estimates
      *
      set h_ustar(root)(ktype)("est") = h_xstates(ktype)(t)(2)
      set h_ustar(root)(ktype)("hi")  = h_xstates(ktype)(t)(2)+sqrt(h_vstates(ktype)(t)(2,2))
      set h_ustar(root)(ktype)("lo")  = h_xstates(ktype)(t)(2)-sqrt(h_vstates(ktype)(t)(2,2))
      *
      * Kappa
      *
      set h_kappa(root)(ktype)("est") = h_xstates(ktype)(t)(3)
      set h_kappa(root)(ktype)("hi")  = h_xstates(ktype)(t)(3)+sqrt(h_vstates(ktype)(t)(3,3))
      set h_kappa(root)(ktype)("lo")  = h_xstates(ktype)(t)(3)-sqrt(h_vstates(ktype)(t)(3,3))
      *
      * Theta
      *
      set h_theta(root)(ktype)("est") = h_xstates(ktype)(t)(4)
      set h_theta(root)(ktype)("hi")  = h_xstates(ktype)(t)(4)+sqrt(h_vstates(ktype)(t)(4,4))
      set h_theta(root)(ktype)("lo")  = h_xstates(ktype)(t)(4)-sqrt(h_vstates(ktype)(t)(4,4))
      *
      * Gammas
      *
      set h_gamma(root)(ktype)("est") = h_xstates(ktype)(t)(5)
      set h_gamma(root)(ktype)("hi")  = h_xstates(ktype)(t)(5)+sqrt(h_vstates(ktype)(t)(5,5))
      set h_gamma(root)(ktype)("lo")  = h_xstates(ktype)(t)(5)-sqrt(h_vstates(ktype)(t)(5,5))
      *
      * Inflation expectations (uses time-varying weights on observable data)
      *
      set h_pie(root)(ktype) = h_theta(root)(ktype)("est")*pibar+(1-h_theta(root)(ktype)("est"))*pistar{1}
   end do ktype
end do root
*
compute title="TVC Unconstrained"
source msgraphs.src