HMM Procedure

Consider the following CRSP NYSE/NYSEMKT/Nasdaq/Arca Value-Weighted Market Index (Center for Research in Security Prices 2018) from January 19, 1926, to December 29, 2017. You can access the CRSP data through the SASECRSP interface engine; for more information, see the chapter "The SASECRSP Interface Engine" in SAS/ETS User's Guide. You can also directly access the data if you install the CRSP data locally. If you install the CRSP data in the folder XXX, the following statements extract the daily stock prices (that is, market indices) into the data set vwmiO:

***************************************************************************;
* DATA                                                                     ;
***************************************************************************;

title 'Discovering the Hidden Market States';
libname crspInd "XXX";
data vwmiO;
   set crspInd.sfz_dind(where=(KYINDNO=1000200
      and caldt>='19JAN1926'd and caldt<='29DEC2017'd));
   format date MMDDYY10.;
   date = caldt; price = aind;
   keep date price;
run;

The stock-price series is nonstationary and not suitable to be directly modeled; hence, the stationary weekly-return series is modeled. The following code generates the weekly returns through the daily stock prices and saves them in two data sets. One data set, vwmi, is the full-sample data set, which contains all weekly returns. The other data set, vwmiIn, is the in-sample data set, which contains weekly returns before December 31, 2000. The weekly returns after January 1, 2001, are in the out-of-sample period and are used for testing the models.

%let ds = vwmi;
%let cutDate =  '31DEC2000'd;
data &ds. &ds.In;
   set &ds.O;
   retain cumReturn 0;
   return = (log(price)-log(lag(price)))*100;
   if(return~=.) then cumReturn + return;
   if(mod(_N_,5)=1 and _N_>1) then do;
      returnw = cumReturn;
      w + 1;
      output &ds.;
      if (date<=&cutDate.) then output &ds.In;
      cumReturn = 0;
   end;
   keep w date returnw price;
   label w="Week Index" date="Date" returnw="Weekly Return"
         price="Market Index";
run;

The following statements plot the stock prices and weekly returns, as shown in Output 20.1.1. The in-sample and out-of-sample periods are separated by a vertical line.

proc sgplot data=&ds.;
   series x=date y=returnw;
   series x=date y=price / y2axis;
   refline 0 / axis=y lineattrs=(color=gray thickness=1);
   refline &cutDate. / axis=x
                       lineattrs=(color=gray thickness=2 pattern=dash);
   yaxis max=20 min=-20;
run;

Output 20.1.1: Stock Prices and Weekly Returns

Before you estimate any model, make a copy of the data in the library, as in the following DATA steps:

   data mylib.&ds.In; set &ds.In; run;
   data mylib.&ds.; set &ds.; run;

Consider a simple case. First, assume that there are two hidden market states, bear market state and bull market state, which can be modeled by exactly two corresponding hidden states in an HMM; that is, consider a bistate HMM. Second, assume that the lagged weekly returns cannot help predict current weekly return at each hidden state; that is, the autoregressive order is zero at each hidden state. Hence, as a start, consider a bistate RS-AR(0) model. The following statements estimate the bistate RS-AR(0) model on weekly returns. The input data table, specified in the required DATA= option in the PROC HMM statement, contains the in-sample data. In order to estimate the model, the ID and MODEL statements must be specified. Since the data are sorted by the variable date, the variable date is assigned to the required TIME= option. However, since the data can also be regarded as being sorted by the variable w, the variable w, instead of the variable date, can be assigned to the TIME= option; if you do so, there is no impact on the model estimation. In the MODEL statement, the dependent variable, the weekly return returnw, is specified first; the TYPE= option is set to AR for the RS-AR model; the NSTATE= option is set to 2 for the bistate model; and the YLAG= option is set to 0 for an AR(0) model. The SMOOTH, DECODE, and EVALUATE statements specify where to output, respectively, the smoothed probabilities, decoded state path, and log likelihood up to each time point.

***************************************************************************;
* ANALYSIS WITH A SIMPLE MODEL                                             ;
***************************************************************************;

proc hmm data=mylib.&ds.In;
   id time=date;
   model returnw / type=ar nstate=2 ylag=0;
   smooth out=mylib.&ds.inSmooth_k2To2_p0To0;
   decode out=mylib.&ds.inDecode_k2To2_p0To0;
   evaluate out=mylib.&ds.inEval_k2To2_p0To0;
run;

The number of observations and model information are shown in Output 20.1.2.

The estimates of the initial state probability vector (ISPV) are shown in Output 20.1.3. In this example, the hidden Markov chain is assumed to be stationary. Hence, the steady-state probability distribution (STPD, which is also called the equilibrium distribution or stationary distribution) of the Markov chain is equal to the ISPV.

The estimates of the transition probability matrix (TPM) are shown in Output 20.1.4.

The estimates of the intercept for each state are shown in Output 20.1.5.

The estimates of the variance for each state are shown in Output 20.1.6.

The estimates of all parameters are shown in Output 20.1.7, which displays columns for the parameter name, estimate value, standard error, t value, and p-value.

