The RATS Software Forum

Code: Select all

compute p(1,1)=%(z=log(%dot(%eqnxvector(p1eq,time-1),v1)),%cdf(z)
compute p(1,2)=%(z=log(%dot(%eqnxvector(p2eq,time-1),v2)),%cdf(-z)

Code: Select all

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)

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