在 R 中使用 MSwM 包复制 Hamilton 的 Markov Switching Model 示例

Question

我正在尝试估计 Hamilton (1989) 的基本马尔可夫切换模型，就像 E-views webpage 中的帖子一样。该模型本身就是 RATS 中现有模型的精确复制。

这是示例的时间序列:

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")

我想使用MSwM包，所以写了如下代码:

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)

我得到的结果与在 Eviews 或 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

主要区别在于截距，因为在两种情况下都获得了正值，而不是 Eviews 或 RATS 中的值。
这种差异是由于使用了最大化算法(MsWm 中的 EM)？或者我在我的 R 代码中犯了一些错误？

非常感谢。

Answer 1

我看到的不同之处在于您定义的模型包含一个切换截距，而 Hamilton (1989) 的模型则指定了一个切换均值。也就是说，您的模型是:



和 Hamilton (1989) 模型定义为:



在 AR 模型中，参数 alpha和 mu通常，将采用不同的值。
正如所讨论的那样，这在 R 中可能有些困惑 here .

通过在您的模型中获取期望(并为简单起见省略切换项 S_t )
我们得出以下关系:



根据这种关系，我们可以期望能够恢复均值。然而，在这种情况下，切换截距不会导致 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

例如，这种映射通常可以应用于 AR(4) 模型:

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

在这种情况下，似乎应该根据平均值定义模型，以便获得与引用文件中报告的接近的开关参数估计值。

如果函数 msmFit允许作为输入 arima 返回的结果，它可以如下使用:

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))

我不知道使用 lm 定义具有均值的 AR 模型的直接方法，这是使用 msmFit 所需的输出.

我认为模型参数化的这种差异更可能解释结果的差异，而不是使用 EM 算法。