zoukankan      html  css  js  c++  java
  • 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)附近的点。

    如果您有任何疑问,请在下面发表评论。

  • 相关阅读:
    Ansible跳板机自动部署
    nginx展示文件目录
    【转】消息钩子注册浅析
    windows临界区
    windbg定位死锁
    Windows工作集内存
    我的spring boot,杨帆、起航!
    CursorFileManager对cursor文件的读写
    eclipse执行maven install命令时跳过test
    bASE--Risk
  • 原文地址:https://www.cnblogs.com/tecdat/p/11543515.html
Copyright © 2011-2022 走看看