The log likelihood and information criteria are shown in Output 20.1.8.

The smoothed probabilities of states for each observation are output in the data table mylib.vwmiInSmooth_k2To2_p0To0. This example discusses one type of inference. The mean and variance (or standard error) of each state can be calculated through the weighted sample, where the smoothed probabilities are the weights. The following statements calculate and print the sample mean and standard error of each state, as shown in Output 20.1.9. In the case of an RS-AR(0) model that has no regressors other than the intercept (which is equivalent to the Gaussian HMM), the sample mean and standard error of each state are almost the same as, respectively, the intercept estimate and the square root of variance estimate of each state.

*print sample mean and standard error for each state
 according to the weighted sample,
 where the smoothed probabilities are the weight;
data &ds.inSmooth_k2to2_p0to0;
   set mylib.&ds.inSmooth_k2to2_p0to0;
run;
proc sort data=&ds.inSmooth_k2to2_p0to0; by date; run;
data &ds.inWeighted;
   merge &ds.inSmooth_k2to2_p0to0(in=a) &ds.in(in=b);
   by date;
   if(a=b);
run;
%macro printMeanAndStdErrOfAState(k);
   data &ds.inMeanStdErr;
      set &ds.inWeighted end=eof;
      array cWeights{&k.} _TEMPORARY_;
      array cMeans[&k.] mean1-mean&k.;
      array cVars[&k.] sd1-sd&k.;
      array weights[&k.] state1-state&k.;
      retain cWeights cMeans cVars;
      if(_N_=1) then do;
         do i = 1 to &k.;
            cWeights[i] = weights[i];
            cMeans[i] = weights[i]*returnw;
            cVars[i] = weights[i]*returnw*returnw;
         end;
      end;
      else do;
         do i = 1 to &k.;
            cWeights[i] = cWeights[i] + weights[i];
            cMeans[i] = cMeans[i] + weights[i]*returnw;
            cVars[i] = cVars[i] + weights[i]*returnw*returnw;
         end;
      end;
      if eof then do;
         do i = 1 to &k.;
            if(cWeights[i]>0) then do;
               cMeans[i] = cMeans[i] / cWeights[i];
               cVars[i] = cVars[i] / cWeights[i] - cMeans[i]*cMeans[i];
               if(cVars[i]>=0) then cVars[i] = sqrt(cVars[i]);
            end;
            else do;
               cMeans[i] = .;
               cVars[i] = .;
            end;
         end;
         %do i = 1 %to &k.;
            call symputx("mu&i.",cMeans[&i.],'G');
            call symputx("sigma&i.",cVars[&i.],'G');
            state = &i.;  mean = mean&i.; stdErr = sd&i.;
            output;
         %end;
      end;
      label state="State" mean="Sample Mean"
            stdErr="Sample Standard Error";
      keep state mean stdErr;
   run;
   proc print data=&ds.inMeanStdErr noobs label; run;
%mend;
%printMeanAndStdErrOfAState(2);

The following statements plot each state as a Gaussian kernel, which shows how each state looks. As shown in Output 20.1.10, it is obvious that state 1 should represent the bull market state (positive expected return and low risk) and state 2 the bear market state (negative expected return and high risk).

%macro plotState(myData,myColumn,ks,ke);
   proc sgplot data=&myData.;
      refline 0 / axis=x lineattrs=(thickness=3);
      histogram &myColumn. / nbins=200;
      %do i = &ks. %to &ke.;
         density &myColumn. / type=normal(mu=&&mu&i. sigma=&&sigma&i.)
            name="state&i."
            legendlabel="State &i.";
      %end;
      keylegend %do i = &ks. %to &ke.; "state&i." %end;;
      xaxis min=-20 max=20 ;
      yaxis min=0 max=15;
   run;
%mend;
%plotState(&ds.In,returnw,1,2);

Output 20.1.10: Gaussian Distribution of Each State

The decoded path of states, which provide a historical view of what happened, given all the available data, is output in the data table mylib.vwmiInDecode_k2To2_p0To0. The decoded path assigns each time point a market state. It seems that the hidden market states have become observable. The following statements plot the decoded market states, as shown in Output 20.1.11:

* plot decoded results;
data &ds.InDecode_k2to2_p0to0; set mylib.&ds.InDecode_k2to2_p0to0; run;
proc sort data=&ds.InDecode_k2to2_p0to0; by date; run;
data decodeMaketStates;
   set &ds.In;
   set &ds.InDecode_k2to2_p0to0;
   array returnws returnw1-returnw2;
   array prices price1-price2;
   if(state~=.) then do;
      if(state=1) then do;
         returnws[1] = returnw; prices[1] = price;
         marketState="Bull Market";
      end;
      if(state=2) then do;
         returnws[2] = returnw; prices[2] = price;
         marketState="Bear Market";
      end;
   end;
   label marketState="Market State";
run;
proc sgplot data=decodeMaketStates;
   styleattrs datacontrastcolors=(green red);
   scatter x=date y=marketState / group=marketState;
   yaxis label="Market State" offsetmin=0.1 offsetmax=0.1 reverse;
