HMM Procedure

This example revisits the model in Hamilton (1989) to illustrate how the regime-switching mean-adjusted autoregression model is used to estimate the business cycle with quarterly US real gross national product (GNP) data from 1951 to 1984.

The following statements import the data and calculate the GNP growth rate (in percentage) to which the model applies:

title 'Analysis of the Business Cycle';
data gnpHamilton;
   input date: anydtdte. gnp;
   format date yyq6.;
   dgnp = 100*(log(gnp)-log(lag(gnp)));
   if cmiss(of _all_) then delete;
   label date="Date" gnp="GNP" dgnp="Growth Rate of GNP";
datalines;
1951q1  1286.6
1951q2  1320.4
1951q3  1349.8

   ... more lines ...   

1982q4  3159.3
1983q1  3190.6
1983q2  3259.3
1983q3  3303.4
1983q4  3357.2
1984q1  3449.4
1984q2  3492.6
1984q3  3510.4
1984q4  3515.6
;

The following statements plot the data in Output 20.3.1:

proc sgplot data=gnpHamilton;
   series x=date y=dgnp;
   series x=date y=gnp / y2axis;
   refline 0 / axis=y lineattrs=(color=gray thickness=1);
run;

Hamilton’s model is a bistate regime-switching mean-adjusted AR(4) model which can be described as

where is the unobserved state.

In Hamilton (1989), the autoregressive parameters ’s and variance parameters ’s are assumed to be state-independent, which means and . The following statements re-create Hamilton’s result by setting the initial parameter values as the parameter estimates in Table 1 in Hamilton (1989) and also setting the MAXITER= option in the OPTIMIZE statement to zero. The FILTER statement conducts the probabilities of state on the basis of data available up to the current period. The DECODE statement outputs the best possible path of states on the basis of all observations. The SMOOTH statement finds the probabilities of states on the basis of all observations.

data mylib.gnpHamilton; set gnpHamilton; run;
proc hmm data=mylib.gnpHamilton labelswitch=(sort=desc(const));
   id time=date;
   model dgnp / armean=adjusted type=ar ylag=4 method=ml nstate=2
                stateIndependent=(ar cov);
   optimize maxiter=0;
   initial tpm={0.9049 0.0951, 0.2450 0.7550},
      ar={0.014 -0.058 -0.247 -0.213, 0.014 -0.058 -0.247 -0.213},
      const={1.1643, -0.3577},
      cov={0.5914, 0.5914};
   filter out=mylib.filter;
   decode out=mylib.decode;
   smooth out=mylib.smooth;
run;

The number of observations and model information are shown in Output 20.3.2.

The fit statistics are shown in Output 20.3.3.

The following statements plot the filtered probabilities of state 2 in Output 20.3.4, which is comparable to Figure 1 in Hamilton (1989). State 2 corresponds to recession because the intercept (and then the expected mean) in state 2 is negative.

data filter; set mylib.filter; run;
proc sort data=filter; by date; run;
proc sgplot data=filter;
   series x=date y=state2;
   yaxis label="Filtered Prob. of Recession";
run;

In order to show how states correspond to the signs of GNP growth rate, the following statements overlay the decoded states on the GNP growth rate. Recession is from peak to trough, and boom is from trough to peak. As shown in Output 20.3.5, the recession includes most periods that have a negative GNP growth rate.

data plot2state;
   merge mylib.decode mylib.gnpHamilton;
   if state=2 then do; recession=1; boom=0; end;
   else if state=1 then do; recession=0; boom=1; end;
   rec = recession*(-3);
   boo = boom*3;
   by date;
   if cmiss(of _all_) then delete;
run;

proc sgplot data=plot2state;
   band x=date upper=0 lower=rec / legendlabel='Recession' type=step
                                   fillattrs=(color=red);
   band x=date upper=boo lower=0 / legendlabel='Boom' type=step
                                   fillattrs=(color=green);
   series x=date y=dgnp / lineattrs=(thickness=2);
run;

Hamilton (1989) also compares his choice of business-cycle dates with the historical business-cycle dates that are reported from the National Bureau of Economic Research (NBER). The following data step inputs the business-cycle dates from NBER:

data mylib.nber;
   input date: anydtdte. recessionNBER;
   format date yyq6.;
   label date="Date" recessionNBER="Recession Reported from NBER";
datalines;
1951q1  0
1951q2  0
1951q3  0
1951q4  0
1952q1  0
1952q2  0

   ... more lines ...   

1982q4  1
1983q1  0
1983q2  0
1983q3  0
1983q4  0
1984q1  0
1984q2  0
1984q3  0
1984q4  0
;

