Dynamic Factor Model with Regime Switching

[TomDoan]

Dear Tom:


my questions continue for next step according you advice.
question 2:
If all of 4 variables use current cycle, then all would have %fill(1,...) in the settings for gamma. I just modify this one, the filter probability of hold similiar to origin setting, however, the coincident indicator is reverse from upwards to downwards.

The model does nothing to fix the sign of the cycle---if you multiply the cycle states by -1 and the gammas by -1, you get the identical model. The only thing aimed at enforcing positive loadings on gamma is the guess values. You might try adjusting them higher.

If all of 4 varibles use current on the cycle, some metric should be sparse, could some row and cols could remove. If so , A metric should be;

dec rect ar2(2,2)
ewise ar2(i,j)=(i==j+1)
compute a=ar~\ar2~\ar2~\ar2~\ar2

F metric should be:

dec rect f2(2,1)
compute f2=%unitv(2,1)
compute f=f2~\f2~\f2~\f2~\f2

C metric should be:
compute c=%zeros(2,4)~~(f2~\f2~\f2~\f2)

Is it right? the question still holds.

Yes. The extra two augmenting states on the cycle aren't necessary if you aren't using extra lags on the loadings.

question3: I use three varibles and current cycle, result is similar to that four varibles and current cycle. The question still holds.

I'm not sure what question you mean, but if you're talking about reducing the number of lag states on the cycle, yes.

[hardmann]

Dear Tom:

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 

[TomDoan]

What you have is correct. There is a reason they used extra lags on employment---the model doesn't really fit well with only a current loading of the cycle since employment lags the cycle.

Also, if the model works correctly, there actually should be no trend in the accumulated "cycle". The data are de-meaned growth rates, and so represent data that are already detrended. What happens when you exclude the extra lags on employment is you get a negative mean in the cycle state and positive means in all the noise terms, while you get the reverse of that when you include the extra lags on employment. As the model is designed, the noises should be zero mean processes, which would force the cycle to be zero mean as well.

[hardmann]

Dear Tom:
In kimnp126, how can forecast coincident index, and filtered probability?


Best regard
hardmann

[TomDoan]

Dear Tom:
In kimnp126, how can forecast coincident index, and filtered probability?

The probabilities are forecast exactly the same way as for any other Markov switching model. Kim and Nelson don't discuss the use of the model for forecasting. If you want to do that, you'll have to figure out what calculations are required. Forecasting the state in a standard state space model is an extension of the Kalman filter out of sample. Since the MS state space model is only approximated using the Kim filter, it might not be a simple matter to extend that out of sample.

[hardmann]

Dear Tom:

I hav known the method of forecast with state space, however, I had not know and found method and example with Markov swithing from user guide and E_course.
Could you tell me a example or procedure?


Best regard
Hardmann

[TomDoan]

Dear Tom:

I hav known the method of forecast with state space, however, I had not know and found method and example with Markov swithing from user guide and E_course.
Could you tell me a example or procedure?

I guess you've noticed that the literature on forecasting with Markov switching models is extremely thin. You might want to think about why that is.

[rawraw]

acc indicator.RGF

acc indicator.RGF

Dear Tom:

I hav known the method of forecast with state space, however, I had not know and found method and example with Markov swithing from user guide and E_course.
Could you tell me a example or procedure?

I guess you've noticed that the literature on forecasting with Markov switching models is extremely thin. You might want to think about why that is.

Hi Tom

My model contains three variables , in estimation result, why are all psi 0,also the t statistics and significance? the indicator contains negtive value,but the original variables in my model are positive,what is the problem? estimation results are as follows .thanks a lot.
MAXIMIZE - Estimation by BFGS
Convergence in 51 Iterations. Final criterion was 0.0000056 <= 0.0000100
Monthly Data From 1998:02 To 2013:08
Usable Observations 187
Function Value -1165.4565

Variable Coeff Std Error T-Stat Signif
************************************************************************************
1. GAMMA(1)(1) 0.28265467 0.08905458 3.17395 0.00150380
2. GAMMA(2)(1) 0.01802997 0.11668133 0.15452 0.87719724
3. GAMMA(3)(1) 0.06906504 0.00777962 8.87769 0.00000000
4. GAMMA(3)(2) 0.06926082 0.00743392 9.31687 0.00000000
5. PSI(1)(1) 0.00000000 0.00000000 0.00000 0.00000000
6. PSI(1)(2) 0.00000000 0.00000000 0.00000 0.00000000
7. PSI(2)(1) 0.00000000 0.00000000 0.00000 0.00000000
8. PSI(2)(2) 0.00000000 0.00000000 0.00000 0.00000000
9. PSI(3)(1) 0.00000000 0.00000000 0.00000 0.00000000
10. PSI(3)(2) 0.00000000 0.00000000 0.00000 0.00000000
11. SIGMA(1) 62.82274663 6.10447175 10.29127 0.00000000
12. SIGMA(2) 98.42792375 9.45054056 10.41506 0.00000000
13. SIGMA(3) 0.00002119 0.00000883 2.39850 0.01646257
14. PHI(1) 1.58589100 0.03887149 40.79831 0.00000000
15. PHI(2) -0.84189232 0.03912040 -21.52055 0.00000000
16. MU(1) -0.87504028 0.17232829 -5.07775 0.00000038
17. MU(2) 0.95409048 0.17176147 5.55474 0.00000003
18. THETA(1,1) 2.35596368 0.45165980 5.21624 0.00000018
19. THETA(1,2) -2.23478259 0.41491415 -5.38613 0.00000007