run;

Output 20.1.11: Decoded Market States

You can incorporate the decoded market states with other variables of interest. The following statements plot the decoded weekly returns, as shown in Output 20.1.12:

proc sgplot data=decodeMaketStates;
   series x=date y=returnw1 / break lineattrs=(color=green)
                              legendlabel="Bull Market";
   series x=date y=returnw2 / break lineattrs=(color=red)
                              legendlabel="Bear Market";
   refline 0 / axis=y lineattrs=(color=gray thickness=1);
   yaxis max=20 min=-20 label="Weekly Return";
run;

Output 20.1.12: Decoded Weekly Returns

The following statements plot the decoded stock prices, as shown in Output 20.1.13:

proc sgplot data=decodeMaketStates;
   series x=date y=price1 / break lineattrs=(color=green)
                            legendlabel="Bull Market";
   series x=date y=price2 / break lineattrs=(color=red)
                            legendlabel="Bear Market";
   yaxis label="Market Index";
run;

Output 20.1.13: Decoded Stock Prices

The following statements plot the log likelihood of each observation conditional on its past; the plot is shown in Output 20.1.14.

* plot conditional log likelihood;
%macro plotEvaluation(inds, outds, qConditional);
   data &outds.;
      set &inds.;
      format date MMDDYY10.;
   run;
   %if (&qConditional.=1) %then %do;
      proc sort data=&outds.;
         by date;
      run;
      data &outds.;
         set &outds.;
         if(_N_=1) then
            condLogLikelihood = logLikelihood;
         else
            condLogLikelihood = logLikelihood - lag(logLikelihood);
      run;
      proc sgplot data=&outds.;
         series x=date y=condLogLikelihood;
         yaxis label='Conditional Log Likelihood' min=-15 max=0;
      run;
   %end;
   %else %do;
      proc sgplot data=&outds.;
         series x=date y=logLikelihood;
         yaxis label='Log Likelihood';
      run;
   %end;
%mend plotEvaluation;

%plotEvaluation(mylib.&ds.inEval_k2to2_p0to0, &ds.inEval_k2to2_p0to0, 1);

Output 20.1.14: Conditional Log Likelihood

As shown in Output 20.1.14, some observations are explained well by the model, but there are still many exceptional observations whose conditional log likelihoods are much smaller than others. This means that there might be model misspecification, which can be addressed by solving the model selection problem.

To model the weekly returns, how many states should be considered? How many AR lags should a state include? These are common questions in the model selection process. In this example, the best model is selected using Akaike’s information criterion (AIC): the smaller the AIC, the better the model.

The following macro estimates the model or models according to the arguments:

***************************************************************************;
* MODEL SELECTION                                                          ;
***************************************************************************;

%macro estimateRSAR(myds, inEstDs, kStart, kEnd, pStart, pEnd, method,
                    maxiter, qMultiStart);
   proc hmm data=&myds. labelSwitch=(sort=none)
      outstat=&myds.&method.Stat_k&kStart.To&kEnd._p&pStart.To&pEnd.;
      id time=date;
      model returnw / type=ar nstate=&kStart.:&kEnd. ylag=&pStart.:&pEnd.
         method=&method.;
      optimize printLevel=3 printIterFreq=1 algorithm=interiorpoint
         maxiter=&maxiter. Multistart=&qMultiStart.;
      score outmodel=&myds&method.Model_k&kStart.To&kEnd._p&pStart.To&pEnd.;
      learn out=&myds.&method.Learn_k&kStart.To&kEnd._p&pStart.To&pEnd.
         %if %length(&inEstDs.)>0 %then %do; in=&inEstDs. %end;
         ;
      display / excludeall;
   run;
%mend;

First, the following statements estimate two- to nine-state RS-AR(0) models. The MAP method is applied to bound and smooth the objective functions. The parameter estimates are saved for further use as the initial values for new estimations. For the best results, multistart mode is strongly recommended, although it might be time consuming.

* (1) MAP + MULTISTART for AR(0);
* be aware that the following macro might take tens of hours to finish;
* uncomment it to run;
* %estimateRSAR(myds=mylib.&ds.In, inEstDs=,
   kStart=2, kEnd=9, pStart=0, pEnd=0,
   method=MAP, maxiter=128, qMultiStart=1);

Then, the following statements estimate two- to nine-state RS-AR(0) to AR(5) models. The ML method is used because the AIC is based on the log likelihood. The parameter estimates of RS AR(0) models from previous MAP method estimations are used to ensure that the estimates in the ML method are reasonable. The k-state RS AR(p) models, , use the corresponding k-state RS AR() models’ parameter estimates as initial values. This is the default implementation in the HMM procedure. The multistart mode is not needed any longer because the initial values are good.

* (2) ML + initial values estimated from MAP's AR(0) for AR(0)
      and initial values estimated from AR(p-1) for AR(p), p>0;
* be aware that the following macro might take tens of minutes to finish;
* uncomment it to run;
* %estimateRSAR(myds=mylib.&ds.In, inEstDs=mylib.&ds.InMAPLearn_k2to9_p0to0,
   kStart=2, kEnd=9, pStart=0, pEnd=5,
   method=ML, maxiter=128, qMultiStart=0);

