Hamilton Susmel JOE 1994

Hamilton and Susmel(1994) is one of two papers independently published in 1994 that proposed different methods of implementing an Markov-Swiitching ARCH model (MS-ARCH) to explain the high degree of persistence in volatility in financial time series. (Too high to be explained by a simple ARCH model; this is a competitor for the GARCH model). The alternative is Cai(1994) which is demonstrated in the SWARCH.RPF example.

In the Hamilton-Susmel formulation, the variance in regime $S(t)$ is $g\left( {S(t)} \right)h(t)$ where $h(t)$ is an ARCH process in the regime-specific standardized residuals $u\left( {t - j} \right)^2 /g\left( {S\left( {t - j} \right)} \right)$. This makes the variance a function of the current and $q$ lags of the regime where $q$ is the number of ARCH lags.

The H-S replication has three program files:

hs_garch.rpf

Standard ARCH/GARCH models

hs_swarch_gl22.rpf

Two regime-two lag model with Gaussian errors and asymmetry

hs_swarch_tl32.rpf

Three regime-two lag model with Student-t errors and asymmetry

The data series is roughly 25 years of weekly value-weighted returns on the US stock market. (Data are already in returns).

The first of these just does various standard GARCH instructions with various combinations of options, such as asymmetry and t-distributed errors.

The SWARCH setup in the GL22 program starts with:

compute q=2

@MSSetup(lags=q,regimes=2)

