Replicating the example of Markov Switching Model of Hamilton using MSwM package in R

Question

I'm trying to estimate the basic Markov Switching Model of Hamilton (1989) as is post in E-views webpage. This model is itself is an exact replication of the existing in RATS.

This is the time series of the example:

gnp <- 
structure(c(2.59316410021381, 2.20217123302681, 0.458275619103479, 
0.968743815568942, -0.241307564718414, 0.896474791426144, 2.05393216767198, 
1.73353647046698, 0.938712869506845, -0.464778333117193, -0.809834082445603, 
-1.39763692441103, -0.398860927649558, 1.1918415768741, 1.4562004729396, 
2.1180822079447, 1.08957867423914, 1.32390272784813, 0.87296368144358, 
-0.197732729861307, 0.45420214345009, 0.0722187603196887, 1.10303634435563, 
0.820974907499614, -0.0579579499110212, 0.584477722838197, -1.56192668045796, 
-2.05041027007508, 0.536371845140342, 2.3367684244086, 2.34014568267516, 
1.23392627573662, 1.88696478737248, -0.459207909351867, 0.84940472194713, 
1.70139850766727, -0.287563102546191, 0.095946277449187, -0.860802907461483, 
1.03447124467041, 1.23685943797014, 1.42004498680119, 2.22410642769683, 
1.3021017302965, 1.0351769691057, 0.925342521818, -0.165599507925585, 
1.3444381723048, 1.37500136316918, 1.73222186043569, 0.716056342342333, 
2.21032138350616, 0.853330335823775, 1.00238777849592, 0.427254413549543, 
2.14368353713136, 1.4378918561536, 1.5795993028646, 2.27469837381376, 
1.95962653201067, 0.2599239932111, 1.01946919515563, 0.490163994319276, 
0.563633789161385, 0.595954621290765, 1.43082852218349, 0.562301244017229, 
1.15388388887095, 1.68722847001462, 0.774382052478202, -0.0964704476805431, 
1.39600141863966, 0.136467982223878, 0.552237133917267, -0.399448716111952, 
-0.61671104590512, -0.0872256083215416, 1.21018349098461, -0.907297546921259, 
2.64916154469762, -0.00806939681695959, 0.511118931407946, -0.00401437145032572, 
2.1682142321342, 1.92586729194597, 1.03504719187207, 1.85897218652101, 
2.32004929969819, 0.255707901889092, -0.0985527428151145, 0.890736834018326, 
-0.55896483237131, 0.283502534230679, -1.31155410054958, -0.882787789285689, 
-1.97454945511993, 1.01275266533046, 1.68264718400186, 1.38271278970291, 
1.86073641586006, 0.444737715592073, 0.414490009766608, 0.992022769383933, 
1.36283572253682, 1.59970527327726, 1.98845814838348, -0.256842316681229, 
0.877869502339381, 3.10956544706826, 0.853244770655281, 1.23337321374495, 
0.0031430232743432, -0.0943336967005583, 0.898833191548979, -0.190366278407953, 
0.997723787687709, -2.39120056095144, 0.0664967330277127, 1.26136016443398, 
1.91637832265846, -0.334802886728505, 0.44207108280265, -1.40664914211265, 
-1.52129894225829, 0.299198686266393, -0.801974492802505, 0.152047924379708, 
0.985850281223592, 2.1303461510993, 1.34397927090998, 1.61550521216825, 
2.70930096486278, 1.24461416484445, 0.508354657516633, 0.148021660957899
), .Tsp = c(1951.25, 1984.75, 4), class = "ts")

I want to use the MSwM package, so I wrote the following code:

library(MSwM) #Load the package  
# Create the model with only an intercept (that after will be switching) 
mod=lm(gnp~1)
# Estimate the Markov Switching Model with only an intercept switching, 
# four lags and two regimes as in Hamilton. 
mod.mswm=msmFit(mod,k=2,p=4,sw=c(T,F,F,F,F,F), control=list(parallel=F)) 
summary(mod.mswm)

I get a result that is very different to obtained in Eviews or RATS:

