MSVAR TVTP Problems

[Jules89]

Hello Tom,
I want to estimate an Hamilton Style MS-AR(1) with TVTP by using EM. Therefore I used the replication code of FIlardo(1994). The Filardo code I use is the following:

Code: Select all

****************************************************************************************
************************* Reading in and Transforming the Data *************************
****************************************************************************************


cal(m) 1948
open data filardo.dat
data(format=prn,org=cols) 1948:1 1992:8 year month ip cli dfi xli sp fed tspread


set ipgr = 100.0*log(ip/ip{1})

stats(noprint) ipgr * 1959:12
compute stdearly=sqrt(%variance)
stats(noprint) ipgr 1960:1 *
compute stdlate=sqrt(%variance)

set ipgradjust = %if(t<=1959:12,ipgr*stdlate/stdearly,ipgr)
set g = ipgradjust
set cligr   = 100.0*log(cli/cli{1})
set spgr    = 100.0*log(sp/sp{1})
set xlidiff = xli-xli{1}
set ffdiff  = fed-fed{1}
set tspdiff = tspread-tspread{1}
set dfismooth = dfi+2*dfi{1}+2*dfi{2}+dfi{3}


* Center variabes to zero means
dofor s = cligr spgr xlidiff ffdiff tspdiff dfismooth
   diff(center) s
end dofor s

*********************************************************************************************
**************************** Setting up the Markov Switching CTP ****************************
*********************************************************************************************

source msvarsetup.src

compute nlags=1

compute gstart=1948:06,gend=1992:8
@MSVARSetup(lags=nlags)
# g

nonlin(parmset=common) mu phi sigma
nonlin(parmset=fixed)  p

frml loglftp = log(%MSVarProb(t))

@MSVARInitial gstart gend


*************************************************************************************
**************************** Setting up the EM-Algorithm ****************************
*************************************************************************************


@msvarEMgeneralsetup

*** Inserting Guess Values
compute mu(1)=-1.7,mu(2)=.3
compute p(1,1)=.3,p(1,2)=.1


@msvartestvector(dop,domu,dovar) baseparms
do emits=1,100
   @msvaremstep gstart gend
   @msvartestvector(dop,domu,dovar) testparms
   compute %cvcrit=%testdiff(baseparms,testparms)
   if %cvcrit<.0001
      break
   compute baseparms=testparms
   disp "Iteration" emits "Log Likelihood" %logl %cvcrit
end do emits

**************************************************************************
**************************** Maximization CTP ****************************
**************************************************************************


maximize(start=(pstar=%MSVARInit()),parmset=fixed+common,$
  method=bhhh,iters=20) loglftp gstart gend



**********************************************************************
******************** Time varying probabilities **********************
**********************************************************************

dec equation p1eq p2eq
dec vector   v1   v2
nonlin(parmset=tv) v1 v2

**********************************************************************
********** Overwriting the standard fixed transition matrix **********
**********************************************************************

function  %MSVARPmat  time
type rect %MSVARPmat
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 %MSVARInitTransition()
*
if nexpand==nstates {
   dim pexpand(nstates,nstates-1)
   ewise pexpand(i,j)=%if(i==nstates,-1.0,i==j)
   compute %MSVARPmat=pexpand*p
   ewise %MSVARPmat(i,j)=%if(i==nstates,1.0+%MSVARPmat(i,j),%MSVARPmat(i,j))
}
else {
   dim %MSVARPmat(nexpand,nexpand)
   ewise %MSVARPmat(i,j)=MSVARTransProbs(MSVARTransLookup(i,j))
}
end

**********************************************************************
******* Initialize the probabilities using only the intercepts *******
**********************************************************************

function TVPInit
type vector TVPInit
*
local rect a
local integer i j
*
compute p=||%logistic(v1(1),1.0),1-%logistic(v2(1),1.0)||
compute TVPInit=%MSVARInit()
end

*******************************************************************************
***** Controlling routine for EM estimation of time-varying probabilities *****
*******************************************************************************

