HMM Procedure

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

where the variable is the (hidden) state; represents the Gaussian mixture distribution; M is the number of components at each state; and the , are mixture component probabilities (MCPs, also called mixture weights), which satisfy the basic requirement for weights ( and ). The and are the mean and covariance parameters for the jth Gaussian component at state , and their values depend on the variable . The state 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 variable is unobservable and follows a Markov chain and because follows a Gaussian mixture distribution conditional on , the model is called the Gaussian mixture hidden Markov model (GM HMM). The GM HMM is sometimes described as , where , and it consists of parameters that define the state-dependent distribution of observable variables.

The GM HMM for cross-sectional time series data is similar to the extension from Gaussian HMM to Gaussian HMM for cross-sectional time series data; for more information, see the section Gaussian Hidden Markov Model for Cross-Sectional Time Series Data.

Consider a bistate bi-component GM HMM that has the following parameter values. The ISPV and TPM are

At state 1, the two components have weights

and the mean and covariance parameters for two components are

At state 2, the two components have weights

and the mean and covariance parameters for two components are

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

%let a11 = 0.98;
%let a22 = 0.95;
%let c1_1 = 0.6;
%let c1_2 = 0.4;
%let c2_1 = 0.7;
%let c2_2 = 0.3;
%let mu1_1_1 = 1;
%let mu1_1_2 = 1;
%let mu1_2_1 = -3;
%let mu1_2_2 = -3;
%let sigma1_1_11 = 1;
%let sigma1_1_21 = 0;
%let sigma1_1_22 = 1;
%let sigma1_2_11 = 1.69;
%let sigma1_2_21 = 1;
%let sigma1_2_22 = 1.69;
%let mu2_1_1 = -1;
%let mu2_1_2 = -1;
%let mu2_2_1 = 3;
%let mu2_2_2 = 3;
%let sigma2_1_11 = 1;
%let sigma2_1_21 = 0;
%let sigma2_1_22 = 1;
%let sigma2_2_11 = 1.69;
%let sigma2_2_21 = 1;
%let sigma2_2_22 = 1.69;
%let seed = 1234;
%let T = 200000;

data gmDGP;
   retain cd1_1_11 cd1_1_21 cd1_1_22
          cd1_2_11 cd1_2_21 cd1_2_22
          cd2_1_11 cd2_1_21 cd2_1_22
          cd2_2_11 cd2_2_21 cd2_2_22;
   do t = 1 to &T.;
      if(t=1) then do;
         * initial probability distribution;
         p = (1-&a22.)/((1-&a11.)+(1-&a22.));;
         * Cholesky decomposition;
         cd1_1_11 = sqrt(&sigma1_1_11.);
         cd1_1_21 = &sigma1_1_21./cd1_1_11;
         cd1_1_22 = sqrt(&sigma1_1_22.-cd1_1_21*cd1_1_21);
         cd1_2_11 = sqrt(&sigma1_2_11.);
         cd1_2_21 = &sigma1_2_21./cd1_2_11;
         cd1_2_22 = sqrt(&sigma1_2_22.-cd1_2_21*cd1_2_21);
         cd2_1_11 = sqrt(&sigma2_1_11.);
         cd2_1_21 = &sigma2_1_21./cd2_1_11;
         cd2_1_22 = sqrt(&sigma2_1_22.-cd2_1_21*cd2_1_21);
         cd2_2_11 = sqrt(&sigma2_2_11.);
         cd2_2_21 = &sigma2_2_21./cd2_2_11;
         cd2_2_22 = sqrt(&sigma2_2_22.-cd2_2_21*cd2_2_21);
      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;
         * choose component;
         u = uniform(&seed.);
         if(u<=&c1_1.) then do;
            * (x,y) ~ N(mu, Sigma) at state 1, component 1;
            x = &mu1_1_1. + cd1_1_11*e1;
            y = &mu1_1_2. + cd1_1_21*e1+cd1_1_22*e2;
         end;
         else do;
            * (x,y) ~ N(mu, Sigma) at state 1, component 2;
            x = &mu1_2_1. + cd1_2_11*e1;
            y = &mu1_2_2. + cd1_2_21*e1+cd1_2_22*e2;
         end;
      end;
      else do;
         * choose component;
         u = uniform(&seed.);
         if(u<=&c2_1.) then do;
            * (x,y) ~ N(mu, Sigma) at state 2, component 1;
            x = &mu2_1_1. + cd2_1_11*e1;
            y = &mu2_1_2. + cd2_1_21*e1+cd2_1_22*e2;
         end;
         else do;
            * (x,y) ~ N(mu, Sigma) at state 2, component 2;
            x = &mu2_2_1. + cd2_2_11*e1;
            y = &mu2_2_2. + cd2_2_21*e1+cd2_2_22*e2;
         end;
      end;
      output;
      lags = s;
   end;
