Replicare l'esempio di Markov Switching Model di Hamilton usando il pacchetto MSwM in R

Question

Sto provando a stimare il Markov Switching Model di base di Hamilton (1989) come post in E-views webpage. Questo modello è di per sé una replica esatta dell'esistente in RATS.

Questa è la serie storica dell'esempio:

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

voglio usare il pacchetto MSwM, così ho scritto il seguente codice:

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) 

ottengo un risultato che è molto diverso da ottenere in Eviews o 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 

La differenza principale si ottiene nell'intercetta, perché in entrambi i casi gime si ottiene un valore positivo anziché valori in Eviews o RATS. Questa differenza è dovuta alla massimizzazione utilizzata da algortihm (EM in MsWm)? o ho fatto qualche errore nel mio codice R?

Grazie mille.

Answer 1

La differenza che vedo è che il modello che si sta definendo contiene un'intercetta di commutazione, mentre il modello di Hamilton (1989) specifica invece una media di commutazione. Cioè, il modello è:

e (1989) il modello di Hamilton è definito come:

In un modello AR parametri alpha e mu prenderanno, in generale, valori diversi. Questo può essere un po 'di confusione in R come discusso here.

Prendendo aspettative nel modello (e omettendo per semplicità il termine di commutazione S_t) arriviamo alla seguente relazione:

Su questo rapporto, ci si poteva aspettare di essere in grado di recuperare il significare. Tuttavia, in questo caso la commutazione intercetta non porta ai mezzi di commutazione trovati 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 

Questa mappatura solito può essere applicato, ad esempio, con un AR (4) modello:

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 

Sembra che in questo caso il modello deve essere definito in termini di media in ordine per ottenere stime del parametro di commutazione vicine a quelle riportate nel documento di riferimento.

Se la funzione msmFit consentito in ingresso il risultato restituito da arima, potrebbe essere utilizzato come segue:

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

non so un modo semplice per definire un modello AR con media utilizzando lm, che è l'output richiesto per utilizzare msmFit.

Penso che questa differenza nella parametrizzazione del modello sia più probabile che spieghi la differenza nei risultati piuttosto che l'utilizzo dell'algoritmo EM.