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

最近我們被客戶要求撰寫關(guān)于馬爾可夫轉(zhuǎn)換模型的研究報告，包括一些圖形和統(tǒng)計輸出。

假設(shè) 有時間序列數(shù)據(jù)，如下所示。經(jīng)驗表明，目標變量y似乎與解釋變量x有關(guān)。然而，乍一看，y在水平中間波動，所以它似乎并不總是有穩(wěn)定的關(guān)系（背后有多個狀態(tài)）

上面的樣本數(shù)據(jù)創(chuàng)建如下。x和y之間的關(guān)系數(shù)據(jù)根據(jù)時間改變。

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之間的關(guān)系如下圖所示。

數(shù)據(jù)

?

在馬爾可夫轉(zhuǎn)換模型中，觀察數(shù)據(jù)被認為是從幾個狀態(tài)生成的，并且如上所示可以很好地分離。

觀察到的數(shù)據(jù)

點擊標題查閱往期內(nèi)容

【視頻】馬爾可夫鏈蒙特卡羅方法MCMC原理與R語言實現(xiàn)|數(shù)據(jù)分享

左右滑動查看更多

01

02

03

04

創(chuàng)建馬爾可夫轉(zhuǎn)換模型

?

模型公式?

?

# 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

?

?參數(shù)的含義是

• k：馬爾可夫轉(zhuǎn)換模型的狀態(tài)數(shù)。在這里，它被指定為后面有兩個狀態(tài)。

• sw：指定每個參數(shù)在狀態(tài)更改時是否更改

• p：AR模型系數(shù)

• family:(在GLM的情況下）概率分布族

# 馬爾可夫轉(zhuǎn)換模型# # ? ? ? ?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

輸出中的區(qū)制1和區(qū)制2表示模型的兩個狀態(tài) 。?

# 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 ?

可以看到區(qū)制2?與y1 <- x * 4 + 20匹配。

從調(diào)整后的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 ? ?

?模型

對于每個狀態(tài)，處于該狀態(tài)的概率以陰影繪制

每個時間點的概率?

每次獲取狀態(tài)和更改點

如果你想知道你在某個特定時間點所在的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...

異常值/變化點是狀態(tài)更改的時間?

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

因此，我們可以看到檢測到在第一次數(shù)據(jù)創(chuàng)建時指定的變化點。

點擊文末?“閱讀原文”

獲取全文完整代碼數(shù)據(jù)資料。

本文選自《R語言如何做馬爾可夫轉(zhuǎn)換模型markov switching model》。

點擊標題查閱往期內(nèi)容

matlab用馬爾可夫鏈蒙特卡羅 (MCMC) 的Logistic邏輯回歸模型分析汽車實驗數(shù)據(jù)
【視頻】馬爾可夫鏈蒙特卡羅方法MCMC原理與R語言實現(xiàn)|數(shù)據(jù)分享
R語言BUGS/JAGS貝葉斯分析: 馬爾科夫鏈蒙特卡洛方法（MCMC）采樣
馬爾可夫Markov區(qū)制轉(zhuǎn)移模型分析基金利率
馬爾可夫區(qū)制轉(zhuǎn)移模型Markov regime switching
時變馬爾可夫區(qū)制轉(zhuǎn)換MRS自回歸模型分析經(jīng)濟時間序列
馬爾可夫轉(zhuǎn)換模型研究交通傷亡人數(shù)事故時間序列預(yù)測
如何實現(xiàn)馬爾可夫鏈蒙特卡羅MCMC模型、Metropolis算法？
Matlab用BUGS馬爾可夫區(qū)制轉(zhuǎn)換Markov switching隨機波動率模型、序列蒙特卡羅SMC、M H采樣分析時間序列
R語言BUGS序列蒙特卡羅SMC、馬爾可夫轉(zhuǎn)換隨機波動率SV模型、粒子濾波、Metropolis Hasting采樣時間序列分析
matlab用馬爾可夫鏈蒙特卡羅 (MCMC) 的Logistic邏輯回歸模型分析汽車實驗數(shù)據(jù)
stata馬爾可夫Markov區(qū)制轉(zhuǎn)移模型分析基金利率
PYTHON用時變馬爾可夫區(qū)制轉(zhuǎn)換（MRS）自回歸模型分析經(jīng)濟時間序列
R語言使用馬爾可夫鏈對營銷中的渠道歸因建模
matlab實現(xiàn)MCMC的馬爾可夫轉(zhuǎn)換ARMA - GARCH模型估計
R語言隱馬爾可夫模型HMM識別不斷變化的股票市場條件
R語言中的隱馬爾可夫HMM模型實例
用機器學(xué)習(xí)識別不斷變化的股市狀況—隱馬爾科夫模型(HMM)
Matlab馬爾可夫鏈蒙特卡羅法（MCMC）估計隨機波動率（SV，Stochastic Volatility） 模型
MATLAB中的馬爾可夫區(qū)制轉(zhuǎn)移(Markov regime switching)模型
Matlab馬爾可夫區(qū)制轉(zhuǎn)換動態(tài)回歸模型估計GDP增長率
R語言馬爾可夫區(qū)制轉(zhuǎn)移模型Markov regime switching
stata馬爾可夫Markov區(qū)制轉(zhuǎn)移模型分析基金利率
R語言如何做馬爾可夫轉(zhuǎn)換模型markov switching model
R語言隱馬爾可夫模型HMM識別股市變化分析報告
R語言中實現(xiàn)馬爾可夫鏈蒙特卡羅MCMC模型