Estimate the TVTP-MS model with MSRegression

[samson]

Dear Tom,

I estimated the FTP-MS model with @MSRegression you suggested previously, and it works!! (Thank you again for your kindly help )

Now I want to estimate the TVTP-MS model rather than FTP-MS model.

It seems that @MSRegression can also estimate the TVTP-MS model according to the description you made at
http://www.estima.com/forum/viewtopic.p ... 1208#p4302.

I found one example demonstrating the estimations with @MSRegression, but the transition probabilities in this example vary only with means.
(example:http://www.estima.com/forum/viewtopic.p ... t=30#p5340)

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

If I want to estimate the Filardo-type TVTP-MS model with @MSRegression, what modifications should I make?
(That is, how to modify the code of the example to make the transition probabilities vary with some exogenous variables? )

Best regards,
Samson

[TomDoan]

The last estimates on that have transition probabilities that vary with logvix{1}. The first set uses the TV setup to estimate a fixed transition model.

[samson]

Dear Tom:

You mean that the code above actually estimates a TVTP-MS model, right?

But something weird occurred when I relicated Filardo(1994) with the @msregression.

The codes and results from the @msregression are :

Code: Select all

    Variable                         Coeff         Std Error       T-Stat       Signif
******************************************************************************************
1.  V1(1)                             28.045926   162464.058048  1.72628e-004  0.99986226
2.  V1(2)                         -14390.350151 26757918.441531 -5.37798e-004  0.99957090
3.  V2(1)                             31.253470        0.000000       0.00000  0.00000000
4.  V2(2)                             29.059143        0.000000       0.00000  0.00000000
5.  GAMMAS(1)                          0.186422        0.023092       8.07286  0.00000000
6.  GAMMAS(2)                          0.150159        0.034618       4.33760  0.00001440
7.  GAMMAS(3)                         -0.008068        0.035597      -0.22665  0.82069296
8.  GAMMAS(4)                         -0.050465        0.031893      -1.58231  0.11357968
9.  BETAS(1)(1)                  -793779.512908        0.000000       0.00000  0.00000000
10. BETAS(2)(1)                        0.168200        0.025593       6.57216  0.00000000
11. SIGSQ                              0.417337        0.000000       0.00000  0.00000000

Code: Select all

compute gstart=1948:6,gend=1992:8
@MSRegression(states=2,NFIX=4, switch = c) g
# g{1} g{2} g{3} g{4} constant
@MSRegInitial gstart gend
@MSFilterSetup(em)


nonlin(parmset=common) gammas betas sigsq
nonlin(parmset=tv) v1 v2
dec equation p1eq   p2eq
dec vector   v1     v2
dec vector   p1mean p2mean


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 cligr{1}
equation p2eq *
# constant cligr{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 = f=%MSRegFVec(t),tvtp=%MSProb(t,f),log(tvtp)

*
@MSFilterInit
maximize(parmset=tv+common,start=(pstar=TVPInit()),method=bhhh,iters=100) logltvtp gstart gend

It seems that it does not converge at the end of the iterations.





In addition, at the last part of the code I use

Code: Select all

frml logltvtp = f=%MSRegFVec(t),tvtp=%MSProb(t,f),log(tvtp) 

rather than

Code: Select all

logltvtp = p=%MSPMat(t),log(%MSRegProb(t)

because RATS returns: ## SX11. Identifier %MSREGPROB is Not Recognizable.

How to solve this problem?



Best regards,
Samson

[samson]

Dear Tom,

I fix the weird problem described at last post. (I forgot sourcing the msregression.src.)

But the estimates of replication of Filardo(1994) are somehow different between the @msvar and the @msreg.

The codes of @msreg are:

Code: Select all

compute gstart=1948:6,gend=1992:8
@MSRegression(states=2,NFIX=4, switch = c) g
# g{1} g{2} g{3} g{4} constant
@MSRegInitial gstart gend
@MSFilterSetup(em)


nonlin(parmset=common) gammas betas sigsq
dec equation p1eq   p2eq
dec vector   v1     v2
dec vector   p1mean p2mean
nonlin(parmset=tv) v1 v2


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 cligr{1}
equation p2eq *
# constant cligr{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=100) logltvtp gstart gend

And the results from both procedures are:

Code: Select all

MSVAR

MAXIMIZE - Estimation by BHHH
Convergence in    36 Iterations. Final criterion was  0.0000076 <=  0.0000100
Monthly Data From 1948:06 To 1992:08
Usable Observations                       531
Function Value                      -595.8451

    Variable                        Coeff      Std Error      T-Stat      Signif
************************************************************************************
1.  V1(1)                         1.330180283  0.453513451      2.93306  0.00335643
2.  V1(2)                        -1.191657247  0.711571594     -1.67468  0.09399634
3.  V2(1)                         4.285992396  0.991396606      4.32319  0.00001538
4.  V2(2)                         1.649574337  0.716078802      2.30362  0.02124391
5.  MU(1)(1)                     -0.881223949  0.121836490     -7.23284  0.00000000
6.  MU(2)(1)                      0.472023399  0.060484834      7.80400  0.00000000
7.  PHI(1)(1,1)                   0.210997467  0.037889959      5.56869  0.00000003
8.  PHI(2)(1,1)                   0.072442917  0.054571339      1.32749  0.18434658
9.  PHI(3)(1,1)                   0.103731081  0.050537956      2.05254  0.04011740
10. PHI(4)(1,1)                   0.102676704  0.048535910      2.11548  0.03438913
11. SIGMA(1,1)                    0.472888920  0.027754232     17.03844  0.00000000


MSREG

MAXIMIZE - Estimation by BHHH
Convergence in    33 Iterations. Final criterion was  0.0000078 <=  0.0000100
Monthly Data From 1948:06 To 1992:08
Usable Observations                       531
Function Value                      -598.3348

    Variable                        Coeff      Std Error      T-Stat      Signif
************************************************************************************
1.  V1(1)                         0.325691305  0.578086749      0.56340  0.57316580
2.  V1(2)                        -2.819968505  1.338370453     -2.10702  0.03511616
3.  V2(1)                         2.297306819  0.872960135      2.63163  0.00849767
4.  V2(2)                         2.309537673  1.014590928      2.27632  0.02282663
5.  GAMMAS(1)                     0.237389211  0.032336567      7.34120  0.00000000
6.  GAMMAS(2)                    -0.015598865  0.043780601     -0.35630  0.72161869
7.  GAMMAS(3)                     0.032617005  0.045796841      0.71221  0.47633425
8.  GAMMAS(4)                     0.032454099  0.038069459      0.85250  0.39393832
9.  BETAS(1)(1)                  -0.367262722  0.076315622     -4.81242  0.00000149
10. BETAS(2)(1)                   0.432266265  0.043772752      9.87524  0.00000000
11. SIGSQ                         0.417336613  0.000000000      0.00000  0.00000000

The standard error and the T-statistic of the variance from @msreg are both zero. (That's really strange.)

Is there something wrong in the codes?


I'd really appreciate your help~~

Best regards,
Samson