Estimating all 48 of the two- to nine-state RS AR(0) to AR(5) models produces AIC values. The following statements print the AIC values, which are shown in Output 20.1.15:

* print fit statistics in matrix form;
data &ds.InMLStat_k2to9_p0to5; set mylib.&ds.InMLStat_k2to9_p0to5; run;
proc sort data=&ds.InMLStat_k2to9_p0to5; by modelIndex; run;
%macro printFitStat(fitStat);
   data &ds.In&fitStat.;
      set &ds.InMLStat_k2to9_p0to5 end=eof;
      array &fitStat.s{8,6} _TEMPORARY_;
      array yLags{6} yLag0-yLag5;
      &fitStat.s[nState-1,yLag+1] = &fitStat.;
      if eof then do;
         do i = 1 to 8;
            nState = i+1;
            do j = 1 to 6; yLags(j) = &fitStat.s[i,j]; end;
            output;
         end;
      end;
      keep nState yLag0-yLag5;
   run;
   proc print data=&ds.In&fitStat. label noobs
      style(header)={textalign=center};
      var nState yLag0-yLag5;
      label nState='k' yLag0='p = 0' yLag1='p = 1' yLag2='p = 2'
         yLag3='p = 3' yLag4='p = 4' yLag5='p = 5';
   run;
%mend;
%printFitStat(AIC);

* Uncomment the following macros if you want to see other fit statistics;
* %printFitStat(logLikelihood);
* %printFitStat(AICC);
* %printFitStat(BIC);
* %printFitStat(HQC);

Output 20.1.15 indicates that the smallest AIC corresponds to the eight-state RS AR(4) model.

The following statements estimate the AIC-selected eight-state RS AR(4) model:

***************************************************************************;
* THE SELECTED MODEL AND IN-SAMPLE INFERENCE                               ;
***************************************************************************;

* select eight-state RS-AR(4) model according to AIC or AICC;
* adjust data to be the same sample size for AR(5);
proc hmm data=mylib.&ds.In(where=(date>='30JAN1926'd));
   id time=date;
   model returnw / type=ar nstate=8 ylag=4;
   optimize printLevel=3 algorithm=interiorpoint maxiter=0;
   score outmodel=mylib.&ds.InMLModel_k8To8_p4To4;
   learn out=mylib.&ds.InMLLearn_k8To8_p4To4
         in=mylib.&ds.InMLLearn_k2to9_p0to5;
   smooth out=mylib.&ds.inSmooth_k8to8_p4to4;
   decode out=mylib.&ds.inDecode_k8to8_p4to4;
   evaluate out=mylib.&ds.inEval_k8to8_p4to4;
   display tpm;
run;

The following statement plots the conditional log likelihoods. Compared to the bistate RS-AR(0) model (as shown in Output 20.1.14), the eight-state RS-AR(0) model (as shown in Output 20.1.16) has much higher conditional log likelihoods, which means a much better in-sample goodness of fit.

%plotEvaluation(mylib.&ds.inEval_k8to8_p4to4, &ds.inEval_k8to8_p4to4, 1);

Output 20.1.16: Conditional Log Likelihoods

The selected model has 112 parameters. 48 of which are related to observation distributions and are often called observation parameters. The following statements print the observation parameters in Output 20.1.17. The analysis of each AR(4) process for each state is beyond this example.

* print observation parameters;
data &ds.InObsParm;
   set mylib.&ds.InMLLearn_k8To8_p4To4(where=(type='EST' and index>=17));
   if(index=17) then parameter='Constant';
   if(index=19) then parameter='AR(1)';
   if(index=21) then parameter='AR(2)';
   if(index=23) then parameter='AR(3)';
   if(index=25) then parameter='AR(4)';
   if(index=27) then parameter='Variance';
run;
proc sort data=&ds.InObsParm; by index; run;
proc print data=&ds.InObsParm noobs label;
   var parameter state1-state8;
   label state1='State 1' state2='State 2' state3='State 3'
         state4='State 4' state5='State 5' state6='State 6'
         state7='State 7' state8='State 8';
run;

The following statements calculate and print the sample mean and standard error for each state in Output 20.1.18 and then plot the corresponding Gaussian kernels in Output 20.1.19:

data &ds.inSmooth_k8to8_p4to4;
   set mylib.&ds.inSmooth_k8to8_p4to4(where=(date>='25FEB1926'd));
run;
proc sort data=&ds.inSmooth_k8to8_p4to4; by date; run;
data &ds.inTruncated; set &ds.in(where=(date>='25FEB1926'd)); run;
data &ds.inWeighted;
   merge &ds.inSmooth_k8to8_p4to4(in=a) &ds.inTruncated(in=b);
   by date;
   if(a=b);
run;
%printMeanAndStdErrOfAState(8);

%plotState(&ds.In,returnw,1,8);

Output 20.1.19: Gaussian Kernels for Eight States