Q is the number of "ARCH" lags, and this sets up the model with 2 regimes. (The other program will use REGIMES=3). Because the variance process depends upon the lagged regimes (since lagged squared u's are normalized by the regime at the time), the @MSSetup needs the LAGS option.

GV will be the variance scales in the regimes:

dec vect gv(nregimes)

The next is the vector of ARCH parameters, the asymmetry parameter (if it's used) and the degrees of freedom for the t errors (if they are used). The constant in the ARCH equation is fixed at 1 as the normalization. (This differs from the normalization in the paper which normalizes one of the relative variances. The two are equivalent, though this is easier to implement).

dec vect a(q)

dec real xi

dec real nu

We have three parts of the parameter set: the mean equation parameters (here, an intercept and an AR lag), the ARCH model parameters, here including an asymmetry term, but not the degrees of freedom parameter, and the Markov switching parameters (which will be done using the logistic parameterization).

nonlin(parmset=meanparms) alpha beta

nonlin(parmset=archparms) a gv xi

nonlin(parmset=msparms) theta

UU and U are used for the series of squared residuals and the series of residuals needed to compute the ARCH variances.

Each of the SWARCH programs includes separate FUNCTION's for evaluating the likelihoods (note: likelihood not log likelihood) for a given time period TIME and computed residual E for Gaussian and for $t$ distributed errors. This is ARCHRegimeGaussF for the Gaussian/normal errors:

function ARCHRegimeGaussF time e

type vector ARCHRegimeTF

type real e

type integer time

local integer i k

local real hi

dim ARCHRegimeGaussF(nexpand)

do i=1,nexpand

compute hi=1.0

do k=1,q

compute hi=hi+a(k)*uu(time-k)/gv(%MSLagRegime(i,k))

end do k

compute hi=hi+xi*%if(u(time-1)<0.0,uu(time-1)/gv(%MSLagRegime(i,1)),0.0)

compute ARCHRegimeGaussF(i)=%if(hi>0,exp(%logdensity(hi*gv(%MSLagRegime(i,0)),e)),0.0)

end do i

end

Note the use of the %MSLagRegime(i,lag) functions. This is needed when dealing with models that use lagged regimes, because all the calculations need to be done over all combinations of current and (in this case) two lagged regimes. The NEXPAND variable is here equal to 8, which is the number of expanded regimes (2 regimes to the power 3 (for current + 2 lags)). %MSLagRegime(i,lag) returns the actual regime at lag number LAG (which includes LAG=0 for the current time period) for expanded regime number I.

The inner DO loop:

compute hi=1.0

do k=1,q

compute hi=hi+a(k)*uu(time-k)/gv(%MSLagRegime(i,k))

end do k

computes the ARCH model on the standardized squared residuals:

\begin{equation} 1 + \sum\limits_{k = 1}^q {a(k)u_{t - k}^2 /g\left( {S\left( {t - k} \right)} \right)} \end{equation}

The next line does an asymmetry adjustment using just the 1st lag. XI is the adjustment coefficient. Given the sign convention in this, a positive value of XI increases volatility for negative residuals (it's a GJR-style adjustment). If XI is pegged as zero, there is no asymmetry adjustment.

compute hi=hi+xi*%if(u(time-1)<0.0,uu(time-1)/gv(%MSLagRegime(i,1)),0.0)

Finally,

compute ARCHRegimeGaussF(i)=%if(hi>0,exp(%logdensity(hi*gv(%MSLagRegime(i,0)),e)),0.0)

computes the likelihood for expanded regime I by scaling the standardized ARCH process value (HI) by the current variance scale factor $g\left( {S(t)} \right)$ to get the regime-specific variance, evaluates the log density for the Normal for that variance and the input residual E. That's exp'ed to give the likelihood. The %IF checks if the HI is positive; if it isn't the probability returns zero, which will force a missing value when converted to logs in doing the log likelihood.

ARCHRegimeTF is almost identical except for that final line which uses %logtdensity (with the added parameter of the NU) rather than %logdensity.

This runs a linear regression to get guess values for the AR1 mean model:

linreg y

# constant y{1}

compute alpha=%beta(1),beta=%beta(2)

As is typically the case with Markov switching models, there is no global identification of the regimes. By defining a fairly wide spread for GV, we'll hope that we'll stay in the zone of the likelihood where regime 1 is low variance and regime 2 is high. In this setup for the SWARCH, the "ARCH" part is scale-free, so the actual variance scale will be in the GV, so these are multiples of the variance from the LINREG.

compute gv=%seesq*||0.2,5.0||

These are our guess values for the P matrix. (We don't have any really strong information, so both regimes as set to be moderately persistent). We have to invert that to get the guess values for the logistic indexes.

compute p=||.8,.2||

compute theta=%msplogistic(p)

This defines the FRML for the residual:

frml uf = y-alpha-beta*y{1}

This gives guess values for the ARCH parameters (small positive numbers), zeros out the asymmetry parameter XI and initializes the UU to the residual variance, and U to the linear regression residuals.

compute a=%const(0.05),xi=0.0

set uu = %seesq

set u = %resids

This is the range over which we can run the SWARCH model: from the linear regression, we lose Q data points for the lagged residuals:

compute gstart=%regstart()+q,gend=%regend()

We need to keep series of the residual (U) and squared residual (UU). Because the mean function is the same across regimes, we can just compute the residual and send it to ARCHRegimeGaussF. %MSProb takes the VECTOR of likelihoods and updates (internally) the probabilities and returns the probability-weighted likelihood. The end calculation in the LOGL FRML is the log of that.

frml logl = u(t)=uf(t),uu(t)=u(t)^2,f=ARCHRegimeGaussF(t,u(t)),$

fpt=%MSProb(t,f),log(fpt)

@MSFilterInit makes sure all the bookkeeping is set. Then the MAXIMIZE estimates the model. The START option here transforms the logistic parameters THETA to create the standard probabilities P, then computes the initial expanded regime probabilities into PSTAR.

@MSFilterInit

maximize(parmset=meanparms+archparms+msparms,$

start=%(p=%mslogisticp(theta),pstar=%MSInit()),$

method=bfgs,iters=400,pmethod=simplex,piters=5) logl gstart gend

And finally, this computes and graphs the smoothed probabilities of the regimes:

@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 Gaussian model",$

style=stacked) 2

# p1

# p2