TVTP Markov Regime switching model

[moderator]

I'd suggest that you email the code and data set to us at

support@estima.com

That will be easier than trying to cut and paste code out of the forum post.

Please include your full name and RATS serial number in the email.

Regards,
Tom Maycock
Estima

[TomDoan]

This will estimate the model using logistic index transitions using EM with a BHHH polishing.

Code: Select all

cal(d) 2006:03:24
open data "one_sample_ap.xls"
data(format=xls,org=cols) 2006:03:24 2012:03:05 apexcess aprm aprmsq logvix
sample(smpl=%valid(apexcess)) apexcess / apexcess_nomissing
sample(smpl=%valid(aprm)) aprm / aprm_nomissing
sample(smpl=%valid(aprmsq)) aprmsq / aprmsq_nomissing
sample(smpl=%valid(logvix)) logvix / logvix_nomissing

calendar(irregular)
compute gstart=1,gend=%nobs
@msregression(states=2,switch=ch) apexcess_nomissing
# constant aprm_nomissing aprmsq_nomissing
@msreginitial gstart gend
@MSFilterSetup(em)
*
nonlin(parmset=common) betas sigsqv
dec equation p1eq   p2eq
dec vector   v1     v2
dec vector   p1mean p2mean
nonlin(parmset=tv) v1 v2
*
* Define the logistic index for the transitions
*
equation p1eq *
# constant
equation p2eq *
# constant
dim p1mean(1) p2mean(1)
compute p1mean(1)=1.0,p2mean(1)=1.0
**************************************************************************
*
* Function to compute the transition probability matrix
*
function %MSPmat time
type rect %MSPmat
type int time

local integer i j
local rect pexpand
local real z
local rect p
dim p(1,2)
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 means
*
function TVPInit
type vector TVPInit
*
local rect a
local integer i j
*
dim p(1,2)
compute p(1,1)=%(z=exp(%dot(p1mean,v1)),z/(1+z))
compute p(1,2)=%(z=exp(%dot(p2mean,v2)),1/(1+z))
compute %MSInitTransition()
compute TVPInit=%MSInit()
compute p=%MSPExpand(p)
end
**************************************************************************
procedure MSTVPEMEstimate start end
type integer start end
option integer iters 100
option real cvcrit 1.e-6

local integer emits time k
local series l11 l22 w11 w22
local vector baseparms testparms
local rect thisEntry
*
dim thisEntry(nstates,nstates)
*
compute baseparms=v1~~v2
do emits=1,iters
   compute %eqnsetcoeffs(p1eq,v1)
   compute %eqnsetcoeffs(p2eq,v2)
   compute TVPInit()
   @MSRegestep start end
   compute %emlogl=%logl
   @MSRegmstep(nodop,nodobeta,dosigma) start end
   set l11 start end = 0.0
   set l22 start end = 0.0
   set w11 start end = 0.0
   set w22 start end = 0.0

   compute tvemlogl=0.0
   do time=start,end
      compute p=%MSPMat(time)
      compute thisEntry=%zeros(nstates,nstates)
      do k=1,MSEMSize
         compute thisEntry(MSEMLagRegime(k,1),MSEMLagRegime(k,2))+=MSEMpt_sm(time)(k)
         compute tvemlogl+=log(p(MSEMLagRegime(k,1),MSEMLagRegime(k,2)))*MSEMpt_sm(time)(k)
      end do k
      compute w11(time)=p(1,1)*(1-p(1,1))*(thisEntry(1,1)+thisEntry(2,1))
      compute w22(time)=p(1,2)*(1-p(1,2))*(thisEntry(2,2)+thisEntry(1,2))
      compute l11(time)=(thisEntry(1,1)- p(1,1)*(thisEntry(1,1)+thisEntry(2,1)))/w11(time)
      compute l22(time)=(thisEntry(2,2)-(1-p(1,2))*(thisEntry(2,2)+thisEntry(1,2)))/w22(time)
   end do time
   *
   * Update the coefficients
   *
   linreg(noprint,weight=w11,equation=p1eq) l11 start end
   compute v1=v1+%beta
   linreg(noprint,weight=w22,equation=p2eq) l22 start end
   compute v2=v2+%beta
   compute tvemlogl=0.0
   do time=start,end
      compute p=%MSPMat(time)
      compute thisEntry=%zeros(nstates,nstates)
      do k=1,MSEMSize
         compute thisEntry(MSEMLagRegime(k,1),MSEMLagRegime(k,2))+=MSEMpt_sm(time)(k)
         compute tvemlogl+=log(p(MSEMLagRegime(k,1),MSEMLagRegime(k,2)))*MSEMpt_sm(time)(k)
      end do k
   end do time
   *
   * Check for convergence
   *
   compute testparms=v1~~v2
   compute %cvcrit=%testdiff(baseparms,testparms)
   if %cvcrit<cvcrit
      break
   compute baseparms=testparms
   disp "Iteration" emits "Log Likelihood" %emlogl %cvcrit