As shown in Output 20.1.19, states 1, 2, and 4 seem like a bull market; they have positive expected returns with low risk. The remaining states seem like a bear market; they have negative expected returns or very high risk (or both). However, the market states and their corresponding states are not that simple. The transition probability matrix of the eight-state AR(4) model is shown in Output 20.1.20.

The eight states seem to be categorized into three groups:

The transition probability matrix tells the more detailed story:

Hence, according to the eight-state RS-AR(4) model, there is one more market state besides the well-known bear market (the first group) and bull market (the third group): the transition market state (the second group). Figure 107 shows the relationship between these market states.

Figure 107: Market States and Their Transitions

The following statement plots three states that correspond to the bull market state, which are shown in Output 20.1.21:

%plotState(&ds.In,returnw,1,3);

Output 20.1.21: Three Gaussian Kernels for the Bull Market State

The following statement plots three states that correspond to the transition market state, which are shown in Output 20.1.22:

%plotState(&ds.In,returnw,4,6);

Output 20.1.22: Three Gaussian Kernels for the Transition Market State

The following statement plots two states that correspond to the bear market state, which are shown in Output 20.1.23:

%plotState(&ds.In,returnw,7,8);

Output 20.1.23: Two Gaussian Kernels for the Bear Market State

The following statements plot the decoded path of states according to the eight-state RS-AR(4) model; the result is shown in Output 20.1.24.

* plot decoded results;
data &ds.InDecode_k8to8_p4to4; set mylib.&ds.InDecode_k8to8_p4to4; run;
proc sort data=&ds.InDecode_k8to8_p4to4; by date; run;
data decodeMaketStates;
   set &ds.In(where=(date>='30JAN1926'd));
   set &ds.InDecode_k8to8_p4to4;
   array returnws returnw1-returnw3;
   array prices price1-price3;
   if(state~=.) then do;
      if(state=1 | state=2 | state=3) then do;
         returnws[1] = returnw; prices[1] = price;
         marketState="Bull Market      ";
      end;
      if(state=4 | state=5 | state=6) then do;
         returnws[2] = returnw; prices[2] = price;
         marketState="Transition Market";
      end;
      if(state=7 | state=8) then do;
         returnws[3] = returnw; prices[3] = price;
         marketState="Bear Market      ";
      end;
   end;
   label marketState="Market State";
run;
proc sgplot data=decodeMaketStates;
   styleattrs datacontrastcolors=(red yellow green);
   scatter x=date y=marketState / group=marketState nomissinggroup;
   yaxis label="Market State" offsetmin=0.1 offsetmax=0.1;
   keylegend / sortorder=reverseauto;
run;

Output 20.1.24: Decoded Market States

The following statements plot the decoded weekly returns, which are shown in Output 20.1.25:

proc sgplot data=decodeMaketStates;
   series x=date y=returnw1 / break lineattrs=(color=green)
                              legendlabel="Bull Market";
   series x=date y=returnw2 / break lineattrs=(color=yellow)
                              legendlabel="Transition Market";
   series x=date y=returnw3 / break lineattrs=(color=red)
                              legendlabel="Bear Market";
   refline 0 / axis=y lineattrs=(color=gray thickness=1);
   yaxis max=20 min=-20 label="Weekly Return";
run;

Output 20.1.25: Decoded Weekly Returns

The following statements plot the decoded stock prices, which are shown in Output 20.1.26:

proc sgplot data=decodeMaketStates;
   series x=date y=price1 / break lineattrs=(color=green)
                            legendlabel="Bull Market";
   series x=date y=price2 / break lineattrs=(color=yellow)
                            legendlabel="Transition Market";
   series x=date y=price3 / break lineattrs=(color=red)
                            legendlabel="Bear Market";
   yaxis label="Market Index";
run;

Output 20.1.26: Decoded Stock Prices

In the model selection process, the AIC, but not the in-sample goodness of fit (that is, the log likelihood), is used as the criterion for selecting the best model because a bias-variance trade-off exists. In this section, three types of out-of-sample forecasts are tested to check whether the selected model has a good forecast performance.

* print Predictive Performance of VaR Forecasts;
%let alpha1 = 0.80;
%let alpha2 = 0.90;
%let alpha3 = 0.98;
%macro forecastVaR(myLib, dsData, dsModel, timeId, varReturn,
   k, p, iStart, iEnd, oosStart, oosEnd,
   dsForecastPrefix);
   %do i = &iStart. %to &iEnd.;
      proc hmm data=&myLib..&dsData.;
         score inmodel=&myLib..&dsModel.;
         forecast out=&myLib..&dsForecastPrefix.&i. alpha=&&alpha&i. online;
      run;
      data &dsForecastPrefix.&i.;
         set &myLib..&dsForecastPrefix.&i.;
      run;
      proc sort data=&dsForecastPrefix.&i.; by &timeId.; run;
      data &dsForecastPrefix.&i.;
         set &dsForecastPrefix.&i. (FIRSTOBS=&oosStart. OBS=&oosEnd.
            keep=&timeId. &varReturn._Q1
            rename=(&varReturn._Q1=&varReturn._Q1_&k.&p.&i.));
         time=_N_; nState = &k; yLag = &p;
      run;
   %end;
