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, 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 regression parameters depend on the state, the follows different regressions in different regimes; hence, this type of model is called the regime-switching regression model.

Consider a univariate regime-switching regression model that has two regimes and two regressors (an intercept and an exogenous variable):

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 x_lb = -8;
%let x_ub = 0;
%let pi1 = 0.5;
%let a11 = 0.95;
%let a22 = 0.95;
%let const1_1 = 0;
%let xl1_0_1_1 = 1;
%let cov1_11 = 2.56;
%let const2_1 = 4;
%let xl2_0_1_1 = 1.5;
%let cov2_11 = 4;
%let T = 4000;
%let seed = 1234;

data rsregdgp;
   retain cd1_11 cd2_11;
   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;
      x = &x_lb. + (&x_ub. - &x_lb.)*uniform(&seed.);
      e = normal(&seed.);
      if(s=1) then do;
         /* y ~ N(beta1_0+beta1_1*x, Sigma1) at state 1 */
         y = &const1_1. + &xl1_0_1_1 * x + cd1_11*e;
      end;
      else do;
         /* y ~ N(beta2_0+beta2_1*x, Sigma2) at state 2 */
         y = &const2_1. + &xl2_0_1_1 * x + cd2_11*e;
      end;
      output;
      lags = s;
   end;
run;

data rsreg;
   set rsregdgp;
   keep t y x;
run;

The following statements create the scatter plot in Figure 43:

proc sgplot data=rsreg;
   scatter y=y x=x;
run;

As shown in Figure 43, when the time dependence is ignored, it is difficult—if not impossible—to distinguish the two regimes.

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

data mylib.rsreg; set rsreg; run;
proc hmm data=mylib.rsreg;
   id time=t;
   model y = x / type=reg nstate=2 method=ml;
   optimize algorithm=activeset maxiter=256 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 44, which displays columns for the 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 regression model is calculated by the following statements:

data rsregCheck;
   merge rsregdgp(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 x y s correct ds;
run;

%let qFlipState = -1;
data rsregAccuracy;
   set rsregCheck;
   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 rsregCheck;
   set rsregCheck;
   if(&qFlipState.=1) then do;
      ds = -ds;
      correct = 1 - correct;
   end;
run;

proc print data = rsregAccuracy noobs; run;

The accuracy of the regime-switching regression model, which is shown in Figure 45, is close to 90%. This means that the regime-switching regression model successfully distinguishes the two regimes.

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

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