run;

data gmTrain;
   set gmDGP(where=(t<=&T./2));
   keep t x y;
run;

data gmTest;
   set gmDGP(where=(t>&T./2));
   keep t x y;
run;

In general, the data generating process is unknown. The data tables gmTrain and gmTest can be seen, but the data table gmDGP cannot be seen.

First, assume that the following statements plot the simulated series x that is contained in the data table gmTrain:

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

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

The following statements plot the simulated series y that is contained in the data set gmTrain:

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

Figure 23 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 that are contained in the data table gmTrain:

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

As shown in Figure 24, it seems that all points come from one bivariate Gaussian distribution that has a high positive correlation. However, according to the true DGP, some points are from state 1, in which there are two clusters, as shown in Figure 25:

proc sgplot data=gmDGP(where=(t<=&T./2 and s=1));
   scatter y=y x=x;
run;

Other points are from state 2 in which there are two clusters, as shown in Figure 26:

proc sgplot data=gmDGP(where=(t<=&T./2 and s=2));
   scatter y=y x=x;
run;

Without considering the sequential order, it is impossible to distinguish the two states.

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

   data mylib.gmTrain; set gmTrain; run;
   data mylib.gmTest; set gmTest; run;

The following statements estimate the bistate bi-component GM HMM. The type of model is specified in the TYPE= option in the MODEL statement. The number of states is specified in the NSTATE= option in the MODEL statement. The number of components for each state is specified in the NCOMPONENT= option in the MODEL statement. The OUTMODEL= option in the SCORE statement outputs the model information for future scoring on new data table. The LABELSWITCH=(SORT=NONE) option in the PROC HMM statement suppresses the relabeling of states and components; for more information about the label switching problem, see the section Label Switching Problem.

proc hmm data=mylib.gmTrain labelSwitch=(sort=none);
   id time=t;
   model x y / type=gaussianMixture nstate=2 ncomponent=2;
   optimize printLevel=3 printIterFreq=1;
   score outmodel=mylib.gmModel;
run;

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

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

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 28.

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

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

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

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

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

The estimates of the mixture component probabilities (MCP) are shown in Figure 34.

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

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

The estimates of all parameters are shown in Figure 37, which displays columns for the parameter name, estimate value, standard error, t value, and p-value. The state labels ("1" and "2" in this example) and the component labels ("1" and "2" in this example) might exchange randomly, compared to the true DGP. In the preceding estimation, the state labels are the same as the labels in the DGP. The component labels for state 1 are also the same as the ones in the DGP, but the component labels for state 2 do exchange! Compared to the parameter values in the DGP, the estimates are very accurate.

The log likelihood and information criteria are shown in Figure 38.

The estimates from the data table gmTrain are used to score the testing data table gmTest to check the predictability in the following statements:

proc hmm data=mylib.gmTest;
   score inmodel=mylib.gmModel;
   decode out=mylib.gmDecode;
run;

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