%mend forecastVaR;

%macro evaluateVaR(dsData,dsForecastPrefix,varReturn,k,p,iStart,iEnd,n,dsOut);
   data &dsData.out;
      set &dsData.;
      if (date>&cutDate.) then output;
   run;
   data forecastData;
      set &dsData.out;
      time=_N_;
   run;
   data forecastData;
      merge &dsForecastPrefix.: forecastData;
      by time;
   run;
   data &dsOut.;
      set forecastData;
      retain %do i = &iStart. %to &iEnd.; cviol&k.&p.&i. 0 %end; ;
      retain %do i = &iStart. %to &iEnd.; c&varReturn.&k.&p.&i. 0 %end; ;
      %do i = &iStart. %to &iEnd.;
         if &varReturn. le &varReturn._Q1_&k.&p.&i. then do;
            cviol&k.&p.&i.=cviol&k.&p.&i.+1;
         end;
         c&varReturn.&k.&p.&i. = c&varReturn.&k.&p.&i.
                                  + &varReturn._Q1_&k.&p.&i.;
      %end;
      if _N_ = &n. then do;
         %do i = &iStart. %to &iEnd.;
            Norminal = (1-&&alpha&i.)/2;
            Viol = cviol&k.&p.&i. / &n.;
            LR = 2*(cviol&k.&p.&i.*log(Viol)+(&n.-cviol&k.&p.&i.)*log(1-Viol)
               -(cviol&k.&p.&i.*log(Norminal)+(&n.-cviol&k.&p.&i.)
                *log(1-Norminal)));
            pValue = 1 - cdf("CHISQUARE", LR, 1);
            meanVaR = c&varReturn.&k.&p.&i. / &n.;
            output;
         %end;
      end;
      label Norminal='Target Prob.' Viol='Violation Ratio' LR='LR Stat.'
            pValue='Pr > ChiSq' meanVaR='Avg. of VaR' nState='Number of States'
            yLag='AR Lag';
      keep nState yLag Norminal Viol LR pValue meanVaR;
   run;
   proc print data=&dsOut. noobs label;  format Viol LR pValue meanVaR 6.4; run;
%mend evaluateVaR;

%macro VaR(k,p);
   %forecastVaR(myLib=mylib, dsData=&ds.,
      dsModel=&ds.inMLModel_k&k.To&k._p&p.To&p.,
      timeId=date, varReturn=returnw,
      k=&k., p=&p., iStart=1, iEnd=3,
      oosStart=%eval(3999-(5-&p.)), oosEnd=%eval(4854-1-(5-&p.)),
      dsForecastPrefix=&ds.Forecastk&k._p&p.);
   %evaluateVaR(dsData=&ds.,dsForecastPrefix=&ds.Forecastk&k._p&p.,
      varReturn=returnw,k=&k.,p=&p.,iStart=1,iEnd=3,n=855,
      dsOut=VaR_outputk&k._p&p.);
%mend VaR;

%VaR(8,4);

* print in-sample and out-of-sample average weekly log likelihood;
proc hmm data=mylib.&ds.
         outstat=mylib.&ds.Stat_k2to9_p0to5;
   score inmodel=mylib.&ds.inmlmodel_k2to9_p0to5;
   evaluate out=mylib.&ds.Eval_k2to9_p0to5;
run;
data inLL;
   set mylib.&ds.Eval_k2to9_p0to5(where=(date='26DEC2000'd)
      rename=(logLikelihood=isLL));
run;
data fullLL;
   set mylib.&ds.Eval_k2to9_p0to5(where=(date='22DEC2017'd)
      rename=(logLikelihood=fullLL));
run;
data oosLL;
   merge inLL(in=a) fullLL(in=b);
   by modelIndex;
   oosLL = (fullLL - isLL)/855;
   if(a=b);
   drop date isLL fullLL;
run;

data inLL;
   set mylib.&ds.InMLStat_k2to9_p0to5;
   isLL = logLikelihood / (3999-5);
   keep modelIndex nState yLag isLL;
run;
proc sort data=inLL; by modelIndex; run;

data avgLL;
   merge inLL(in=a) oosLL(in=b);
   by modelIndex;
   if(a=b and (modelIndex=1 or modelIndex=41 or modelIndex=48));
   drop modelIndex;
run;

proc print data=avgLL noobs label;
   label nState='k' yLag='p' isLL='In-Sample' oosLL='Out-Of-Sample';
run;

Trading Strategies and Portfolio Values

As mentioned at the beginning of this example, if you knew it would be a bull market, you might choose to buy the stock; if you knew it would be a bear market, you might choose to sell the stock. In this section, the market states are forecasted and then several trading strategies are compared.

The following statements forecast and plot the probability distribution of market states in Output 20.1.29. The forecast of market states is in probabilities, which is similar to a weather forecast. For example, a forecast of 60% chance of bear market, 39% chance of transition market, and 1% chance of bull market is comparable to a forecast of 60% chance of raining, 39% chance of being cloudy, and 1% chance of being sunny.

