R语言如何做马尔可夫转换模型markov switching model|附代码数据

模式流星
• 阅读 843

全文链接:http://tecdat.cn/?p=6962

最近我们被客户要求撰写关于马尔可夫转换模型的研究报告,包括一些图形和统计输出。
假设 有时间序列数据,如下所示。经验表明,目标变量y似乎与解释变量x有关。然而,乍一看,y在水平中间波动,所以它似乎并不总是有稳定的关系(背后有多个状态)

R语言如何做马尔可夫转换模型markov switching model|附代码数据

上面的样本数据创建如下。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之间的关系如下图所示。

数据

R语言如何做马尔可夫转换模型markov switching model|附代码数据 

在马尔可夫转换模型中,观察数据被认为是从几个状态生成的,并且如上所示可以很好地分离。

观察到的数据

R语言如何做马尔可夫转换模型markov switching model|附代码数据


点击标题查阅往期内容

R语言如何做马尔可夫转换模型markov switching model|附代码数据

【视频】马尔可夫链蒙特卡罗方法MCMC原理与R语言实现|数据分享

R语言如何做马尔可夫转换模型markov switching model|附代码数据

左右滑动查看更多

R语言如何做马尔可夫转换模型markov switching model|附代码数据

01

R语言如何做马尔可夫转换模型markov switching model|附代码数据

02

R语言如何做马尔可夫转换模型markov switching model|附代码数据

03

R语言如何做马尔可夫转换模型markov switching model|附代码数据

04

R语言如何做马尔可夫转换模型markov switching model|附代码数据

创建马尔可夫转换模型

 

模型公式 

 

# 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的情况下)概率分布族

<!---->

 # 马尔可夫转换模型
# 
#        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   

可以看到区制2 与y1 <- x * 4 + 20匹配。

从调整后的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    

 模型

对于每个状态,处于该状态的概率以阴影绘制

R语言如何做马尔可夫转换模型markov switching model|附代码数据

每个时间点的概率 

R语言如何做马尔可夫转换模型markov switching model|附代码数据

每次获取状态和更改点

如果你想知道你在某个特定时间点所在的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
...

异常值/变化点是状态更改的时间 

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

因此,我们可以看到检测到在第一次数据创建时指定的变化点。


R语言如何做马尔可夫转换模型markov switching model|附代码数据

点击文末 “阅读原文”

获取全文完整代码数据资料。

本文选自《R语言如何做马尔可夫转换模型markov switching model》。

点击标题查阅往期内容

matlab用马尔可夫链蒙特卡罗 (MCMC) 的Logistic逻辑回归模型分析汽车实验数据\
【视频】马尔可夫链蒙特卡罗方法MCMC原理与R语言实现|数据分享\
R语言BUGS/JAGS贝叶斯分析: 马尔科夫链蒙特卡洛方法(MCMC)采样\
马尔可夫Markov区制转移模型分析基金利率\
马尔可夫区制转移模型Markov regime switching\
时变马尔可夫区制转换MRS自回归模型分析经济时间序列\
马尔可夫转换模型研究交通伤亡人数事故时间序列预测\
如何实现马尔可夫链蒙特卡罗MCMC模型、Metropolis算法?\
Matlab用BUGS马尔可夫区制转换Markov switching随机波动率模型、序列蒙特卡罗SMC、M H采样分析时间序列\
R语言BUGS序列蒙特卡罗SMC、马尔可夫转换随机波动率SV模型、粒子滤波、Metropolis Hasting采样时间序列分析\
matlab用马尔可夫链蒙特卡罗 (MCMC) 的Logistic逻辑回归模型分析汽车实验数据\
stata马尔可夫Markov区制转移模型分析基金利率\
PYTHON用时变马尔可夫区制转换(MRS)自回归模型分析经济时间序列\
R语言使用马尔可夫链对营销中的渠道归因建模\
matlab实现MCMC的马尔可夫转换ARMA - GARCH模型估计\
R语言隐马尔可夫模型HMM识别不断变化的股票市场条件\
R语言中的隐马尔可夫HMM模型实例\
用机器学习识别不断变化的股市状况—隐马尔科夫模型(HMM)\
Matlab马尔可夫链蒙特卡罗法(MCMC)估计随机波动率(SV,Stochastic Volatility) 模型\
MATLAB中的马尔可夫区制转移(Markov regime switching)模型\
Matlab马尔可夫区制转换动态回归模型估计GDP增长率\
R语言马尔可夫区制转移模型Markov regime switching\
stata马尔可夫Markov区制转移模型分析基金利率\
R语言如何做马尔可夫转换模型markov switching model\
R语言隐马尔可夫模型HMM识别股市变化分析报告\
R语言中实现马尔可夫链蒙特卡罗MCMC模型

