HMM Procedure

Let , denote a p-dimensional time series vector of random variables. The can be modeled in the regression

where , the regressor is observable and contains the lagged dependent variables, the latent variable is the so-called state, and and are mean and covariance parameters whose values depend on the state . The variable follows the first-order Markov chain; that is,

where denotes the conditional probability. The range of is a finite set, . The transition probability from state i to state j is expressed as

The matrix is called the transition probability matrix (TPM). The last element in the model is the initial state probability vector (ISPV), , of the first state :

Because the lagged dependent variables are included in the regressors and the autoregression parameters depend on the state, the follows different autoregressions in different regimes; hence, this type of model is called the regime-switching autoregression model.

Consider a univariate regime-switching autoregression model that has two regimes:

The initial state probability vector and transition probability matrix are as follows:

The following statements simulate the time series from the previous model to provide test data for the HMM procedure:

%let pi1 = 0.5;
%let a11 = 0.95;
%let a22 = 0.95;
%let ar1_1_1_1 = 0.8;
%let cov1_11 = 2.56;
%let ar2_1_1_1 = -0.7;
%let cov2_11 = 4;
%let T = 1000;
%let seed = 1234;

data rsardgp;
   retain cd1_11 cd2_11;
   ylag = 0;
   do t = 1 to &T.;
      if(t=1) then do;
         /* initial probability distribution */
         p = &pi1.;
         /* Cholesky decomposition of COV1 */
         cd1_11 = sqrt(&cov1_11.);
         /* Cholesky decomposition of COV2 */
         cd2_11 = sqrt(&cov2_11.);
      end;
      else do;
         /* transition probability matrix */
         if(lags=1) then p = &a11.;
         else p = 1-&a22.;
      end;
      u = uniform(&seed.);
      if(u<=p) then s=1;
      else s = 2;
      e = normal(&seed.);
      if(s=1) then do;
         /* y ~ N(beta1*ylag, Sigma1) at state 1 */
         y = &ar1_1_1_1 * ylag + cd1_11*e;
      end;
      else do;
         /* y ~ N(beta2*ylag, Sigma2) at state 2 */
         y = &ar2_1_1_1 * ylag + cd2_11*e;
      end;
      output;
      lags = s;
      ylag = y;
   end;
run;

data rsar;
   set rsardgp;
   keep t y;
run;

The following statements create the series plot in Figure 47. It is difficult to tell directly from the plot that the series are from two different regimes.

proc sgplot data=rsar;
   series y=y x=t;
run;

The following code uploads the data table and estimates the regime-switching autoregression model:

data mylib.rsar; set rsar; run;
proc hmm data=mylib.rsar;
   id time=t;
   model y / type=ar noint ylag=1 nstate=2 method=ml;
   optimize algorithm=activeset printLevel=3 printIterFreq=1;
   learn out=mylib.mylearn;
   filter out=mylib.myfilter;
   smooth out=mylib.mysmooth;
   decode out=mylib.mydecode;
   evaluate out=mylib.myeval;
run;

These statements assume that your engine libref is named mylib, but you can substitute any appropriately defined engine libref.

The estimates of all parameters are shown in Figure 48, which displays columns for parameter name, estimate value, standard error, t value, and p-value. The parameter estimates are very close to the true parameter values that are used in the data generating process.

The accuracy of classification through the regime-switching autoregression model is calculated by the following statements:

data rsarCheck;
   merge rsardgp(in=a) mylib.mydecode(in=b);
   by t;
   if(s=state) then do; correct=1; ds = s; end;
   else do; correct=0; ds = -s; end;
   if(a=b);
   keep t y s correct ds;
run;

%let qFlipState = -1;
data rsarAccuracy;
   set rsarCheck;
   retain count 0 correctCount 0;
   count = count + 1;
   correctCount = correctCount + correct;
   if(count>&T.-0.5) then do;
      accuracy = correctCount / count;
      if(accuracy<0.5) then call symputx("qFlipState",1,'G');
      else call symputx("qFlipState",0,'G');
      if(accuracy<0.5) then accuracy = 1 - accuracy;
      output;
   end;
   keep accuracy;
run;

data rsarCheck;
   set rsarCheck;
   if(&qFlipState.=1) then do;
      ds = -ds;
      correct = 1 - correct;
   end;
run;

proc print data = rsarAccuracy noobs; run;

The accuracy of the regime-switching autoregression model, which is shown in Figure 49, is more than 92%. This means that the regime-switching autoregression model successfully distinguishes the two regimes.

The following statements create the scatter plot in Figure 50 to show how well the data points are classified:

proc sgplot data=rsarCheck;
   scatter y=y x=t / group=ds;
run;