proc hmm data=mylib.&ds.;
   score inmodel=mylib.&ds.InMLModel_k8To8_p4To4;
   forecast out=mylib.&ds.Forecast_k8to8_p4to4 online;
   display / excludeall;
run;

%let forecastStartDate = '26DEC2000'd;
%let forecastEndDate = '22DEC2017'd;
data ForecastProb;
   set mylib.&ds.Forecast_k8to8_p4to4
          (where=(date>=&forecastStartDate. and date<=&forecastEndDate.));
run;
proc sort data=ForecastProb; by date; run;
data ForecastProb;
   set ForecastProb(keep=state1 - state8) end=last;
   set &ds.(where=(date>=&forecastStartDate. and date<=&forecastEndDate.)
            keep=date price);
   * probability distribution of market states;
   bullProb = state1+state2+state3;
   transitionProb = state4+state5+state6;
   bearProb = state7+state8;
   bearTransitionProb = bearProb + transitionProb;
run;
* plot forecasted prob. dist. of market states and trading strategies;
%macro plotForecastAndTradingStrategies(ds, bearThreshold, bullThreshold,
                                        oneMinusBullThreshold);
   proc sgplot data=&ds.;
      band x=date lower=0 upper=bearProb /
         fillattrs=(color=red) name="bear" legendlabel="Bear Market";
      band x=date lower=bearProb upper=bearTransitionProb /
         fillattrs=(color=yellow)
         name="transition" legendlabel="Transition Market";
      band x=date lower=bearTransitionProb upper=1 /
         fillattrs=(color=green) name="bull" legendlabel="Bull Market";
      series x=date y=price / y2axis lineattrs=(color=blue)
         name="index" legendlabel="Market Index";
      %if %length(&bearThreshold.)>0 %then %do;
         refline &bearThreshold. / axis=y lineattrs=(color=red thickness=3)
            label="Bear Threshold=&bearThreshold." labelloc=inside
            labelattrs=(weight=bold);
      %end;
      %if %length(&bearThreshold.)>0 %then %do;
         refline &oneMinusBullThreshold. / axis=y
            lineattrs=(color=green thickness=3)
            label="Bull Threshold=&bullThreshold." labelloc=inside
            labelattrs=(weight=bold);
      %end;
      xaxis values=("01JAN2001"d to "31DEC2017"d by year) display=(nolabel)
         tickvalueformat=year4. offsetmin=0.03 offsetmax=0.15;
      yaxis label="Forecasted Prob. Dist. of Market States";
      y2axis label="Market Index";
      keylegend "bull" "transition" "bear" "index"/
         location=outside position=bottom;
   run;
%mend;

%plotForecastAndTradingStrategies(ForecastProb);

Output 20.1.29: Forecasted Probabilities of Market States

In this example, four simple trading strategies are considered:

• Risk-free strategy: hold the bond for the entire time period. This strategy does not depend on any forecast, and it is the first benchmark strategy.

• Market strategy: hold the stock for the entire time period. This strategy does not depend on any forecast, either. It reflects the market performance, and it is the second benchmark strategy.

• Bear strategy: sell the stock if the chance of a bear market in the next week is no less than 30%; otherwise, buy the stock.

• Bull strategy: buy the stock if the chance of a bull market in the next week is no less than 10%; otherwise, sell the stock.

The thresholds in bear and bull strategies, 30% and 10% respectively, are ad hoc, which reflect the investor’s risk attitude. The trading strategy of switching between bond and stock can also be analogous to the decision of traveling with or without umbrella. The risk-free strategy is an analogous to a decision to always travel with an umbrella; the market strategy is to analogous to never traveling with umbrella; the bear strategy is analogous to traveling with an umbrella if the chance of raining tomorrow is no less than 30%; and the bull strategy is analogous to traveling without an umbrella if the chance of being sunny tomorrow is no less than 10%. The following statements plot the thresholds for the bear and bull trading strategies:

%let bearThreshold=0.30;
%let bullThreshold=0.10;
%let oneMinusBullThreshold=0.90;
%plotForecastAndTradingStrategies(ForecastProb, &bearThreshold.,
                                  &bullThreshold., &oneMinusBullThreshold.);

Output 20.1.30: Thresholds for Bull and Bear Trading Strategies

The following statements calculate the portfolio values and plot the wealth curves of four trading strategies:

%let strategy1=RiskFree;
%let strategy2=Market;
%let strategy3=Bear;
%let strategy4=Bull;
%let strategy1Label=Risk Free;
%let strategy2Label=Market;
%let strategy3Label=Bear;
%let strategy4Label=Bull;
* assume that 1 year = 252 trading days = 50.4 weeks;
* weekly risk free rate = exp(ln(1+APY)/50.4)-1;
%let riskFreeRate = 0.0002661; * APY=1.35%;
%let initialMoney=1;