Coefficients:
Regime 1 
---------
               Estimate Std. Error t value  Pr(>|t|)    
(Intercept)(S)   0.5747     1.0044  0.5722 0.5671865    
gnp_1            0.3097     0.0903  3.4297 0.0006042 ***
gnp_2            0.1273     0.0900  1.4144 0.1572445    
gnp_3           -0.1213     0.0867 -1.3991 0.1617830    
gnp_4           -0.0892     1.6918 -0.0527 0.9579709    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.98316
Multiple R-squared: 0.1437

Standardized Residuals:
        Min          Q1         Med          Q3         Max 
-1.86974671 -0.37107376  0.03466299  0.39090950  1.67876663 

Regime 2 
---------
               Estimate Std. Error t value  Pr(>|t|)    
(Intercept)(S)   0.5461     1.0044  0.5437 0.5866479    
gnp_1            0.3097     0.0903  3.4297 0.0006042 ***
gnp_2            0.1273     0.0900  1.4144 0.1572445    
gnp_3           -0.1213     0.0867 -1.3991 0.1617830    
gnp_4           -0.0892     1.6918 -0.0527 0.9579709    
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.98316
Multiple R-squared: 0.1431

Standardized Residuals:
        Min          Q1         Med          Q3         Max 
-2.51219057 -0.46185366  0.06749067  0.52368275  2.11071358 

Transition probabilities:
          Regime 1  Regime 2
Regime 1 0.3879799 0.3651762
Regime 2 0.6120201 0.6348238

The main difference is obtained in the intercept, because in both regimes a positive value is obtained instead of values in Eviews or RATS. This difference is due to maximization algortihm used (EM in MsWm)? or I have done some mistake in my R-Code?

Thanks a lot.

Answer 1

The difference that I see is that the model that you are defining contains a switching intercept, while the model of Hamilton (1989) specifies a switching mean instead. That is, your model is:

and Hamilton's (1989) model is defined as:

In an AR model the parameters alpha and mu will take, in general, different values. This may be somewhat confusing in R as discussed here.

By taking expectations in your model (and omitting for simplicity the switching term S_t) we arrive to the following relationship:

Upon this relationship, we could expect to be able to recover the mean. However, in this case the switching intercepts does not lead to the switching means found in Hamilton (1989).

0.5747 / (1 - sum(c(0.3097, 0.1273, -0.1213, -0.0892)))
#[1] 0.7429864
0.5461 / (1 - sum(c(0.3097, 0.1273, -0.1213, -0.0892)))
#[1] 0.7060116

This mapping can usually be applied, for example, with an AR(4) model:

fit <- lm(gnp[5:135] ~ 1 + gnp[4:134] + gnp[3:133] + gnp[2:132] + gnp[1:131])
fit
# Coefficients:
# (Intercept)   gnp[4:134]   gnp[3:133]   gnp[2:132]   gnp[1:131]  
#     0.55679      0.30974      0.12726     -0.12126     -0.08923 
#
# the mapping from the intercept to mean leads to a value close to the sample mean
coef(fit)[1]/(1 - sum(coef(fit)[-1]))
# 0.7198458 
mean(gnp)
# 0.7445979
# or close to the mean in an AR(4) model, (labelled as intercept)
arima(gnp, order = c(4,0,0), include.mean = TRUE)
# Coefficients:
#          ar1     ar2      ar3      ar4  intercept
#       0.3188  0.1226  -0.1191  -0.0895     0.7441
# s.e.  0.0860  0.0900   0.0898   0.0872     0.1108

It seems that in this case the model should be defined in terms of the mean in order to get estimates of the switching parameter close to those reported in the reference paper.

If the function msmFit allowed as input the result returned by arima, it could be used as follows:

fit <- arima(gnp, order = c(4,0,0), include.mean = TRUE)
msmFit(fit, k = 2, p = 0, sw = c(T,F,F,F,F,F))

I don't know a straightforward way to define an AR model with mean using lm, which is the output required to use msmFit.

I think that this difference in the parameterization of the model is more likely to explain the difference in the results rather than the use of the EM algorithm.