OPEN DATA "C:\Users\Sam\Desktop\quarterlydata.xlsx"
CALENDAR(Q) 1920:1
DATA(FORMAT=XLSX,ORG=COLUMNS,SHEET="q2") 1920:01 2015:02 quarter indgrowth rep yr34 gdp gnp indgrowths
set demo = 1 -rep
set yr12 = 1- yr34
@nbercycles(down=contract)


* 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
*
*********************************************************************
*
*
* Fill in the fixed coefficients in the A, C and F matrices
*
compute ndlm=5
dec rect a(ndlm,ndlm)
input a
. . 0 0 0
1 0 0 0 0
0 1 0 0 0
0 0 1 0 0
0 0 0 1 0

dec vect phi(2)

dec rect c(1,ndlm)
input c
1 . . . .

dec vect tau(4)

dec rect f(ndlm,1)
compute f=%unitv(ndlm,1)
dec symm sw(ndlm,ndlm)
*
dec vect mu(nstates)
compute sv=0.0
*
dec vect x0(ndlm)
compute x0=%zeros(ndlm,1)
*
nonlin(parmset=initparms) x0
nonlin(parmset=dlmparms) mu sigsq phi tau
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*xstar(j)
      compute swork(i,j)=a*sstar(j)*tr(a)+sw
      *
      * Do the prediction error for y under state i, and compute the
      * density function for the prediction error.
      *
      compute yerr=indgrowth(time)-(%dot(c,xwork(i,j))+mu(i))
      compute vhat=c*swork(i,j)*tr(c)+sv
      compute gain=swork(i,j)*tr(c)*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*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,1,1,tr(phi))
*
* Fill in the C matrix
*
compute %psubmat(c,1,2,tr(tau))
*
* Compute the full matrix for the transition variance
*
compute sw=f*tr(f)*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=2,constant) indgrowth
compute mu(1)=%beta(1)+1.0,mu(2)=%beta(1)-1.0
compute phi=||1.2,-.3||
compute sigsq=%seesq
compute tau=||0.1,0.1,0.1,-0.4||

*
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 1945:04 *
maximize(parmset=msparms+dlmparms+initparms,start=DLMStart(),method=bfgs,trace) kimf  1945:04 *