MSVARIANCES.RPF

MSVARIANCES.RPF is based upon an example from Kim and Nelson(1999), which does a Markov switching model of the variance of a series of returns. This is a (possible) alternative to a GARCH model, as high volatility periods will persist under the control of the switching model, thus you will see the same type of volatility clustering that you would with a GARCH model.

Kim and Nelson do a 3 regime model on the de-meaned data. The model is

${x_t} \sim N\left( {0,{\sigma ^2}({S_t})} \right)$

where ${{S_t}}$ is the variance regime. We will "label" these as 1 for the low, 2 for the medium and 3 for the high. The basic Markov Switching procedures are setup with

@MSSetup(states=3)

There are 3 parameters which govern the likelihood calculations themselves (the three variances) and the ones governing the Markov switching process—in this case, we'll use the logistic (THETA) form. The following initializes the various parameters:

dec vect sigmas(nstates)

input p

.6 .2 .1

.3 .6 .3

compute theta=%msplogistic(p)

stats ew_excs

compute sigmas(1)=0.2*%variance

compute sigmas(2)= %variance

compute sigmas(3)=5.0*%variance

Note that these are just guess values. We aren't saying that the high variance is 5 x the middle variance which is 5 x the low variance—those are designed to be relatively distinct values that are probably in a reasonable range (since they are based upon the unconditional variance of the series). The guess values for the transition probabilities have all regimes being somewhat persistent (all are probability .6 of staying). It's likely that the maximizers will be more persistent than that, but it usually works better to guess too low on the persistence than too high.

The key step in any MS estimation is a function which returns a VECTOR of (not-logged) likelihoods of the data for a single time period. That's done here with

function HStateF time e

type vector HStateF

type real e

type integer time

local integer i j

dim HStateF(nstates)

do i=1,nstates

compute HStateF(i)=exp(%logdensity(sigmas(i),e))

end do i

end

This takes as input the entry at which the calculation is being done (time, which, in this case isn't used) and the residual value (e). This is written to work with any number of states, since it uses NSTATES rather than a hard-coded count of the regimes. Note that this computes the density using %LOGDENSITY, then exp's it—most of the density functions return logged values since that's most convenient in most situations. You generally don't have to worry about this underflowing (which it probably will in some cases), because the probability will naturally be (effectively) zero in those cases. It will only be a problem if all the variances are so far out of line that all the branches end up underflowing.

With this, the log likelihood can be computed using the FRML:

frml logl = f=HStateF(t,ew_excs(t)),fpt=%MSProb(t,f),log(fpt)

f=HStateF(t,ew_excs(t)) returns and puts into F the VECTOR of likelihoods of the entry T value for the ew_excs series using the FUNCTION we just defined. %MSProb (one of the functions pulled in by @MSSetup) does a whole set of calculations for entry T:

1.it computes the predictive probabilities of the regimes (filling in the value of PT_T1(T))

2.it uses Bayes formula to combine the predictive probabilities and the likelihood elements to get the filtered probabilities of the regimes (filling in the value of PT_T(T) and of PSTAR)

3.it returns the (non-logged) likelihood of entry T

We have chosen to use the logistic index method for parameterizing the transition probabilities. With a two-regime model, there is little to be gained in doing this because the likelihood maximum will almost always be in the interior. With three-regimes, and particularly with three variance-based regimes, the maximum likelihood probabilities to jumping from the lowest to the highest or vice versa will often be 0. The theta parameterization tends to work better in those cases because large positive or negative values are able to represent the near zero probabilities without having to use inequality-constrained estimation.

nonlin(parmset=msparms) theta

nonlin(parmset=others) sigmas

@MSFilterInit

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

parmset=msparms+others,$

method=bfgs,iters=400,pmethod=simplex,piters=5) logl * 1986:12

This uses SUMMARIZE to convert the logistics into standard transition probabilities, using the delta method to produce standard errors of the estimates.

summarize(parmset=msparms+others,numerical,title="P(1<-1)") %mslogisticp(theta)(1,1)

summarize(parmset=msparms+others,numerical,title="P(2<-1)") %mslogisticp(theta)(2,1)

summarize(parmset=msparms+others,numerical,title="P(3<-1)") 1-%mslogisticp(theta)(1,1)-%mslogisticp(theta)(2,1)

There are similar sets for the other probabilities.

This computes and graphs the smoothed probabilities of the regimes.

@MSSmoothed %regstart() %regend() psmooth

set p1 = psmooth(t)(1)

set p2 = psmooth(t)(2)

set p3 = psmooth(t)(3)

graph(style=stacked,maximum=1.0,picture="##.##",header="Smoothed Probabilities of Variance States",$

key=below,klabels=||"Low Variance","Medium Variance","High Variance"||) 3

# p1

# p2

# p3

This computes the (smoothed) estimate of the variances. Ex-ante estimates of the variances (which would have a similar interpretation of the conditional variance of a GARCH process) would use the probabilities computed from pt_t1 rather than psmooth.

set variance = p1*sigmas(1)+p2*sigmas(2)+p3*sigmas(3)

graph(footer="Figure 4.10 Estimated variance of historical stock returns")

# variance

Full Program

