R語言如何做馬爾科夫轉換模型markov switching model


 原文鏈接:http://tecdat.cn/?p=6962

 

假設 有時間序列數據,如下所示。經驗表明,目標變量y似乎與解釋變量x有關。然而,乍一看,y的水平在中間移動,所以它似乎並不總是有固定的關系(背后有多個狀態)。


 


上面的樣本數據創建如下。數據根據時間改變x和y之間的關系。

 

x <- rpois(500, lambda = 10)  
y1 <- x * 4 + 20     
y2 <- x * 2 + 60    

 
noise <- rnorm(1:500, mean = 10, sd = 5)
y1 <- y1 + noise
y2 <- y2 + noise

 y <- c(y1[1:200], y2[201:400], y1[401:500])
 observed <- data.frame(x = x, y = y)

x和y1,y2之間的關系如下圖所示。如果您知道x和y有兩種狀態,則x和y看起來像這樣。

數據 

​ 在馬爾可夫轉換模型中,觀察數據被認為是從幾個狀態生成的,並且如上所示很好地分離。

觀察到的數據


創建馬爾可夫轉換模型

 

模型公式 

 

 
# Call:
# lm(formula = y ~ x, data = observed)
# 
# Residuals:
#     Min      1Q  Median      3Q     Max 
# -24.303  -9.354  -1.914   9.617  29.224 
# 
# Coefficients:
#             Estimate Std. Error t value Pr(>|t|)    
# (Intercept)  45.7468     1.7202   26.59   <2e-16 ***
# x             3.2262     0.1636   19.71   <2e-16 ***
# ---
# Signif. codes:  
# 0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 11.51 on 498 degrees of freedom
# Multiple R-squared:  0.4383, Adjusted R-squared:  0.4372 
# F-statistic: 388.7 on 1 and 498 DF,  p-value: < 2.2e-16

 

 參數的含義是

  • k:馬爾可夫轉換模型的狀態數。在這里,它被指定為后面有兩個狀態。
  • sw:使用邏輯指定每個參數在狀態更改時是否更改
  • p:AR模型系數
  • family:(在GLM的情況下)概率分布族
 # Markov Switching Model
 
# 
#        AIC      BIC    logLik
#   3038.846 3101.397 -1513.423
# 
# Coefficients:
# 
# Regime 1 
# ---------
#                Estimate Std. Error t value Pr(>|t|)    
# (Intercept)(S)  69.3263     4.0606 17.0729   <2e-16 ***
# x(S)             2.1795     0.1187 18.3614   <2e-16 ***
# y_1(S)          -0.0103     0.0429 -0.2401   0.8103    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 4.99756
# Multiple R-squared: 0.6288
# 
# Standardized Residuals:
#           Min            Q1           Med            Q3           Max 
# -1.431396e+01 -2.056292e-02 -1.536781e-03 -1.098923e-05  1.584478e+01 
# 
# Regime 2 
# ---------
#                Estimate Std. Error t value Pr(>|t|)    
# (Intercept)(S)  30.2820     1.7687 17.1210   <2e-16 ***
# x(S)             3.9964     0.0913 43.7722   <2e-16 ***
# y_1(S)          -0.0045     0.0203 -0.2217   0.8245    
# ---
# Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
# 
# Residual standard error: 4.836684
# Multiple R-squared: 0.8663
# 
# Standardized Residuals:
#           Min            Q1           Med            Q3           Max 
# -13.202056966  -0.771854514   0.002211602   1.162769110  12.417873232 
# 
# Transition probabilities:
#             Regime 1    Regime 2
# Regime 1 0.994973376 0.003347279
# Regime 2 0.005026624 0.996652721

輸出中的制度1和制度2表示后面的兩個狀態 。 

# Regime 1 
# ---------
#                Estimate Std. Error t value Pr(>|t|)    
# (Intercept)(S)  69.3263     4.0606 17.0729   <2e-16 ***
# x(S)             2.1795     0.1187 18.3614   <2e-16 ***
# y_1(S)          -0.0103     0.0429 -0.2401   0.8103   

y1 <- x * 4 + 20 可以看到Regime 2 與之兼容。

 可以說從調整后的R平方值整體上有所改善。

# Regime 2 
# ---------
#                Estimate Std. Error t value Pr(>|t|)    
# (Intercept)(S)  30.2820     1.7687 17.1210   <2e-16 ***
# x(S)             3.9964     0.0913 43.7722   <2e-16 ***
# y_1(S)          -0.0045     0.0203 -0.2217   0.8245    

 模型

對於每個regime,目標變量+指定的解釋變量和處於該狀態的概率以陰影繪制


每個時間點的概率 


每次獲取狀態和更改點

如果你想知道你在某個特定時間點所在的regime,那么就選擇那個時刻概率最高的 。

> probable
  [1] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
 [30] 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2 2
...

異常值/變化點是Regime更改的時間 

c(FALSE, diff(probable) != 0)
  [1] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
 [11] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
...
[181] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[191] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE
[201] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
...
[381] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
[391] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE  TRUE
[401] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE
...
[491] FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE FALSE

因此,我們可以看到檢測到在第一次數據創建時指定的變化點(201,401th)附近的點。

如果您有任何疑問,請在下面發表評論。

 


免責聲明!

本站轉載的文章為個人學習借鑒使用,本站對版權不負任何法律責任。如果侵犯了您的隱私權益,請聯系本站郵箱yoyou2525@163.com刪除。



 
粵ICP備18138465號   © 2018-2025 CODEPRJ.COM