end do emits
end MSTVPEMEstimate


compute v1=log(.7/.3)
compute v2=log(.95/.05)

@MSTVPEMEstimate(cvcrit=.0001) gstart gend

equation p1eq *
# constant logvix{1}
equation p2eq *
# constant logvix{1}
dim p1mean(%eqnsize(p1eq)) p2mean(%eqnsize(p2eq))
do i=1,%eqnsize(p1eq)
   sstats(mean) gstart gend %eqnxvector(p1eq,t)(i)>>p1mean(i)
end do i
do i=1,%eqnsize(p2eq)
   sstats(mean) gstart gend %eqnxvector(p2eq,t)(i)>>p2mean(i)
end do i

compute v1=v1~~0.0
compute v2=v2~~0.0
@MSTVPEMEstimate(cvcrit=.0001,iters=50) gstart gend
*
* Polish estimates using BHHH
*
frml logltvtp = p=%MSPMat(t),log(%MSRegProb(t))
*
@MSFilterInit
maximize(parmset=tv+common,start=(pstar=TVPInit()),method=bhhh,iters=10) logltvtp gstart gend

[superper2008]

Dear Tom:

Appreciate you very much for your help.


Sincerely,
Feiyu

[superper2008]

Hello Tom:

Now I have got all the results by means of the code you help me figure out. Thanks again. However, when I am analyzing the model, I still need to know which state of two states in my model is recession or expansion period. In the code, you use p(1.1) and p(1.2), I am not sure if I can modify into p(1.1) and p(2,2) to identify which state is recession or expansion period. If I only get the result shown below, I can't determine the recession or expansion according to the signs of V1(1), V2(1) and V1(2), V2(2). If I can write the code P(1,1) and P(2.2), then I can identify the recession or expansion period according to the signs of them.

Appreciate you very much for your help.

Sincerely,
Feiyu

MAXIMIZE - Estimation by BHHH
Convergence in 29 Iterations. Final criterion was 0.0000081 <= 0.0000100
Usable Observations 1493
Function Value 4818.4560

Variable Coeff Std Error T-Stat Signif
************************************************************************************
1. V1(1) 33.09651919 7.41884567 4.46114 0.00000815
2. V1(2) -20.10919941 5.01057196 -4.01335 0.00005986
3. V2(1) 9.48357495 4.25329673 2.22970 0.02576738
4. V2(2) -4.71040158 2.78992367 -1.68836 0.09134173
5. BETAS(1)(1) 0.00022131 0.00024208 0.91420 0.36061336
6. BETAS(1)(2) 0.96252328 0.00499307 192.77178 0.00000000
7. BETAS(1)(3) -0.08981680 0.02266712 -3.96243 0.00007419
8. BETAS(2)(1) 0.00146803 0.00175547 0.83626 0.40301048
9. BETAS(2)(2) 0.98587491 0.01025490 96.13696 0.00000000
10. BETAS(2)(3) -0.04525605 0.02520002 -1.79587 0.07251454
11. SIGSQV(1) 0.00006727 0.00000265 25.37280 0.00000000
12. SIGSQV(2) 0.00041067 0.00004017 10.22425 0.00000000

[superper2008]

Hello Tom:

I have solved the problem. I made one mistake previously.

Feiyu

[TomDoan]

In more complicated models, you may not be able to define the regimes with simple labels like "recession" and "expansion". In the case of Hamilton's original model, the only thing that switched was the process mean for the growth rate---no matter what, one would be bigger than the other, so you could distinguish higher growth from lower growth. Instead, you have a 3 coefficient linear regression and variances which switch, so there is not necessarily any simple ordering for the two. If you're fortunate, you'll get readily identifiable regimes, but you can't count on it. (Krolzig's VAR's often have two or more local modes which are quite different in properties, but have almost identical log likelihoods).

[superper2008]

Hello Tom;
Appreciate very much for your help.

Sincerely,
Feiyu