The following statements plot together the recession dates that are found by the model and those that are reported from NBER. As shown in Output 20.3.6, the recession dates that are found by Hamilton’s model (the upper part) is extremely close to the dates that are reported from NBER (the lower part).

data gnpsmooth;
   merge mylib.smooth mylib.nber mylib.gnpHamilton;
   if state2>=0.5 then recession=1; else recession=0;
   by date;
   gnp2=gnp;
   if cmiss(of _all_) then delete;
   label recession="Recession Found by the Model" gnp2="GNP";
run;
proc transpose data=gnpsmooth  out=gnpsmooth_t;
   by date;
   var recessionNBER recession;
run;
proc transpose data=gnpsmooth  out=gnpsmooth_t2;
   by date;
   var gnp gnp2;
run;
data gnpplot;
   merge gnpsmooth_t(rename=(_name_=name1 _label_=label1 col1=recession))
         gnpsmooth_t2(rename=(_name_=name2 _label_=label2 col1=gnp));
   by date;
   rec=recession*3520;
run;
proc sgpanel data=gnpplot;
   panelby label1 / layout=rowlattice uniscale=column novarname;
   band x=date upper=rec lower=0 / legendlabel='Recession' type=step;
   series x=date y=gnp / legendlabel='GNP';
   rowaxis display=(nolabel) values=(1000 to 3520 by 500);
run;

The following statements use PROC HMM to estimate Hamilton’s model while using the STATEINDEPENDENT=(AR COV) option in the MODEL statement so that the autoregression parameter and the variance are fixed across the states:

proc hmm data=mylib.gnpHamilton labelswitch=(sort=desc(const));
   id time=date;
   model dgnp / armean=adjusted type=ar ylag=4 method=ml nstate=2
             stateIndependent=(ar cov);
   optimize  printIterFreq=1 printlevel=3 ;
   filter out=mylib.filter;
   decode out=mylib.decode;
   smooth out=mylib.smooth;
run;

The parameter estimates are shown in Output 20.3.7, and the fit statistics are shown in Output 20.3.8. The log likelihood value and the information criterion in Output 20.3.8 show that this solution is better than the one in Output 20.3.3 when MAXITER=0.

The following statements plot the filtered probabilities of state 2 (recession). Output 20.3.9 shows the results.

data filter; set mylib.filter; run;
proc sort data=filter; by date; run;
proc sgplot data=filter;
   series x=date y=state2;
   yaxis label="Filtered Prob. of Recession";
run;

The following statements overlay the decoded states on the GNP growth rate. Recession is from peak to trough, and boom is from trough to peak.

data plot2state;
   merge mylib.decode mylib.gnpHamilton;
   if state=2 then do; recession=1; boom=0; end;
   else if state=1 then do; recession=0; boom=1; end;
   rec = recession*(-3);
   boo = boom*3;
   by date;
   if cmiss(of _all_) then delete;
run;

proc sgplot data=plot2state;
   band x=date upper=0 lower=rec / legendlabel='Recession' type=step
                                   fillattrs=(color=red);
   band x=date upper=boo lower=0 / legendlabel='Boom' type=step
                                   fillattrs=(color=green);
   series x=date y=dgnp / lineattrs=(thickness=2);
run;

Output 20.3.10 shows that more states switch than are shown in Hamilton’s decoded plot. The reason is that the estimates indicate less persistence of each state.

The following statements plot the recession dates that are found by the model and those that are reported from NBER:

data gnpsmooth;
   merge mylib.smooth mylib.nber mylib.gnpHamilton;
   if state2>=0.5 then recession=1; else recession=0;
   by date;
   gnp2=gnp;
   if cmiss(of _all_) then delete;
   label recession="Recession Found by the Model" gnp2="GNP";
run;
proc transpose data=gnpsmooth  out=gnpsmooth_t;
   by date;
   var recessionNBER recession;
run;
proc transpose data=gnpsmooth  out=gnpsmooth_t2;
   by date;
   var gnp gnp2;
run;
data gnpplot;
   merge gnpsmooth_t(rename=(_name_=name1 _label_=label1 col1=recession))
         gnpsmooth_t2(rename=(_name_=name2 _label_=label2 col1=gnp));
   by date;
   rec=recession*3520;
run;
proc sgpanel data=gnpplot;
   panelby label1 / layout=rowlattice uniscale=column novarname;
   band x=date upper=rec lower=0 / legendlabel='Recession' type=step;
   series x=date y=gnp / legendlabel='GNP';
   rowaxis display=(nolabel) values=(1000 to 3520 by 500);
run;

As shown in Output 20.3.11, the recession dates that are found by the model (the upper part) are extremely close to the dates that are reported from NBER (the lower part).