HMM Procedure

Let , , denote a random discrete variable whose sample space is a finite set of nominal classes. These values from this sample space can be represented by the finite integers . can be modeled in a categorical distribution,

where the category probability is the probability of observing , conditional on the latent variable . The category probability matrix (CPM) is used to define the state-dependent categorical distribution of observable variables. The variable is the so-called state and 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 is from a finite support, the model is called the finite hidden Markov model (finite HMM). The finite HMM is sometimes described as .

Consider a standard die with six faces, numbered 1 through 6. A fair die has an equal probability of displaying 1 to 6; that is, . In order to tell whether a die is fair, checking whether the opportunity of displaying each number is equal might not be enough.

The following data record the results from the roll of a die 1,000 times:

data mylib.finite1;
   input t  y ;
datalines;
1   6
2   3
3   6
4   5
5   1
6   2
7   1
8   6

   ... more lines ...   

993 2
994 3
995 6
996 6
997 5
998 5
999 2
1000    6
;

The following statements estimate the unistate finite HMM:

proc hmm data=mylib.finite1;
   id time=t;
   model y / type=finite nstate=1 nCategory=6;
   estimate outall=mylib.myEstimate;
run;

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

The estimates of the initial state probability vector (ISPV) and the transition probability matrix are trivial. The category probability matrix is shown in Figure 52, which shows the probabilities of occurrence of each number.

The estimates of all parameters are shown in Figure 53, which displays columns for parameter name, estimated value, standard error, t value, and p-value.

The following statements run a Wald test under the null hypothesis that all category probabilities are equal to 1/6:

proc iml;
  use mylib.myEstimate;
  read all var {Estimate cov_2 cov_3 cov_4 cov_5 cov_6} into m;
  close mylib.myEstimate;

cov = m[2:6,2:6];
invcov = inv(cov);
v = m[2:6,1];
c = 1/6 // 1/6 // 1/6 // 1/6 // 1/6 ;
v1 = v-c;

stat = v1`*invcov *v1;
df = 5;
pvalue = 1- cdf('CHISQ',stat,df);
varname={"stat" "df" "pvalue"};
create table var varname;
append;
close table;
quit;

proc print data=table noobs label;
   label stat="Stat" df = "DF" pvalue = "Pr > ChiSq";
run;

The p-value shown in Figure 54 is very large, which means that the test cannot reject the null hypothesis that the category probabilities are equal to 1/6 at the significance level. The die looks fair.

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

The following code estimates a bistate finite HMM:

proc hmm data=mylib.finite1;
   id time=t;
   model y / type=finite nstate=2 nCategory=6;
   score outmodel=mylib.myScore;
run;

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

The estimates of the ISPV are shown in Figure 57.

The estimates of the TPM are shown in Figure 58.

The estimates of the CPM are shown in Figure 59.

The estimates of all parameters are shown in Figure 60, which displays columns for parameter name, estimated value, standard error, t value, and p-value.

In Figure 60, you see that the category probabilities for small numbers (1, 2, 3) in the first state and the category probabilities for large numbers (4, 5, 6) in the second state are not significant at the significance level. This indicates that the small numbers appear only in one state and the large numbers appear only in the other state.

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

All the information criteria from the bistate finite HMM are smaller than those from the unistate finite HMM. This indicates that the bistate finite HMM is a better model. The tristate finite HMM, which is not shown here, is also beaten by the bistate finite HMM for this case.

The casino game sic bo involves whether you will get a small or big number when you roll a die. A small bet is when you predict that the die will display a 1, 2, or 3. A big bet is when you predict that the die will display a 4, 5, or 6. If your guess is correct, you get a reward; otherwise you lose the bet.

You can use PROC HMM to test the bistate finite HMM’s ability to forecast the result of this game. The following data record the results from the roll of a die 100 times:

data mylib.finite2;
   input t  y ;
datalines;
1   6
2   4
3   6
4   1
5   1
6   1
7   1
8   4

   ... more lines ...   

93  6
94  6
95  5
96  2
97  1
98  6
99  4
100 6
;

The following statements perform one-step online (continuous) forecasting of betting big with the estimated parameters from the preceding PROC HMM statements and compare the prediction with the finite2 data set. If the forecasting is correct, you gain one unit. Otherwise you lose one unit.

proc hmm data=mylib.finite2;
   id time=t;
   score inmodel=mylib.myScore;
   forecast out=mylib.myforecast online;
run;
data mylib.forecast;
   set mylib.myforecast;
   if (Category1+Category2+Category3>0.5) then bet = 0;
   else bet = 1;
run;
data mylib.forecast;
   set mylib.forecast;
   t=t+step; *so t equals the time when forecast is on;
 run;
data mylib.finite2;
   set mylib.finite2;
 if(y<=3) then big=0;
 else big=1;
run;
data fnHmmCheck;
   merge mylib.finite2(in=a) mylib.forecast(in=b);
   by t;
run;
data fnHmmCheck;
   set fnHmmCheck;
   *for first round bet according to inital probability;
   if(t=1) then bet=1;
   if(big=bet) then gain=1;
   else gain=-1;
   keep t y gain;
run;
data totalgain;
   set fnHmmCheck;
   Totalgain + gain;
   if (t=100) then do
      Zscore = Totalgain/10;
      pvalue = 2*(1- cdf('NORMAL',Zscore));
      output totalgain;
   end;
   keep Totalgain Zscore pvalue;
run;
proc print data = totalgain noobs label;
   label totalgain = "Total gain" Zscore="Z-score" pvalue="Pr > |z|";
run;

The total gain from rolling the die 100 times is shown in Figure 62. Under the null hypothesis that the die is fair, according to the central limit theorem the expected return of betting 100 times is 0 with a standard error of 10 units. The values shown in Figure 62 indicate that the null hypothesis is rejected even at the significance level. This means that the finite HMM successfully detects the loaded die.

The data generating process behind the observed data sets is as follows:

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

The following statements simulate the original time series data for this example:

%let pi1 = 0.5;*initial probabilty in state 1;
%let a11 = 0.75;
%let a22 = 0.75;

%macro simData(name, T, seed);
data mylib.&name.(keep = t y);
  p = &pi1.;
  call streaminit(&seed);
  do t = 1 to &T.;
   if(t=1) then
      /* initial probability distribution */
     p = &pi1.;
   else
      /* transition probability matrix */
     if(lagx=1) then p = &a11.;
     else p = 1-&a22.;
     u = rand("Uniform");
     if(u<=p) then x=1;
     *Using the random value u to control regimes, set x=number of regimes;
     else x = 2;
     if(x=1) then
     do;
      /* at state 1 */
     y = rand("Table", 1/3, 1/3, 1/3, 0 , 0 , 0);
   end;
   else do;
      /* at state 2 */
     y = rand("Table", 0 , 0 , 0 , 1/3, 1/3, 1/3);
   end;
   output;
   lagx=x;
end;
run;
%mend;

%simData(finite1,1000, 3);
%simData(finite2,100, 4);

PROC HMM assumes that values of the response variables are consecutive natural numbers that start from 1. If this is not the case, you might need to preprocess the input data set before running PROC HMM. You can use the SAS macro %FiniteHMM for automatic preprocessing and postprocessing. For more information about %FiniteHMM, see the section The %FiniteHMM Macro.