The RATS Software Forum

Code: Select all

*
* Replication file for Dueker(1997),  "Markov Switching in GARCH
* Processes and Mean-Reverting Stock-Market Volatility," J of Business &
* Economic Statistics, vol. 15, no 1, 26-34.
*
* Data set is a reconstruction.
*
* SW-GARCH-NF model with switching normalization factors and means with
* conditionally t distributed errors.
*
open data sp500.xls
data(format=xls,org=columns,top=17,left=2,nolabels) 1 2529 spchange
*
* Copy the data over to y to make it easier to change the program.
*
set y = spchange
*
* This sets the number of regimes. The analysis requires one lagged
* regime.
*
@MSSetup(states=2,lags=1)
*
* GV will be the relative variances in the regimes.
*
dec vect gv(nstates)
*
* This will be the vector of GARCH parameters. The constant in the GARCH
* equation is fixed at 1 as the normalization.
*
dec vect msg(2)
dec real nu
*
* We have three parts of the parameter set: the mean equation parameters
* (separate mu's for each regime), the GARCH model parameters, and the
* Markov switching parameters.
*
dec vect mu(nstates)
nonlin(parmset=meanparms)  mu
nonlin(parmset=garchparms) msg gv nu
nonlin(parmset=msparms)    p
*
* uu and u are used for the series of squared residuals and the series of
* residuals needed to compute the GARCH variances.
*
clear uu u h
*
* These are used for the lagged variance for each regime.
*
compute DFFilterSize=nstates^nlags
dec vect[series] hs(DFFilterSize)
*
* These for the lagged residuals and lagged squared residuals, which are
* regime-dependent because of the regime-dependent means.
*
dec vect[series] uus(nstates) us(nstates)
************************************************************************
*
* GARCHRegimeH returns the VECTOR of the H values across augmented
* regime settings.
*
function GARCHRegimeH time
type vector GARCHRegimeH
type integer time
*
local integer i
*
dim GARCHRegimeH(nexpand)
do i=1,nexpand
   compute GARCHRegimeH(i)=1.0+msg(2)*hs(%MSLagState(i,1))(time-1)+$
       msg(1)*uus(%MSLagState(i,1))(time-1)/gv(%MSLagState(i,1))
end do i
end
**************************************************************************
function GARCHRegimeResids time
type vector 	GARCHRegimeResids
type integer	time
*
* Compute the residuals for each current regime
*
local integer i
*
dim GARCHRegimeResids(nstates)
do i=1,nstates
   compute GARCHRegimeResids(i)=y(time)-mu(i)
end do i
end
**************************************************************************
*
* Do a GARCH model to get initial guess values and matrices for random
* walk M-H.
*
linreg y
# constant y{1}
set u  = %resids
set uu = %sigmasq
do i=1,nstates
   set uus(i) = %sigmasq
   set us(i)  = %resids
end do i
*
* We'll start this at entry 2 to match the MS model
*
garch(p=1,q=1,distrib=t) 2 * y
compute gstart=%regstart(),gend=%regend()
*
ewise mu(i)=%beta(1)
compute msg=%xsubvec(%beta,3,4)
compute nu=%beta(5)
*
* The "H" series used in the Hamilton-Susmel formulation need scaling by
* the normalization factor to become variances. So we scale the
* pre-sample values down to match.
*
set h  = %sigmasq/%beta(2)
do i=1,nstates
   set hs(i) = %sigmasq/%beta(2)
end do i
*
* Start with guess values for GV that scale the original constant down
* and up.
*
compute gv(1)=%beta(3)/2.0
compute gv(2)=%beta(3)*2.0
*
compute p=||.9,.1||
**************************************************************************
*
* This does a single step of the Dueker (approximate) filter
*
function DuekerFilter time
type integer 		time
*
local real       	e
local vector      fvec evec hvec ph
*
dim fvec(nexpand) evec(nstates)
*
* Compute the "H" matrices for each expanded regime
*
compute hvec=GARCHRegimeH(time)
*
* Compute the residuals for each current regime
*
compute evec=GARCHRegimeResids(time)
*
* Compute the likelihood for each expanded regime.
*
do i=1,nexpand
   compute e=evec(MSLagState(i,1))
   compute fvec(i)=%if(hvec(i)>0,exp(%logtdensity(hvec(i)*gv(%MSLagState(i,0)),e,nu)),0.0)
end do i
*
* Do a filter step to update the probabilities. The log likelihood
* element is returned in the LOGF option.
*
@MSFilterStep(logf=DuekerFilter) time fvec
*
* Take probability weighted averages of the h's for each regime.
*
compute ph=pt_t(time).*hvec
compute pdiv  =%sumc(%reshape(pt_t(time),nstates,DFFilterSize))
compute phstar=%sumc(%reshape(ph,nstates,DFFilterSize))./pdiv
*
* Push them out into the HS vector
*
compute %pt(hs,t,phstar)
compute %pt(uus,t,evec.^2)
compute %pt(us,t,evec)
end
**************************************************************************
function DFStart
@MSFilterInit
end
**************************************************************************
frml LogDFFilter = DuekerFilter(t)
*
@MSFilterSetup(noem)
maximize(start=DFStart(),parmset=meanparms+garchparms+msparms,$
  pmethod=simplex,piters=5,method=bfgs) logDFFilter gstart gend
@MSSmoothed gstart gend psmooth
set p1 = psmooth(t)(1)
set p2 = psmooth(t)(2)
graph(max=1.0,footer="Probabilities of Regimes in Two-Regime GARCH-NF model",$
  style=stacked) 2
# p1
# p2
*
* Same model with means not switching
*
nonlin(parmset=nomeans) mu(2)=mu(1)
compute p=||.8,.2||
@MSFilterSetup(noem)
maximize(start=DFStart(),parmset=meanparms+garchparms+msparms+nomeans,$
  pmethod=simplex,piters=5,method=bfgs) logDFFilter gstart gend
@MSSmoothed gstart gend psmooth
set p1 = psmooth(t)(1)
set p2 = psmooth(t)(2)
graph(max=1.0,footer="Probabilities of Regimes in Two-Regime GARCH-NF model with fixed means",$
  style=stacked) 2
# p1
# p2

Code: Select all

*
* Replication file for Dueker(1997),  "Markov Switching in GARCH
* Processes and Mean-Reverting Stock-Market Volatility," J of Business &
* Economic Statistics, vol. 15, no 1, 26-34.
*
* Data set is a reconstruction.
*
* SW-GARCH-NF model with switching normalization factors and means with
* conditionally t distributed errors.
*
open data sp500.xls
data(format=xls,org=columns,top=17,left=2,nolabels) 1 2529 spchange
*
* Copy the data over to y to make it easier to change the program.
*
set y = spchange
*
* This sets the number of regimes. The analysis requires one lagged
* regime.
*
@MSSetup(states=3,lags=1)
*
* GV will be the relative variances in the regimes.
*
dec vect gv(nstates)
*
* This will be the vector of GARCH parameters. The constant in the GARCH
* equation is fixed at 1 as the normalization.
*
dec vect msg(2)
dec real nu
*
* We have three parts of the parameter set: the mean equation parameters
* (separate mu's for each regime), the GARCH model parameters, and the
* Markov switching parameters.
*
dec vect mu(nstates)
nonlin(parmset=meanparms)  mu
nonlin(parmset=garchparms) msg gv nu
nonlin(parmset=msparms)    p
*
* uu and u are used for the series of squared residuals and the series of
* residuals needed to compute the GARCH variances.
*
clear uu u h
*
* These are used for the lagged variance for each regime.
*
compute DFFilterSize=nstates^nlags
dec vect[series] hs(DFFilterSize)
*
* These for the lagged residuals and lagged squared residuals, which are
* regime-dependent because of the regime-dependent means.
*
dec vect[series] uus(nstates) us(nstates)
************************************************************************
*
* GARCHRegimeH returns the VECTOR of the H values across augmented
* regime settings.
*
function GARCHRegimeH time
type vector GARCHRegimeH
type integer time
*
local integer i
*
dim GARCHRegimeH(nexpand)
do i=1,nexpand
   compute GARCHRegimeH(i)=1.0+msg(2)*hs(%MSLagState(i,1))(time-1)+$
       msg(1)*uus(%MSLagState(i,1))(time-1)/gv(%MSLagState(i,1))
end do i
end
**************************************************************************
function GARCHRegimeResids time
type vector 	GARCHRegimeResids
type integer	time
*
* Compute the residuals for each current regime
*
local integer i
*
dim GARCHRegimeResids(nstates)
do i=1,nstates
   compute GARCHRegimeResids(i)=y(time)-mu(i)
end do i
end
**************************************************************************
*
* Do a GARCH model to get initial guess values and matrices for random
* walk M-H.
*
linreg y
# constant y{1}
set u  = %resids
set uu = %sigmasq
do i=1,nstates
   set uus(i) = %sigmasq
   set us(i)  = %resids
end do i
*
* We'll start this at entry 2 to match the MS model
*
garch(p=1,q=1,distrib=t) 2 * y
compute gstart=%regstart(),gend=%regend()
*
ewise mu(i)=%beta(1)
compute msg=%xsubvec(%beta,3,4)
compute nu=%beta(5)
*
* The "H" series used in the Hamilton-Susmel formulation need scaling by
* the normalization factor to become variances. So we scale the
* pre-sample values down to match.
*
set h  = %sigmasq/%beta(2)
do i=1,nstates
   set hs(i) = %sigmasq/%beta(2)
end do i
*
* Start with guess values for GV that scale the original constant down
* and up.
*
compute gv(1)=%beta(3)/4.0
compute gv(2)=%beta(3)*1.0
compute gv(3)=%beta(3)*4.0
*
compute p=|| .4,.2,.1 | .1,.2,.4 ||
**************************************************************************
*
* This does a single step of the Dueker (approximate) filter
*
function DuekerFilter time
type integer 		time
*
local real       	e
local vector      fvec evec hvec ph
*
dim fvec(nexpand) evec(nstates)
*
* Compute the "H" matrices for each expanded regime
*
compute hvec=GARCHRegimeH(time)
*
* Compute the residuals for each current regime
*
compute evec=GARCHRegimeResids(time)
*
* Compute the likelihood for each expanded regime.
*
do i=1,nexpand
   compute e=evec(%MSLagState(i,1))
   compute fvec(i)=%if(hvec(i)>0,exp(%logtdensity(hvec(i)*gv(%MSLagState(i,0)),e,nu)),0.0)
end do i
*
* Do a filter step to update the probabilities. The log likelihood
* element is returned in the LOGF option.
*
@MSFilterStep(logf=DuekerFilter) time fvec
*
* Take probability weighted averages of the h's for each regime.
*
compute ph=pt_t(time).*hvec
compute pdiv  =%sumc(%reshape(pt_t(time),nstates,DFFilterSize))
compute phstar=%sumc(%reshape(ph,nstates,DFFilterSize))./pdiv
*
* Push them out into the HS vector
*
compute %pt(hs,t,phstar)
compute %pt(uus,t,evec.^2)
compute %pt(us,t,evec)
end
**************************************************************************
function DFStart
@MSFilterInit
end
**************************************************************************
frml LogDFFilter = DuekerFilter(t)
*
@MSFilterSetup(noem)
maximize(start=DFStart(),parmset=meanparms+garchparms+msparms,$
  pmethod=simplex,piters=5,method=bfgs) logDFFilter gstart gend
@MSSmoothed gstart gend psmooth
set p1 = psmooth(t)(1)
set p2 = psmooth(t)(2)
set p3 = psmooth(t)(3)
graph(max=1.0,footer="Probabilities of Regimes in Two-Regime GARCH-NF model",$
  style=stacked) 3
# p1
# p2
# p3

Code: Select all

    Variable                        Coeff      Std Error      T-Stat      Signif
************************************************************************************
1.  MU(1)                         0.042105278  0.081061000      0.51943  0.60346295
2.  MU(2)                         0.068270121  0.048009604      1.42201  0.15502345
3.  MU(3)                         0.056200900  0.033263509      1.68957  0.09111098
4.  MSG(1)                       0.041850534  0.002716916    15.40369  0.00000000
5.  MSG(2)                       0.893396446  0.003275907   272.71724 0.00000000
6.  GV(1)                         0.023086475  0.005840585      3.95277  0.00007725
7.  GV(2)                         0.083183166  0.010420948      7.98230  0.00000000
8.  GV(3)                         0.047956776  0.004598765     10.42819  0.00000000
9.  NU                             6.419945642  0.564231852     11.37820  0.00000000
10. P(1,1)                       -0.000010000  0.003678399     -0.00272  0.99783089
11. P(2,1)                        0.017597064  0.418682882      0.04203  0.96647512
12. P(1,2)                        0.288007594  0.137790592      2.09018  0.03660134
13. P(2,2)                        0.002260262  0.235553809      0.00960  0.99234400
14. P(1,3)                        0.089021276  0.090389649      0.98486  0.32469216
15. P(2,3)                        0.339809062  0.180221421      1.88551  0.05936119

Code: Select all

compute gv(1)=%beta(3)/4.0
compute gv(2)=%beta(3)*1.0
compute gv(3)=%beta(3)*4.0
*
compute p=|| .4,.2,.1 | .1,.2,.4 ||