end(reset)

OPEN DATA
CALENDAR(Q) 1976:1
DATA(FORMAT=free,NOLABELS,ORG=COLUMNS,TOP=2) 1976:01 2015:04 ctgdp	aux

* Get initial values - Standardised HP filter settings

@localdlm(type=trend,shocks=trend,a=ahp,c=chp,f=fhp)
compute lambda = 400000

dlm(a=ahp,c=chp,f=fhp,sv=1.0,sw=1.0/lambda,presample=diffuse,$
  type=filter,var=concentrate,y=ctgdp) / hpstates
set hp_ctgdp_t   = hpstates(t)(1)
set hp_ctgdp_c = ctgdp - hp_ctgdp_t

linreg hp_ctgdp_c
# hp_ctgdp_c{1 2}
compute ph1=%beta(1),ph2=%beta(2),se=sqrt(%seesq)



* UC-2 decomposition with AR(2) cycle, Random walk with discrete shock trend or $
*   markov-switching trend growth, auxilliary observable.

*********************************************************************
*
* This is common to any MS state space model analyzed using the
* Kim filter.
*
@MSSetup(states=2)
*
* These are the collapsed SSM mean and variance from the previous time
* period.
*
dec vect[vect] xstar(nstates)
dec vect[symm] sstar(nstates)
*
* This is the filtered best estimate for the state vector
*
dec series[vect] xstates
*
*********************************************************************

* This is specific to this model

* Number of lags in the AR

compute q=2
dec vect phi(q)

* Fill in the fixed coefficients in the A, C and F matrices
set dctgdp = ctgdp-ctgdp{1}
set daux = aux-aux{1}

compute ndlm=q

dec rect a_uc2(ndlm,ndlm)
ewise a_uc2(i,j)=(i==j+1)

linreg daux
# hp_ctgdp_c{0 1}
compute a3=%beta(1),a4=%beta(2),sigec2=sqrt(%seesq)

function %%DLMStart
compute [rect] c_uc2 =||1.0,a3|$
                       -1.0,a4||
compute [rect] sv_uc2=||0.0|$
                        sigec2^2||
end %%DLMStart

dec rect f_uc2(ndlm,1)
compute [rect] f_uc2=||1.0|$
                       0.0||

dec symm sw_uc2(ndlm,ndlm)


dec vect mu_uc2(nstates)

dec vect x0(ndlm)
compute x0=%zeros(ndlm,1)

nonlin(parmset=initparms) x0
nonlin(parmset=dlmparms) mu_uc2 sigsq sigec2 phi a3 a4
nonlin(parmset=msparms) theta
compute theta=%const(0.0)

*********************************************************************
*
* This does a single step of the Kim (approximate) filter
*
function KimFilter time
type integer 		time
*
local integer    	i j
local real       	yerr likely
local symm       	vhat
local rect       	gain
local rect       	phat(nstates,nstates)
local rect    	  	fwork(nstates,nstates)
local rect[vect] 	xwork(nstates,nstates)
local rect[symm] 	swork(nstates,nstates)
*
do i=1,nstates
   do j=1,nstates
      *
      * Do the SSM predictive step
      *
      compute xwork(i,j)=a_uc2*xstar(j)
      compute swork(i,j)=a_uc2*sstar(j)*tr(a_uc2)+sw_uc2
      *
      * Do the prediction error for y under state i, and compute the
      * density function for the prediction error.
      *
      compute yerr=ctgdp(time)-(%dot(c_uc2(i,j),xwork(i,j))+mu_uc2(i))
      compute vhat=c_uc2(i,j)*swork(i,j)*tr(c_uc2)+sv_uc2
      compute gain=swork(i,j)*tr(c_uc2)*inv(vhat)
      compute fwork(i,j)=exp(%logdensity(vhat,yerr))
      *
      * Do the SSM update step
      *
      compute xwork(i,j)=xwork(i,j)+gain*yerr
      compute swork(i,j)=swork(i,j)-gain*c_uc2(i,j)*swork(i,j)
   end do j
end do i
*
* Compute the Hamilton filter likelihood
*
compute likely=0.0
do i=1,nstates
   do j=1,nstates
      compute phat(i,j)=p(i,j)*pstar(j)*fwork(i,j)
      compute likely=likely+phat(i,j)
   end do j
end do i
*
* Compute the updated probabilities of the state combinations
*
compute phat=phat/likely
compute pstar=%sumr(phat)
compute pt_t(time)=pstar
*
* Collapse the SSM matrices down to one per state
*
compute xstates(time)=%zeros(ndlm,1)
do i=1,nstates
   compute xstar(i)=%zeros(ndlm,1)
   do j=1,nstates
      compute xstar(i)=xstar(i)+xwork(i,j)*phat(i,j)/pstar(i)
   end do j
   compute sstar(i)=%zeros(ndlm,ndlm)
   do j=1,nstates
      compute sstar(i)=sstar(i)+phat(i,j)/pstar(i)*$
         (swork(i,j)+%outerxx(xstar(i)-xwork(i,j)))
   end do j
   *
   * This is the overall best estimate of the filtered state.
   *
   compute xstates(time)=xstates(time)+xstar(i)*pstar(i)
end do i
compute KimFilter=likely
end
*********************************************************************
*
* This is the start up code for each function evaluation
*
function DLMStart
*
local integer i
*
* Fill in the top row in the A matrix
*
compute %psubmat(a_uc2,1,1,tr(phi))
*
* Compute the full matrix for the transition variance
*
compute sw_uc2=f_uc2*tr(f_uc2)*sigsq
*
* Transform the theta to transition probabilities.
*
compute p=%MSLogisticP(theta)
compute pstar=%mcergodic(p)
compute p=%mspexpand(p)
*
* Initialize the KF state and variance. This is set up for estimating
* the pre-sample states.
*
ewise xstar(i)=x0
ewise sstar(i)=%zeros(ndlm,ndlm)
end
*********************************************************************
*
* Get guess values for the means from running an AR(2) on the growth
* rate. However, the guess values for the AR coefficients need to be
* input.
*
boxjenk(ar=q,constant) dctgdp
compute mu_uc2(1)=%beta(1)+1.0,mu_uc2(2)=%beta(1)-1.0
compute phi=||ph1,ph2||
compute sigsq=%seesq
*
frml kimf = kf=KimFilter(t),log(kf)
compute x0=||5.224,2.699||
*********************************************************************
gset xstates = %zeros(ndlm,1)
@MSFilterInit
maximize(parmset=msparms+dlmparms,start=DLMStart(),method=bfgs,trace) kimf 1976:2 2015:4
maximize(parmset=msparms+dlmparms+initparms,start=DLMStart(),method=bfgs,trace) kimf 1976:2 2015:4
*
set cycle_uc2 1976:2 2015:4 = xstates(t)(1)
set trend_uc2 1976:2 2015:4 = ctgdp-cycle
*
set phigh 1976:2 2015:4 = pt_t(t)(1)
graph(footer="Figure 5.2a Filtered Probabilities of High-Growth State")
# phigh
graph(footer="Figure 5.3 PSC/GDP ratio Markov-switching trend component") 2
# ctgdp 1976:2 2015:4
# trend_uc2
graph(footer="Figure 5.4 Cyclical component")
# cycle_uc2