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



* 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=dlmparms) mu sigsq phi tau
*
* Time varying probabilities
*
dec equation p1eq p2eq
dec vector   v1   v2
nonlin(parmset=tv) v1 v2
*********************************************************************
*
* This overrides the standard fixed definition
*
function  %MSPmat  time
type rect %MSPmat
type int  time
*
local integer i j
local rect    pexpand
local real    z
*
compute p(1,1)=%(z=exp(%dot(%eqnxvector(p1eq,time),v1)),z/(1+z))
compute p(1,2)=%(z=exp(%dot(%eqnxvector(p2eq,time),v2)),1/(1+z))
compute %MSInitTransition()
*
if nexpand==nstates {
   dim pexpand(nstates,nstates-1)
   ewise pexpand(i,j)=%if(i==nstates,-1.0,i==j)
   compute %MSPmat=pexpand*p
   ewise %MSPmat(i,j)=%if(i==nstates,1.0+%MSPmat(i,j),%MSPmat(i,j))
}
else {
   dim %MSPmat(nexpand,nexpand)
   ewise %MSPmat(i,j)=MSTransProbs(MSTransLookup(i,j))
}
end
*********************************************************************
*
* Initialize the probabilities using only the intercepts
*
function TVPInit
type vector TVPInit
*
local rect a
local integer i j
*
compute pstar=||%logistic(v1(1),1.0),1-%logistic(v2(1),1.0)||
compute TVPInit=%MSInit()
end

*********************************************************************
*
* 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 p=%MSPMat(time)

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 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||
*

* Define the logistic index for the transitions
*
equation p1eq *
# demo rep
equation p2eq *
# demo rep
*
* Standard guess values aiming to make state 1 the low growth state.
*
compute v1=log(.7/.3)~~0.0
compute v2=log(.95/.05)~~0.0




frml kimf = kf=KimFilter(t),log(kf)
*********************************************************************
gset xstates = %zeros(ndlm,1)
@MSFilterInit
maximize(parmset=tv+dlmparms,start=DLMStart(),method=bfgs,trace) kimf 1945:04 *