*
* Based upon Application 3 from pp 86-92 of Kim and Nelson, "State-space
* Models with Regime Switching"
*
open data ew_excs.prn
calendar(m) 1926:1
data(format=free,org=columns) 1926:01 1986:12 ew_excs
*
* Extract the mean from the returns.
*
diff(center) ew_excs
*
@MSSetup(states=3)
dec vect sigmas(nstates)
*
input p
.6 .2 .1
.3 .6 .3
compute theta=%msplogistic(p)
*
stats ew_excs
compute sigmas(1)=0.2*%variance
compute sigmas(2)= %variance
compute sigmas(3)=5.0*%variance
*********************************************************************
* HStateF returns a vector of likelihoods for the various states at
* time given residual e. The likelihoods differ in the states based upon
* the values of sigmas.
*
function HStateF time e
type vector HStateF
type real e
type integer time
*
local integer i j
*
dim HStateF(nstates)
do i=1,nstates
compute HStateF(i)=exp(%logdensity(sigmas(i),e))
end do i
end
*********************************************************************
frml logl = f=HStateF(t,ew_excs(t)),fpt=%MSProb(t,f),log(fpt)
*
nonlin(parmset=msparms) theta
nonlin(parmset=others) sigmas
@MSFilterInit
maximize(start=%(p=%mslogisticp(theta),pstar=%MSInit()),$
parmset=msparms+others,$
method=bfgs,iters=400,pmethod=simplex,piters=5) logl * 1986:12
*
* Use delta method to convert theta values into P's with standard errors
*
summarize(parmset=msparms+others,numerical,title="P(1<-1)") %mslogisticp(theta)(1,1)
summarize(parmset=msparms+others,numerical,title="P(2<-1)") %mslogisticp(theta)(2,1)
summarize(parmset=msparms+others,numerical,title="P(3<-1)") 1-%mslogisticp(theta)(1,1)-%mslogisticp(theta)(2,1)
*
summarize(parmset=msparms+others,numerical,title="P(1<-2)") %mslogisticp(theta)(1,2)
summarize(parmset=msparms+others,numerical,title="P(2<-2)") %mslogisticp(theta)(2,2)
summarize(parmset=msparms+others,numerical,title="P(3<-2)") 1-%mslogisticp(theta)(1,2)-%mslogisticp(theta)(2,2)
*
summarize(parmset=msparms+others,numerical,title="P(1<-3)") %mslogisticp(theta)(1,3)
summarize(parmset=msparms+others,numerical,title="P(2<-3)") %mslogisticp(theta)(2,3)
summarize(parmset=msparms+others,numerical,title="P(3<-3)") 1-%mslogisticp(theta)(1,3)-%mslogisticp(theta)(2,3)
*
@MSSmoothed %regstart() %regend() psmooth
set p1 = psmooth(t)(1)
set p2 = psmooth(t)(2)
set p3 = psmooth(t)(3)
graph(style=stacked,maximum=1.0,picture="##.##",header="Smoothed Probabilities of Variance States",$
key=below,klabels=||"Low Variance","Medium Variance","High Variance"||) 3
# p1
# p2
# p3
*
set variance = p1*sigmas(1)+p2*sigmas(2)+p3*sigmas(3)
graph(footer="Figure 4.10 Estimated variance of historical stock returns")
# variance
*
set stdu = ew_excs/sqrt(variance)
graph(footer="Figure 4.11a Plot of standardized stock returns")
# stdu

Output

MAXIMIZE - Estimation by BFGS

Convergence in 65 Iterations. Final criterion was 0.0000013 <= 0.0000100

Monthly Data From 1926:01 To 1986:12

Usable Observations 732

Function Value 1001.8949

Variable Coeff Std Error T-Stat Signif

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

1. THETA(1,1) 66.67235534 12.97392298 5.13895 0.00000028

2. THETA(2,1) 63.02286265 13.06145554 4.82510 0.00000140

3. THETA(1,2) 0.50932677 0.86531018 0.58861 0.55612560

4. THETA(2,2) 4.41636496 0.63821800 6.91984 0.00000000

5. THETA(1,3) -58.05123899 26.62872889 -2.18002 0.02925576

6. THETA(2,3) -2.95954238 0.67969578 -4.35422 0.00001335

7. SIGMAS(1) 0.00123003 0.00017970 6.84486 0.00000000

8. SIGMAS(2) 0.00399421 0.00047418 8.42343 0.00000000

9. SIGMAS(3) 0.03111257 0.00569573 5.46244 0.00000005

P(1<-1)

Value 0.97465477 t-Statistic 61.92307

Standard Error 0.01573977 Signif Level 0.0000000

P(2<-1)

Value 0.02534523 t-Statistic 1.61027

Standard Error 0.01573977 Signif Level 0.1073395

P(3<-1)

Value 4.163e-17 t-Statistic 3.07783e-05

Standard Error 1.353e-12 Signif Level 0.9999754

P(1<-2)

Value 0.01947333 t-Statistic 1.66493

Standard Error 0.01169621 Signif Level 0.0959274

P(2<-2)

Value 0.96882514 t-Statistic 69.67795

Standard Error 0.01390433 Signif Level 0.0000000

P(3<-2)

Value 0.01170153 t-Statistic 1.58644

Standard Error 0.00737596 Signif Level 0.1126394

P(1<-3)

Value 5.844e-26 t-Statistic 0.03759

Standard Error 1.555e-24 Signif Level 0.9700154

P(2<-3)

Value 0.04928744 t-Statistic 1.54752

Standard Error 0.03184932 Signif Level 0.1217380

P(3<-3)

Value 0.95071256 t-Statistic 29.85033

Standard Error 0.03184932 Signif Level 0.0000000

Graphs