I am trying to extend the 2-state model of Dueker(1997), "Markov Switching in GARCH Processes and Mean-Reverting Stock-Market Volatility," J of Business & Economic Statistics, vol. 15, no 1, 26-34, specifically the "dueker_swgarch_nf.rpf" (GARCH model with switching variance) code given below to 3-state:
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
Thank you very much in advance.
Best regards,
Nik