procedure MSVARTVPEMEstimate 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
*
@msvartestvector(domu,dovar) baseparms
compute baseparms=baseparms~~v1~~v2
do emits=1,iters
   compute %eqnsetcoeffs(p1eq,v1)
   compute %eqnsetcoeffs(p2eq,v2)
   @msvarestep start end
   compute %emlogl=%logl
   @msvarmstep(nodop) start end
   *
   * Compute the weights and dependent variable for the one-step
   * correction on the logistics.
   *
   set l11 start end = 0.0
   set l22 start end = 0.0
   set w11 start end = 0.0
   set w22 start end = 0.0
   do time=start,end
      compute %MSVARPMat(time)
      compute thisEntry=%zeros(nstates,nstates)
      do k=1,MSVARFilterSize
         compute thisEntry(MSVARLagState(k,1),MSVARLagState(k,2))+=MSVARPTSM(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
   *
   * Check for convergence
   *
   @msvartestvector(domu,dovar) testparms
   compute testparms=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 MSVARTVPEMEstimate

**********************************************************************
********************* Calculating the Likelihood *********************
**********************************************************************

frml logltvtp = %MSVARPMat(t),log(%MSVarProb(t))

**********************************************************************
****************** Initialize the other Parameters *******************
**********************************************************************

equation p1eq *
# constant cligr{1}
equation p2eq *
# constant cligr{1}

@MSVARInitial gstart gend

compute v1=log(.7/.3)~~0.0
compute v2=log(.95/.05)~~0.0
compute mu(1)=-1.7,mu(2)=.3

@MSVARTVPEMEstimate(cvcrit=.0001) gstart gend

maximize(parmset=tv+common,start=(pstar=TVPInit()),$
 method=bhhh) logltvtp gstart gend

**********************************************************************
********************* Various other combinations *********************
**********************************************************************

equation p1eq *
# constant cligr{1 2}
equation p2eq *
# constant cligr{1 2}
*
compute v1=v1~~0.0
compute v2=v2~~0.0

@MSVARTVPEMEstimate(cvcrit=.0001) gstart gend

maximize(parmset=tv+common,start=(pstar=TVPInit()),$
 method=bhhh) logltvtp gstart gend


*********************************************************************

equation p1eq *
# constant xli{1}
equation p2eq *
# constant xli{1}

compute mu(1)=-1.7,mu(2)=.33
compute v1=log(.7/.3)~~0.0
compute v2=log(.95/.05)~~0.0

@MSVARTVPEMEstimate(cvcrit=.0001) gstart gend

maximize(parmset=tv+common,start=(pstar=TVPInit()),$
 method=bhhh,iters=20) logltvtp gstart gend

*********************************************************************

equation p1eq *
# constant xli{1 2}
equation p2eq *
# constant xli{1 2}

compute mu(1)=-1.7,mu(2)=.33
compute v1=log(.7/.3)~~0.0~~0.0
compute v2=log(.95/.05)~~0.0~~0.0

@MSVARTVPEMEstimate(cvcrit=.0001) gstart gend

maximize(parmset=tv+common,start=(pstar=TVPInit()),$
 method=bhhh,iters=20) logltvtp gstart gend

*********************************************************************

equation p1eq *
# constant spgr{1}
equation p2eq *
# constant spgr{1}

compute mu(1)=-1.7,mu(2)=.33
compute v1=log(.7/.3)~~0.0
compute v2=log(.95/.05)~~0.0

@MSVARTVPEMEstimate(cvcrit=.0001) gstart gend

maximize(parmset=tv+common,start=(pstar=TVPInit()),$
 method=bhhh) logltvtp gstart gend

*********************************************************************

equation p1eq *
# constant ffdiff{1}
equation p2eq *
# constant ffdiff{1}

compute mu(1)=-1.7,mu(2)=.33
compute v1=log(.7/.3)~~0.0
compute v2=log(.95/.05)~~0.0

@MSVARTVPEMEstimate(cvcrit=.0001) gstart gend

maximize(parmset=tv+common,start=(pstar=TVPInit()),$
 method=bhhh,iters=20) logltvtp gstart gend

*********************************************************************

equation p1eq *
# constant tspdiff{1 2}
equation p2eq *
# constant tspdiff{1 2}

compute mu(1)=-1.7,mu(2)=.33
compute v1=log(.7/.3)~~0.0~~0.0
compute v2=log(.95/.05)~~0.0~~0.0

@MSVARTVPEMEstimate(cvcrit=.0001) gstart gend

maximize(parmset=tv+common,start=(pstar=TVPInit()),$
 method=bhhh,iters=20) logltvtp gstart gend

It works fine. Now I want to use it to estimate my model. It is a AR(1) with switching mean and variance. Here is the code I use:
(Its basicly the Filardo code with switching mean and varaince and other data)

Code: Select all

*********************************************************
*********** Reading in Data and Renaming Them ***********
*********************************************************
OPEN DATA "my directory"
CALENDAR(D) 2013:01:01
DATA(FORMAT=XLSX,ORG=COLUMNS) 2013:01:01 2015:12:15 skew s



*********************************************************************************************
**************************** Setting up the Markov Switching CTP ****************************
*********************************************************************************************

source msvarsetup.src

compute nlags=1

compute gstart=2013:01:01
compute gend=2015:12:15


display gstart
display gend
@MSVARSetup(lags=nlags,switch=mh)
# s

nonlin(parmset=common) mu phi sigmav
nonlin(parmset=fixed)  p

frml loglftp = log(%MSVarProb(t))

@MSVARInitial gstart gend

*************************************************************************************
**************************** Setting up the EM-Algorithm ****************************
*************************************************************************************

@msvarEMgeneralsetup

@msvartestvector(dop,domu,dovar) baseparms
do emits=1,100
   @msvaremstep gstart gend
   @msvartestvector(dop,domu,dovar) testparms
   compute %cvcrit=%testdiff(baseparms,testparms)
   if %cvcrit<.00001
      break
   compute baseparms=testparms
   disp "Iteration" emits "Log Likelihood" %logl %cvcrit
end do emits

**************************************************************************
**************************** Maximization CTP ****************************
**************************************************************************


maximize(start=(pstar=%MSVARInit()),parmset=fixed+common,$
  method=bhhh) loglftp gstart gend



**********************************************************************
******************** Time varying probabilities **********************
**********************************************************************

dec equation p1eq p2eq
dec vector   v1   v2
nonlin(parmset=tv) v1 v2

**********************************************************************
********** Overwriting the standard fixed transition matrix **********
**********************************************************************

function  %MSVARPmat  time
type rect %MSVARPmat
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 %MSVARInitTransition()
*
if nexpand==nstates {
   dim pexpand(nstates,nstates-1)
   ewise pexpand(i,j)=%if(i==nstates,-1.0,i==j)
   compute %MSVARPmat=pexpand*p
   ewise %MSVARPmat(i,j)=%if(i==nstates,1.0+%MSVARPmat(i,j),%MSVARPmat(i,j))
}
else {
   dim %MSVARPmat(nexpand,nexpand)
   ewise %MSVARPmat(i,j)=MSVARTransProbs(MSVARTransLookup(i,j))
}
end

**********************************************************************
******* Initialize the probabilities using only the intercepts *******
**********************************************************************

function TVPInit
type vector TVPInit
*
local rect a
local integer i j
*
compute p=||%logistic(v1(1),1.0),1-%logistic(v2(1),1.0)||
compute TVPInit=%MSVARInit()
end

*******************************************************************************
***** Controlling routine for EM estimation of time-varying probabilities *****
*******************************************************************************

procedure MSVARTVPEMEstimate 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
*
@msvartestvector(domu,dovar) baseparms
compute baseparms=baseparms~~v1~~v2
do emits=1,iters
   compute %eqnsetcoeffs(p1eq,v1)
   compute %eqnsetcoeffs(p2eq,v2)
   @msvarestep start end
   compute %emlogl=%logl
   @msvarmstep(nodop) start end
   *
   * Compute the weights and dependent variable for the one-step
   * correction on the logistics.
   *
   set l11 start end = 0.0
   set l22 start end = 0.0
   set w11 start end = 0.0
   set w22 start end = 0.0
   do time=start,end
      compute %MSVARPMat(time)
      compute thisEntry=%zeros(nstates,nstates)
      do k=1,MSVARFilterSize
         compute thisEntry(MSVARLagState(k,1),MSVARLagState(k,2))+=MSVARPTSM(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
   *
   * Check for convergence
   *
   @msvartestvector(domu,dovar) testparms
   compute testparms=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 MSVARTVPEMEstimate

**********************************************************************
********************* Calculating the Likelihood *********************
**********************************************************************

frml logltvtp = %MSVARPMat(t),log(%MSVarProb(t))

**********************************************************************
****************** Initialize the other Parameters *******************
**********************************************************************

equation p1eq *
# constant skew{1}
equation p2eq *
# constant skew{1}

@MSVARInitial gstart gend

compute v1=log(0.95/(1-0.95))~~0.0
compute v2=log(0.95/(1-0.95))~~0.0

@MSVARTVPEMEstimate(cvcrit=.000001) gstart gend

maximize(parmset=tv+common,start=(pstar=TVPInit()),$
 method=bhhh) logltvtp gstart gend

I get the following two errors:
1) When I start

Code: Select all

@msvarEMgeneralsetup

@msvartestvector(dop,domu,dovar) baseparms
do emits=1,100
   @msvaremstep gstart gend
   @msvartestvector(dop,domu,dovar) testparms
   compute %cvcrit=%testdiff(baseparms,testparms)
   if %cvcrit<.00001
      break
   compute baseparms=testparms
   disp "Iteration" emits "Log Likelihood" %logl %cvcrit
end do emits

nothing happens (Why is that). When I then maximize I get

Code: Select all

## NL6. NONLIN Parameter P(1,1) Has Not Been Initialized. Trying 0
## NL6. NONLIN Parameter P(1,2) Has Not Been Initialized. Trying 0
## NL6. NONLIN Parameter MU(1)(1) Has Not Been Initialized. Trying 0
## NL6. NONLIN Parameter MU(2)(1) Has Not Been Initialized. Trying 0
## NL6. NONLIN Parameter PHI(1)(1,1) Has Not Been Initialized. Trying 0
## NL6. NONLIN Parameter SIGMAV(1)(1,1) Has Not Been Initialized. Trying 0
## NL6. NONLIN Parameter SIGMAV(2)(1,1) Has Not Been Initialized. Trying 0

2) When I want to estimate the model with TVTP I get:

Code: Select all

## SR10. Missing Values And/Or SMPL Options Leave No Usable Data Points
The Error Occurred At Location 961, Line 41 of MSVARTVPEMESTIMA
35445160 Position 3539

I really don't understand why all of that happens, I use basicly the same code as in the replication code. Even If the model is not well suited it should produce something...

Can you please have a look at all those codes and tell me what I did wrong? I am lost in the forest^^

The data are attached .

Thank you very much in advance

[TomDoan]

Your GSTART isn't allowing for the data point lost to a lag, so the function doesn't exist at the start, and that cascades through. If you correct that, it will at least try to estimate something.

However, your major problem is that your series is completely inappropriate for a MS model. It's almost a textbook case of what should be an intervention model. It has a single obvious break at an obvious location---the only question is the dynamics surrounding the adaptation to the new situation. That's what an intervention model is for.

[Jules89]

What do you mean by "GSTART isn't allowing for the data point lost to a lag" ?
I tried many gstart values. The first error does not show up, but the second when I want to estimate the tvtp produces the same error.
It looks like this

Code: Select all

Iteration 1 Log Likelihood    2842.50573       0.35750
Iteration 2 Log Likelihood    3212.36367       0.28460
Iteration 3 Log Likelihood    3222.43631       0.53260
Iteration 4 Log Likelihood    3225.21527       0.54840
Iteration 5 Log Likelihood    3226.15429       0.18753
Iteration 6 Log Likelihood    3226.39767       0.12465
Iteration 7 Log Likelihood    3226.45390       0.11081
Iteration 8 Log Likelihood    3226.47204       0.09975
Iteration 9 Log Likelihood    3226.47870       0.09070
Iteration 10 Log Likelihood    3226.48115       0.08315
Iteration 11 Log Likelihood    3226.48062       0.07677
Iteration 12 Log Likelihood    3226.48043       0.07130
Iteration 13 Log Likelihood    3226.48036       0.06655
Iteration 14 Log Likelihood    3226.48033       0.06240
Iteration 15 Log Likelihood    3226.48032       0.05873
Iteration 16 Log Likelihood    3226.48032       0.05548
Iteration 17 Log Likelihood    3226.48032       0.05256
Iteration 18 Log Likelihood    3226.48032       0.04994
Iteration 19 Log Likelihood    3226.48032       0.04756
Iteration 20 Log Likelihood    3226.48032       0.04540
Iteration 21 Log Likelihood    3226.48032       0.04343
Iteration 22 Log Likelihood    3226.48032       0.04162
Iteration 23 Log Likelihood    3226.48032       0.03996
Iteration 24 Log Likelihood    3226.48032       0.03842
Iteration 25 Log Likelihood    3226.48032       0.03700
Iteration 26 Log Likelihood    3226.48032       0.03568
Iteration 27 Log Likelihood    3226.48032       0.03445
Iteration 28 Log Likelihood    3226.48032       0.03330
Iteration 29 Log Likelihood    3226.48032       0.03223
Iteration 30 Log Likelihood    3226.48032       0.03122
Iteration 31 Log Likelihood    3226.48032       0.03028
Iteration 32 Log Likelihood    3226.48032       0.02939
Iteration 33 Log Likelihood    3226.48032       0.02855
Iteration 34 Log Likelihood    3226.48032       0.02776
Iteration 35 Log Likelihood    3226.48032       0.02701
## SR10. Missing Values And/Or SMPL Options Leave No Usable Data Points
The Error Occurred At Location 961, Line 41 of MSVARTVPEMESTIMA
11655592 Position 3539

[TomDoan]

You had GSTART set to observation 1, which doesn't leave data for the lag. I just explained that if you fix that, it will try to estimate something, but will fail. It tries, but eventually discovers that regime 1 is absorbing. You can't fit a MS model to that data set.