点赞
收藏
评论区
推荐文章
blmius blmius
3年前
MySQL:[Err] 1292 - Incorrect datetime value: ‘0000-00-00 00:00:00‘ for column ‘CREATE_TIME‘ at row 1
文章目录问题用navicat导入数据时,报错:原因这是因为当前的MySQL不支持datetime为0的情况。解决修改sql\mode:sql\mode:SQLMode定义了MySQL应支持的SQL语法、数据校验等,这样可以更容易地在不同的环境中使用MySQL。全局s
Wesley13 Wesley13
3年前
MySQL部分从库上面因为大量的临时表tmp_table造成慢查询
背景描述Time:20190124T00:08:14.70572408:00User@Host:@Id:Schema:sentrymetaLast_errno:0Killed:0Query_time:0.315758Lock_
美凌格栋栋酱 美凌格栋栋酱
6个月前
Oracle 分组与拼接字符串同时使用
SELECTT.,ROWNUMIDFROM(SELECTT.EMPLID,T.NAME,T.BU,T.REALDEPART,T.FORMATDATE,SUM(T.S0)S0,MAX(UPDATETIME)CREATETIME,LISTAGG(TOCHAR(
皕杰报表之UUID
​在我们用皕杰报表工具设计填报报表时,如何在新增行里自动增加id呢?能新增整数排序id吗?目前可以在新增行里自动增加id,但只能用uuid函数增加UUID编码,不能新增整数排序id。uuid函数说明:获取一个UUID,可以在填报表中用来创建数据ID语法:uuid()或uuid(sep)参数说明:sep布尔值,生成的uuid中是否包含分隔符'',缺省为
深度学习技术开发与应用
关键点1.强化学习的发展历程2.马尔可夫决策过程3.动态规划4.无模型预测学习5.无模型控制学习6.价值函数逼近7.策略梯度方法8.深度强化学习DQN算法系列9.深度策略梯度DDPG,PPO等第一天9:0012:0014:0017:00一、强化学习概述1.强化学习介绍2.强化学习与其它机器学习的不同3.强化学习发展历史4.强化学习典
Jacquelyn38 Jacquelyn38
4年前
2020年前端实用代码段,为你的工作保驾护航
有空的时候,自己总结了几个代码段,在开发中也经常使用,谢谢。1、使用解构获取json数据let jsonData  id: 1,status: "OK",data: 'a', 'b';let  id, status, data: number   jsonData;console.log(id, status, number )
Stella981 Stella981
3年前
PhoneGap设置Icon
参考:http://cordova.apache.org/docs/en/latest/config\_ref/images.html通过config.xml中的<icon标签来设置Icon<iconsrc"res/ios/icon.png"platform"ios"width"57"height"57"densi
Wesley13 Wesley13
3年前
mysql设置时区
mysql设置时区mysql\_query("SETtime\_zone'8:00'")ordie('时区设置失败,请联系管理员!');中国在东8区所以加8方法二:selectcount(user\_id)asdevice,CONVERT\_TZ(FROM\_UNIXTIME(reg\_time),'08:00','0
Wesley13 Wesley13
3年前
00:Java简单了解
浅谈Java之概述Java是SUN(StanfordUniversityNetwork),斯坦福大学网络公司)1995年推出的一门高级编程语言。Java是一种面向Internet的编程语言。随着Java技术在web方面的不断成熟,已经成为Web应用程序的首选开发语言。Java是简单易学,完全面向对象,安全可靠,与平台无关的编程语言。
Stella981 Stella981
3年前
Django中Admin中的一些参数配置
设置在列表中显示的字段,id为django模型默认的主键list_display('id','name','sex','profession','email','qq','phone','status','create_time')设置在列表可编辑字段list_editable
Python进阶者 Python进阶者
1年前
Excel中这日期老是出来00:00:00,怎么用Pandas把这个去除
大家好,我是皮皮。一、前言前几天在Python白银交流群【上海新年人】问了一个Pandas数据筛选的问题。问题如下:这日期老是出来00:00:00,怎么把这个去除。二、实现过程后来【论草莓如何成为冻干莓】给了一个思路和代码如下:pd.toexcel之前把这