If all of 4 variables use current cycle,its codes as follows. however, the coincident indicator is negative and downwards, and is so different comparing with origin.
why? could you give me right code.
Code: Select all
*
* KIMNP126.RPF
* Application 3 from pp 126-133 of Kim and Nelson, "State-space Models
* with Regime Switching".
*
* Dynamic factor model with Markov switching component
*
open data s&w.xls
cal(m) 1959:1
data(format=xls,org=columns) 1959:01 1995:01 ip gmyxpq mtq lpnag dcoinc
*
set ipgrow = 100.0*log(ip/ip{1})
set incgrow = 100.0*log(gmyxpq/gmyxpq{1})
set salesgrow = 100.0*log(mtq/mtq{1})
set empgrow = 100.0*log(lpnag/lpnag{1})
*
* Parameters for index
*
dec vect phi(2)
compute phi=%zeros(2,1)
dec real sigmac
compute sigmac=1.0
*
* Parameters for the series. sigma is initialized based upon the sample
* means.
*
dec vect beta(4)
dec vect[vect] gamma(4)
dec vect sigma(4)
dec vect[vect] psi(4)
ewise psi(i)=%zeros(2,1)
*
stats(noprint) ipgrow
compute beta(1)=%mean,sigma(1)=%variance*.5
compute gamma(1)=%fill(1,1,0.6)
stats(noprint) incgrow
compute beta(2)=%mean,sigma(2)=%variance*.5
compute gamma(2)=%fill(1,1,0.5)
stats(noprint) salesgrow
compute beta(3)=%mean,sigma(3)=%variance*.5
compute gamma(3)=%fill(1,1,0.4)
stats(noprint) empgrow
compute beta(4)=%mean,sigma(4)=%variance*.5
compute gamma(4)=%fill(1,1,0.2)
*
* Construct fixed parts of the A matrix. We need four lags for the cycle
* (first block) because of the required lags for the loadings on the
* employment series.
*
dec rect ar2(2,2)
ewise ar2(i,j)=(i==j+1)
compute a=ar2~\ar2~\ar2~\ar2~\ar2
compute ndlm=%rows(a)
*
* These are the positions in the state vector of the noise terms for the
* series.
*
dec vect[integer] spos(4)
compute spos=||3,5,7,9||
*
* Construct the F matrix (always fixed). Again, we need a size 4 block
* for the cycle.
*
dec rect f2(2,1)
compute f2=%unitv(2,1)
compute f=f2~\f2~\f2~\f2~\f2
*
* Construct the fixed parts of the C matrix
*
compute c=%zeros(2,4)~~(f2~\f2~\f2~\f2)
*
dec symm sw(4,4)
*********************************************************************
*
* Setup for Kim filter
*
*********************************************************************
*
* 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)
*
* These are the work SSM mean and variance for the current time period
*
dec rect[vect] xwork(nstates,nstates)
dec rect[symm] swork(nstates,nstates)
*
dec rect fwork(nstates,nstates)
*
dec series[vect] xstates
*
* Transition probabilities (in logistic form)
*
compute theta=%zeros(nstates-1,nstates)
*********************************************************************
*
* Set up the switching component (the "z" shift in the transition
* equation).
*
dec vect[vect] zswitch(nstates)
dec vect mu(nstates)
compute zswitch(1)=%zeros(ndlm,1)
compute zswitch(2)=%zeros(ndlm,1)
*
compute mu(1)=-1.0,mu(2)=1.0
*********************************************************************
function DLMStart
*
* Insert the lag coefficients for the cycle
*
local symm sx0
*
compute %psubmat(a,1,1,tr(phi))
*
* Insert the lag coefficients for the noise terms into the A matrix and
* the loading coefficients into the C matrix.
*
do i=1,4
compute %psubmat(a,spos(i),spos(i),tr(psi(i)))
compute %psubmat(c,1,i,gamma(i))
end do i
compute sw=f*%diag(sigmac~~sigma)*tr(f)
*
compute zswitch(1)(1)=mu(1)
compute zswitch(2)(1)=mu(2)
*
* Kim filter initialization
*
* Transform the theta to transition probabilities.
*
compute p=%MSLogisticP(theta)
compute pstar=%mcergodic(p)
compute p=%mspexpand(p)
*
compute sx0=%psdinit(a,sw)
ewise xstar(i)=%solve(%identity(ndlm)-a,zswitch(i))
ewise sstar(i)=sx0
end DLMStart
*
* The observables are the growth rates minus their estimated means. The
* paper uses the sample means rather than estimating these, so we leave
* beta out of the parameter set.
*
diff(standardize) ipgrow / cip
diff(standardize) incgrow / cinc
diff(standardize) salesgrow / csales
diff(standardize) empgrow / cemp
*
* The observables are the growth rates minus their estimated means. The
* paper uses the sample means rather than estimating these, so we leave
* beta out of the parameter set.
*
dec frml[vect] yf
frml yf = ||ipgrow-beta(1),incgrow-beta(2),$
salesgrow-beta(3),empgrow-beta(4)||
*********************************************************************
*
* This does a single step of the Kim (approximate) filter
*
function KimFilter time
type integer time
*
local integer i j
local vect yerr
local real likely
local symm vhat
local rect gain
local rect phat(nstates,nstates)
local rect fwork(nstates,nstates)
*
* Any "switch" code will be in this loop
*
do i=1,nstates
do j=1,nstates
*
* Do the SSM predictive step. In this model, the shift is in the
* "Z" component.
*
compute xwork(i,j)=a*xstar(j)+zswitch(i)
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=yf(time)-tr(c)*xwork(i,j)
compute vhat=tr(c)*swork(i,j)*c
compute gain=swork(i,j)*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*tr(c)*swork(i,j)
end do j
end do i
*
* From here until the end of the function is the same for all models
*
* 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 xstates(time)=xstates(time)+xstar(i)
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
end do i
compute KimFilter=likely
end
*
* Estimate the parameters by (approximate) maximum likelihood. This has
* multiple modes, one of which has a substantially higher likelihood
* value (once corrected for integrating constants) than the one cited,
* but with very odd parameter values. This might indicate that the
* approximation has failed to work adequately.
*
compute mu(1)=-1.0,mu(2)=.25
nonlin(parmset=dlmparms) gamma psi sigma phi mu
nonlin(parmset=msparms) theta
compute theta(1,1)=log(.8/.2),theta(1,2)=log(.2/.8)
*
frml kimf = kf=KimFilter(t),log(kf)
*********************************************************************
gset xstates = %zeros(ndlm,1)
@MSFilterInit
maximize(parmset=dlmparms+msparms,start=DLMStart(),method=bfgs) kimf 1959:2 *
disp "Log likelihood without integrating constants" %logl-4.0*%nobs*-.5*log(2*%pi)
* other result: Coincident Index, Filtered Probabilities,transition Probabilities
set coincident = xstates(t)(1)
acc coincident
graph(footer="Coincident Index of 4 indicator of USA")
# coincident
set phigh %regstart() %regend() = pt_t(t)(2)
graph(footer="Filtered Probabilities of High-Growth Regime")
# phigh
computer pp = %mslogisticp(theta)
computer p11 = pp(1,1)
computer p22 = pp(2,2)
disp "p11= " p11 " , p22=" p22