%macro port_val_cal(stockPosition, stockPrice, riskFreeRate, strategy,
                    initialMoney, rowNo);
   if (&rowNo. = 1) then do;
      &strategy.CashValue = &initialMoney.-&initialMoney.*&stockPosition.;
      &strategy. = &initialMoney.;
      &strategy.StockAmount = &strategy.*&stockPosition./&stockPrice.;
      &strategy.CurrentCashValue = &strategy.CashValue;
  end;
  else do;
      &strategy.CurrentCashValue = &strategy.CashValue*(1+&riskFreeRate.);
      &strategy. = &strategy.CurrentCashValue
                  + &strategy.StockAmount*&stockPrice.;
      &strategy.StockAmount = &strategy.*&stockPosition./&stockPrice.;
      &strategy.CashValue = &strategy.*(1-&stockPosition.);
  end;
  &strategy.ExcessReturn = &strategy./lag(&strategy.)-1-&riskFreeRate.;
  &strategy.ExcessReturn2 = &strategy.ExcessReturn**2;
%mend port_val_cal;

%macro port_val_gen(ForecastProb, bearThreshold, bullThreadhold, riskFreeRate,
                    initialMoney, startDate, endDate);
   data PortValPath;
      set &ForecastProb. end=last;
      array stockPos[4] %do i = 1 %to 4; &&strategy&i..StockPos %end; ;
      retain
         %do i = 1 %to 4;
            &&strategy&i..CashValue 0 &&strategy&i..StockAmount 0
            &&strategy&i..CurrentCashValue 0 &&strategy&i..CumExcessReturn 0
            &&strategy&i..CumExcessReturn2 0
         %end; ;
      * stock position for each strategy;
      stockPos[1] = 0;
      stockPos[2] = 1;
      if(bearProb>=&bearThreshold.) then stockPos[3]=0;
      else stockPos[3]=1;
      if(bullProb>=&bullThreshold.) then stockPos[4]=1;
      else stockPos[4]=0;
      * wealth for each strategy;
      %do i = 1 %to 4;
         %port_val_cal(stockPosition=stockPos[&i.], stockPrice=price,
                       riskFreeRate=&riskFreeRate., strategy=&&strategy&i..,
                       initialMoney=&initialMoney., rowNo=_N_);
         &&strategy&i..CumExcessReturn+&&strategy&i..ExcessReturn;
         &&strategy&i..CumExcessReturn2+&&strategy&i..ExcessReturn2;
      %end;
      * Sharpe ratio and portfolio final wealth;
      if last then do;
         &strategy1.LastPorVal = &strategy1.;
         %do i = 2 %to 4;
            &&strategy&i..LastPorVal = &&strategy&i..;
            &&strategy&i..ExcessReturnMean =
               &&strategy&i..CumExcessReturn/(_N_-1);
            &&strategy&i..ExcessReturnStd =
               sqrt(&&strategy&i..CumExcessReturn2/(_N_-1)
                  - &&strategy&i..ExcessReturnMean**2);
            &&strategy&i..SharpeRatio = round(&&strategy&i..ExcessReturnMean
               /&&strategy&i..ExcessReturnStd, 0.0001);
            call symputx("&&strategy&i..SharpeRatio",
               &&strategy&i..SharpeRatio, 'G');
         %end;
      end;
   run;
%mend;

* plot portfolio values of different trading strategies;
%let color1=black;
%let color2=blue;
%let color3=red;
%let color4=green;
%macro perf_comp_plot_gen(nStrategies);
   proc sgplot data=PortValPath cycleattrs;
      %do i = 1 %to &nStrategies.;
         %let strategy=&&strategy&i..;
         series x=date y=&&strategy&i.. / name="&strategy."
            lineattrs=(color=&&color&i..)
            %if (&i.=1) %then %do;
               legendlabel="&&strategy&i.Label.";
            %end;
            %else %do;
               legendlabel="&&strategy&i.Label. (SR=&&&strategy.SharpeRatio.)";
            %end;
         scatter x=date y=&&strategy&i..LastPorVal / datalabelpos=topright
            markerattrs=(color=&&color&i..)
            datalabel=&&strategy&i..LastPorVal;
      %end;
      xaxis values=("01JAN2001"d to "31DEC2017"d by year) display=(nolabel)
         tickvalueformat=year4. offsetmin=0.03 offsetmax=0.15;
      yaxis max=2.6 min=0.5 label="Wealth Curve";
      keylegend %do i = 1 %to &nStrategies.; "&&strategy&i.." %end; /
         location=outside position=bottom;
   run;
%mend perf_comp_plot_gen;

%port_val_gen(ForecastProb=ForecastProb,
              bearThreshold=&bearThreshold., bullThreadhold=&bullThreshold.,
              riskFreeRate=&riskFreeRate., initialMoney=&initialMoney.,
              startDate=&forecastStartDate., endDate=&forecastEndDate.);
%perf_comp_plot_gen(nStrategies=4);

Output 20.1.31: Wealth Curves of Four Trading Strategies

As shown in Output 20.1.31, thanks to the forecast of market states by the eight-state RS-AR(4) model, both bear and bull strategies beat the market: compared to the market strategy, they increase the Sharpe ratio more than 50% and get higher portfolio values at the end of out-of-sample period.