data gmHmmCheck;
   merge gmDGP(in=a where=(t>&T./2)) mylib.gmDecode(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 gmHmmAccuracy;
   set gmHmmCheck;
   retain count 0 correctCount 0;
   count = count + 1;
   correctCount = correctCount + correct;
   if(count>&T./2-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 gmHmmCheck;
   set gmHmmCheck;
   if(&qFlipState.=1) then do;
      ds = -ds;
      correct = 1 - correct;
   end;
run;

proc print data = gmHmmAccuracy noobs; run;

The accuracy of the GM HMM, which is shown in Figure 39, is about 98%. This means that the GM HMM successfully takes advantage of the sequential information.

Now, consider the GM HMM for cross-sectional time series data. Each section’s data follow the same GM HMM as described previously. The following statements generate the data table for analysis. There are 200 sections in the data table, and each section of time series consist of 1,000 observations.

%let a11 = 0.98;
%let a22 = 0.95;
%let c1_1 = 0.6;
%let c1_2 = 0.4;
%let c2_1 = 0.7;
%let c2_2 = 0.3;
%let mu1_1_1 = 1;
%let mu1_1_2 = 1;
%let mu1_2_1 = -3;
%let mu1_2_2 = -3;
%let sigma1_1_11 = 1;
%let sigma1_1_21 = 0;
%let sigma1_1_22 = 1;
%let sigma1_2_11 = 1.69;
%let sigma1_2_21 = 1;
%let sigma1_2_22 = 1.69;
%let mu2_1_1 = -1;
%let mu2_1_2 = -1;
%let mu2_2_1 = 3;
%let mu2_2_2 = 3;
%let sigma2_1_11 = 1;
%let sigma2_1_21 = 0;
%let sigma2_1_22 = 1;
%let sigma2_2_11 = 1.69;
%let sigma2_2_21 = 1;
%let sigma2_2_22 = 1.69;
%let seed = 1234;
%let T = 200000;
%let nSections = 200;

data gmDGP;
   retain cd1_1_11 cd1_1_21 cd1_1_22
          cd1_2_11 cd1_2_21 cd1_2_22
          cd2_1_11 cd2_1_21 cd2_1_22
          cd2_2_11 cd2_2_21 cd2_2_22;
   do t = 1 to &T.;
      sec = ceil(t*&nSections./&T.);
      if(t=1) then do;
         * initial probability distribution;
         p = (1-&a22.)/((1-&a11.)+(1-&a22.));;
         * Cholesky decomposition;
         cd1_1_11 = sqrt(&sigma1_1_11.);
         cd1_1_21 = &sigma1_1_21./cd1_1_11;
         cd1_1_22 = sqrt(&sigma1_1_22.-cd1_1_21*cd1_1_21);
         cd1_2_11 = sqrt(&sigma1_2_11.);
         cd1_2_21 = &sigma1_2_21./cd1_2_11;
         cd1_2_22 = sqrt(&sigma1_2_22.-cd1_2_21*cd1_2_21);
         cd2_1_11 = sqrt(&sigma2_1_11.);
         cd2_1_21 = &sigma2_1_21./cd2_1_11;
         cd2_1_22 = sqrt(&sigma2_1_22.-cd2_1_21*cd2_1_21);
         cd2_2_11 = sqrt(&sigma2_2_11.);
         cd2_2_21 = &sigma2_2_21./cd2_2_11;
         cd2_2_22 = sqrt(&sigma2_2_22.-cd2_2_21*cd2_2_21);
      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;
         * choose component;
         u = uniform(&seed.);
         if(u<=&c1_1.) then do;
            * (x,y) ~ N(mu, Sigma) at state 1, component 1;
            x = &mu1_1_1. + cd1_1_11*e1;
            y = &mu1_1_2. + cd1_1_21*e1+cd1_1_22*e2;
         end;
         else do;
            * (x,y) ~ N(mu, Sigma) at state 1, component 2;
            x = &mu1_2_1. + cd1_2_11*e1;
            y = &mu1_2_2. + cd1_2_21*e1+cd1_2_22*e2;
         end;
      end;
      else do;
         * choose component;
         u = uniform(&seed.);
         if(u<=&c2_1.) then do;
            * (x,y) ~ N(mu, Sigma) at state 2, component 1;
            x = &mu2_1_1. + cd2_1_11*e1;
            y = &mu2_1_2. + cd2_1_21*e1+cd2_1_22*e2;
         end;
         else do;
            * (x,y) ~ N(mu, Sigma) at state 2, component 2;
            x = &mu2_2_1. + cd2_2_11*e1;
            y = &mu2_2_2. + cd2_2_21*e1+cd2_2_22*e2;
         end;
      end;
      output;
      lags = s;
   end;
run;

data gmTrain;
   set gmDGP(where=(sec<=&nSections./2));
   keep sec t x y;
run;

data gmTest;
   set gmDGP(where=(sec>&nSections./2));
   keep sec t x y;
run;

The following statements estimate the GM HMM for the cross-sectional time series data after making the copies of the data tables in the library. The SECTION= option in the ID statement requests that the cross-sectional time series data be analyzed.

data mylib.gmTrain; set gmTrain; run;
data mylib.gmTest; set gmTest; run;

proc hmm data=mylib.gmTrain labelSwitch=(sort=none);
   id time=t section=sec;
   model x y / type=gaussianMixture nstate=2 ncomponent=2;
   optimize printLevel=3 printIterFreq=1;
   score outmodel=mylib.gmModel;
   filter out=mylib.gmFilter;
   evaluate out=mylib.gmEval;
run;

The model information is shown in Figure 40. There are 100 sections in the data.

The estimates of all parameters are shown in Figure 41; the estimates are very close to the parameter values in the DGP.

You can apply the estimates from the data table gmTrain to the testing data table gmTest to check the predictability in the following statements:

proc hmm data=mylib.gmTest;
   score inmodel=mylib.gmModel;
   decode out=mylib.gmDecode;
   filter out=mylib.gmFilter2;
   evaluate out=mylib.gmEval2;
run;

data gmHmmCheck;
   merge gmDGP(in=a where=(sec>&nSections./2)) mylib.gmDecode(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 sec t x y s correct ds;
run;

%let qFlipState = -1;
data gmHmmAccuracy;
   set gmHmmCheck;
   retain count 0 correctCount 0;
   count = count + 1;
   correctCount = correctCount + correct;
   if(count>&T./2-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 gmHmmCheck;
   set gmHmmCheck;
   if(&qFlipState.=1) then do;
      ds = -ds;
      correct = 1 - correct;
   end;
run;

proc print data = gmHmmAccuracy noobs; run;

The accuracy of the GM HMM for cross-sectional time series data, which is shown in Figure 42, is also around 98%.