HMM Procedure

Let , denote a p-dimensional time series vector of random variables. Conditional on a latent variable , follows a normal distribution,

where and are mean and covariance parameters whose values depend on the variable . The variable , the so-called state, 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 variable is unobservable and follows a Markov chain and because follows a Gaussian distribution conditional on , the model is called the Gaussian hidden Markov model (Gaussian HMM). The Gaussian HMM is sometimes described as , where , consisting of parameters that define the state-dependent distribution of observable variables.

Consider a bivariate Gaussian HMM that has the following parameter values:

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

%let pi1 = 0.5;
%let a11 = 0.001;
%let a22 = 0.001;
%let mu1_1 = 1;
%let mu1_2 = 1;
%let sigma1_11 = 1;
%let sigma1_21 = 0;
%let sigma1_22 = 1;
%let mu2_1 = -1;
%let mu2_2 = -1;
%let sigma2_11 = 1;
%let sigma2_21 = 0;
%let sigma2_22 = 1;
%let T = 10000;
%let seed = 1234;

data DGPone;
   retain cd1_11 cd1_21 cd1_22 cd2_11 cd2_21 cd2_22;
   do t = 1 to &T.;
      if(t=1) then do;
         /* initial probability distribution */
         p = &pi1.;
         /* Cholesky decomposition of sigma1 */
         cd1_11 = sqrt(&sigma1_11.);
         cd1_21 = &sigma1_21./sqrt(&sigma1_11.);
         cd1_22 = sqrt(&sigma1_22.-&sigma1_21.*&sigma1_21./&sigma1_11.);
         /* Cholesky decomposition of sigma2 */
         cd2_11 = sqrt(&sigma2_11.);
         cd2_21 = &sigma2_21./sqrt(&sigma2_11.);
         cd2_22 = sqrt(&sigma2_22.-&sigma2_21.*&sigma2_21./&sigma2_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;
      e1 = normal(&seed.); e2 = normal(&seed.);
      if(s=1) then do;
         /* (x,y) ~ N(mu1, Sigma1) at state 1 */
         x = &mu1_1. + cd1_11*e1;
         y = &mu1_2. + cd1_21*e1+cd1_22*e2;
      end;
      else do;
         /* (x,y) ~ N(mu2, Sigma2) at state 2 */
         x = &mu2_1. + cd2_11*e1;
         y = &mu2_2. + cd2_21*e1+cd2_22*e2;
      end;
      output;
      lags = s;
   end;
run;

data One;
   set DGPone;
   keep t x y;
run;

In general, the data generating process is unknown. What can be seen is the data table One, but not the data table DGPone.

The following statements plot the simulated series x:

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

Figure 1 shows the plot of x, whose points seem to fall randomly around 0.

The following statements plot the simulated series y:

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

Figure 2 shows the plot of y, whose points also seem to fall randomly around 0.

Ignoring the sequential information, the following statements plot the observed points :

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

As shown in Figure 3, it seems that the points come from two different cluster centers, and many of them lie between two cluster centers.

It is very difficult for any cluster analysis methods that cannot use the sequential information to find the correct states for the points that lie between two cluster centers. The simplest cluster analysis method to classify the points, because of the symmetry shown in Figure 3, is to draw a line , classify all points above the line to one state, and classify all points below the line to another state.^[6] The following statements implement this simplest cluster analysis method and calculate the accuracy, based on the data table DGPone:

/* classify the data points by the line y = -x */
data Cluster;
   set One;
   if(y>-x) then state = 1;
   else state = 2;
   keep t state;
run;

data clusterCheck;
   merge DGPone(in=a) Cluster(in=b);
   by t;
   if(s=state) then correct = 1;
   else correct = 0;
   if(a=b);
   keep correct;
run;

data clusterAccuracy;
   set clusterCheck;
   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 accuracy = 1 - accuracy;
      output;
   end;
   keep accuracy;
run;

proc print data = clusterAccuracy noobs; run;

The accuracy of the cluster analysis, which is shown in Figure 4, is 92.02%. This implies that the cluster analysis cannot choose the correct states for about 800 data points!

What is the big difference that you get by using the sequential information through HMMs? Before estimating any model, make a copy of the data in the library:

   data mylib.One; set One; run;

The following statements estimate the two-state Gaussian HMM and then classify the data points through decoding; here decoding means to find the best possible state path for the observed data. The SORT=DESC(MU) option requests that the states be labeled in descending order of the mean parameter values; for more information about the label switching problem, see the section Label Switching Problem.

proc hmm data=mylib.One labelSwitch=(sort=desc(mu));
   id time=t;
   model x y / type=gaussian nstate=2;
   optimize printLevel=3 printIterFreq=1;
   decode out=mylib.oneDecode;
run;

The number of observations and model information are shown in Figure 5.

The default method, which is the maximum likelihood method, is applied; this is a nonlinear optimization problem. The initial parameter values and the objective function value at the start of the optimization are shown in Figure 6.

The details of the optimization algorithm are shown in Figure 7.

The iterations of the optimization process are shown in Figure 8.

The final parameter values and the objective function value at the end of the optimization are shown in Figure 9.

The estimates of the initial state probability vector (ISPV) are shown in Figure 10.

The estimates of the transition probability matrix (TPM) are shown in Figure 11.

The estimates of the mean vector of the Gaussian distribution for each state are shown in Figure 12.

The estimates of the covariance matrix of the Gaussian distribution for each state are shown in Figure 13.

The estimates of all parameters are shown in Figure 14, 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 Figure 15.

The accuracy of classification through the Gaussian HMM is calculated by the following statements:

data gHmmCheck;
   merge DGPone(in=a) mylib.oneDecode(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 gHmmAccuracy;
   set gHmmCheck;
   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 gHmmCheck;
   set gHmmCheck;
   if(&qFlipState.=1) then do;
      ds = -ds;
      correct = 1 - correct;
   end;
run;

proc print data = gHmmAccuracy noobs; run;

The accuracy of the Gaussian HMM, which is shown in Figure 16, is 99.99%. This means that the Gaussian HMM successfully takes advantage of the sequential information: only 1 data point cannot be correctly classified, compared to about 800 points that cannot be correctly classified by the cluster analysis method without the sequential information.

To show how well the data are classified, the previous statements label any incorrectly classified data points with negative state values. In this example, only 1 point is not correctly classified. You plot the classified data points by using the following statements. The plot is shown in Figure 17.

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