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 processin 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