Log likelihood without integrating constants -478.09046
0.91341 0.09667
0.08659 0.90333

[TomDoan]

That looks like the PSI's never enter the model so you must have set it up wrong. You'll have to post the program, not the output.

[rawraw]

That looks like the PSI's never enter the model so you must have set it up wrong. You'll have to post the program, not the output.

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 levle-three.rat
cal(m) 1998:1
data(format=rat,org=columns) 1998:01 2013:08 cd	ra sc bb
*
set cdgrow    = 100.0*log(cd/cd{1})
set ragrow   = 100.0*log(ra/ra{1})
set scgrow = 100.0*log(sc/sc{1})
set bbgrow = 100.0*log(bb/bb{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(3)
dec vect[vect] gamma(3)
dec vect sigma(3)
dec vect[vect] psi(3)
ewise psi(i)=%zeros(2,1)
*
stats(noprint) cdgrow
compute beta(1)=%mean,sigma(1)=%variance*.5
compute gamma(1)=%fill(1,1,0.6)
stats(noprint) scgrow
compute beta(2)=%mean,sigma(2)=%variance*.5
compute gamma(2)=%fill(1,1,0.5)
stats(noprint) ragrow
compute beta(3)=%mean,sigma(3)=%variance*.5
compute gamma(3)=%fill(2,1,0.4)

*
* 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 ar4(3,3)
dec rect ar2(2,2)
ewise ar4(i,j)=(i==j+1)
ewise ar2(i,j)=(i==j+1)
compute a=ar4~\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(3)
compute spos=||3,5,7||
*
* Construct the F matrix (always fixed). Again, we need a size 4 block
* for the cycle.
*
dec rect f4(3,1)
dec rect f2(2,1)
compute f4=%unitv(3,1)
compute f2=%unitv(2,1)
compute f=f4~\f2~\f2~\f2
*
* Construct the fixed parts of the C matrix
*
compute c=%zeros(3,3)~~(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,3
   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.
*
dec frml[vect] yf
frml yf = ||cdgrow-beta(1),scgrow-beta(2),ragrow-beta(3)||
*********************************************************************
*
* 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(.9/.1),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 1998:2 *
disp "Log likelihood without integrating constants" %logl-4.0*%nobs*-.5*log(2*%pi)

set indicator 1998:02 * = xstates(t)(1)
acc indicator
graph(footer="Coincident Indicator, Smoothed Version")
# indicator

[TomDoan]

With the rest of your set up, you need:

compute spos=||3,5,7|

to be 4,6,8. That looks like the only problem, but you might want to check out the A, C and F matrices that get generated to make sure they make sense.

[rawraw]

With the rest of your set up, you need:

compute spos=||3,5,7|

to be 4,6,8. That looks like the only problem, but you might want to check out the A, C and F matrices that get generated to make sure they make sense.

why 4,6,8? according to gamma(3)=%fill(2,1,0.4) ,spos should be 3,5,7,thanks!

[TomDoan]

You set up the A matrix like this:

dec rect ar4(3,3)
dec rect ar2(2,2)
ewise ar4(i,j)=(i==j+1)
ewise ar2(i,j)=(i==j+1)
compute a=ar4~\ar2~\ar2~\ar2

which means that the common block takes states 1-3, so the others start at 4, 6 and 8.

[rawraw]

You set up the A matrix like this:

dec rect ar4(3,3)
dec rect ar2(2,2)
ewise ar4(i,j)=(i==j+1)
ewise ar2(i,j)=(i==j+1)
compute a=ar4~\ar2~\ar2~\ar2

which means that the common block takes states 1-3, so the others start at 4, 6 and 8.

would you please explain the meaning of the code in every line as follows,many thanks!
*
* 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 ar4(3,3)
dec rect ar2(2,2)
ewise ar4(i,j)=(i==j+1)
ewise ar2(i,j)=(i==j+1)
compute a=ar4~\ar2~\ar2~\ar2
compute ndlm=%rows(a)
*
* Construct the F matrix (always fixed). Again, we need a size 4 block
* for the cycle.
*
dec rect f4(3,1)
dec rect f2(2,1)
compute f4=%unitv(3,1)
compute f2=%unitv(2,1)
compute f=f4~\f2~\f2~\f2
*
* Construct the fixed parts of the C matrix
*
compute c=%zeros(3,3)~~(f2~\f2~\f2)
*
dec symm sw(4,4

[TomDoan]

As was suggested above, the best way to understand this model with regime switching is to start with the same model without it.