【原】R語(yǔ)言馬爾可夫體制轉(zhuǎn)換模型Markov regime switching

拓端數(shù)據(jù) 2020-05-08

展開全文

原文鏈接：http:///?p=12280

總覽

本文簡(jiǎn)要介紹了一種簡(jiǎn)單的狀態(tài)切換模型，該模型構(gòu)成了隱馬爾可夫模型（HMM）的特例。這些模型適應(yīng)時(shí)間序列數(shù)據(jù)中的非平穩(wěn)性。從應(yīng)用的角度來(lái)看，這些模型在評(píng)估經(jīng)濟(jì)/市場(chǎng)狀態(tài)時(shí)非常有用。這里的討論主要圍繞使用這些模型的科學(xué)性。

基本案例

HMM的主要挑戰(zhàn)是預(yù)測(cè)隱藏部分。我們?nèi)绾巫R(shí)別“不可觀察”的事物？HMM的想法是從可觀察的事物來(lái)預(yù)測(cè)潛在的事物。

模擬數(shù)據(jù)

為了演示，我們準(zhǔn)備一些數(shù)據(jù)并嘗試進(jìn)行反向推測(cè)。通過(guò)構(gòu)造，我強(qiáng)加了一些假設(shè)來(lái)創(chuàng)建我們的數(shù)據(jù)。每個(gè)狀態(tài)都具有不同的均值和波動(dòng)率。

library(knitr)library(kableExtra)library(dplyr)theta_v <- data.frame(t(c(2.00,-2.00,1.00,2.00,0.95,0.85)))names(theta_v) <- c("$\\mu_1$","$\\mu_2$","$\\sigma_1$","$\\sigma_2$","$p_{11}$","$p_{22}$")kable(theta_v, "html", booktabs = F,escape = F) %>% kable_styling(position = "center")

mu_1	mu_2	sigma_1	sigma_2	p_ {11}	p_ {22}
2	-2	1個(gè)	2	0.95	0.85

如上表所示，狀態(tài)\（s = 2 \）變成“壞”狀態(tài)，其中過(guò)程\（x_t \）表現(xiàn)出較高的變化性。因此，停留在狀態(tài)2的可能比停留在狀態(tài)1的可能性小。

馬爾可夫過(guò)程

為了模擬過(guò)程\（x_t \），我們從模擬馬爾可夫過(guò)程\（s_t \）開始。為了模擬\（T \）期間的過(guò)程，首先，我們需要構(gòu)建給定\（p_ {11} \）和\（p_ {22} \）的轉(zhuǎn)換矩陣。其次，我們需要從給定狀態(tài)\（s_1 = 1 \）開始。從\（s_1 = 1 \）開始，我們知道有95％的概率停留在狀態(tài)1，有5％的概率進(jìn)入狀態(tài)2。

p11 <- theta_v[1,5] p22 <- theta_v[1,6]P <- matrix(c(p11,1-p22,1-p11,p22),2,2)P[1,]## [1] 0.95 0.05

模擬\（s_t \）是遞歸的，因?yàn)樗惹跋惹暗臓顟B(tài)。因此，我們需要構(gòu)造一個(gè)循環(huán)：

set.seed(13)T_end <- 10^2

s0 <- 1 st <- function(i) sample(1:2,1,prob = P[i,])

s <- st(s0)for(t in 2:T_end) { s <- c(s,st(s[t-1]))}plot(s, pch = 20,cex = 0.5)

上圖說(shuō)明了過(guò)程\（s_t \）的持久性。在大多數(shù)情況下，狀態(tài)1的“實(shí)現(xiàn)”多于狀態(tài)2。實(shí)際上，這可以由固定概率確定，該固定概率由下式表示：

P_stat[1,]## [1] 0.75 0.25

因此，經(jīng)歷時(shí)間的流逝，有15％的概率處于1狀態(tài)，而有25％的概率處于狀態(tài)2。這應(yīng)該反映在模擬過(guò)程中 s，從而

mean(s==1)## [1] 0.69

由于我們使用的是100個(gè)周期的小樣本，因此我們觀察到穩(wěn)定概率為69％，接近但不完全等于75％。

結(jié)果

給定模擬的馬爾可夫過(guò)程，結(jié)果過(guò)程的模擬非常簡(jiǎn)單。一個(gè)簡(jiǎn)單的技巧是模擬\（x_t \ mid s_t = 1 \）的\（T \）周期和\ （x_t \ mid s_t = 2 \）的\（T \）周期。然后，給定\（s_t \）的模擬，我們針對(duì)每個(gè)狀態(tài)創(chuàng)建結(jié)果變量\（x_t \）。

plot(x~t_index, pch = 20)points(x[s == 2]~t_index[s==2],col = 2)

雖然總體而言時(shí)間序列看起來(lái)是平穩(wěn)的，但我們觀察到一些周期（以紅色突出顯示）顯示出較高的波動(dòng)。有人可能會(huì)建議說(shuō)，數(shù)據(jù)存在結(jié)構(gòu)性中斷，或者體制發(fā)生了變化，過(guò)程\（x_t \）變得越來(lái)越大，帶有更多的負(fù)值。雖然如此，事后解釋總是比較容易的。主要的挑戰(zhàn)是識(shí)別這種情況。

估計(jì)參數(shù)

在本節(jié)中，我將使用R軟件手動(dòng)（從頭開始）和非手動(dòng)進(jìn)行統(tǒng)計(jì)分解。在前者中，我將演示如何構(gòu)造似然函數(shù)，然后使用約束優(yōu)化問(wèn)題來(lái)估計(jì)參數(shù)。我將說(shuō)明如何在不經(jīng)歷解析推導(dǎo)的情況下進(jìn)行復(fù)制。

似然函數(shù)-數(shù)值部分

首先，我們需要?jiǎng)?chuàng)建一個(gè)以\（\ Theta \）向量為主要輸入的函數(shù)。其次，我們需要設(shè)置一個(gè)返回MLE的優(yōu)化問(wèn)題。

在優(yōu)化似然函數(shù)之前。讓我們看一下工作原理。假設(shè)我們知道參數(shù)\（\ Theta \）的向量，并且我們有興趣使用\（x_t \）上的數(shù)據(jù)評(píng)估隱藏狀態(tài)隨時(shí)間的變化。

t	xi_ {t mid t，1}	xi_ {t mid t，2}
1	0.878	0.122
2	0.982	0.018
3	0.887	0.113
4	0.875	0.125
5	0.318	0.682
6	0.000	1.000

顯然，這兩種狀態(tài)的每次過(guò)濾器的總和應(yīng)為1?？雌饋?lái)，我們可以處于狀態(tài)1或狀態(tài)2。

all(round(apply(Filter[,-1],1,sum),9) == 1)## [1] TRUE

由于我們?cè)O(shè)計(jì)了此數(shù)據(jù)，因此我們知道狀態(tài)2的時(shí)期。

plot(Filter[,3]~t_index, type = "l", ylab = expression(xi[2]))points(Filter[s==2,3]~t_index[s==2],pch = 20, col = 2)

過(guò)濾器背后的優(yōu)點(diǎn)是僅使用\（x_t \）上的信息來(lái)識(shí)別潛在狀態(tài)。我們觀察到，狀態(tài)2的過(guò)濾器主要在狀態(tài)2發(fā)生時(shí)增加。這可以通過(guò)發(fā)出紅點(diǎn)的概率增加來(lái)證明，紅點(diǎn)表示狀態(tài)2發(fā)生的時(shí)間段。盡管如此，上述還是存在一些重大問(wèn)題。首先，它假定我們知道參數(shù)\（\ Theta \），而實(shí)際上我們需要對(duì)此進(jìn)行估計(jì)，然后在此基礎(chǔ)上進(jìn)行推斷。其次，所有這些都是在樣本中構(gòu)造的。從實(shí)際的角度來(lái)看，決策者對(duì)預(yù)測(cè)的概率及其對(duì)未來(lái)投資的影響感興趣。

手動(dòng)估算

為了優(yōu)化上面定義的 HMM_Lik 功能，我將需要執(zhí)行兩個(gè)附加步驟。首先是建立一個(gè)初始估計(jì)值，作為搜索算法的起點(diǎn)。其次，我們需要設(shè)置約束條件以驗(yàn)證估計(jì)的參數(shù)是否一致，即非負(fù)波動(dòng)性和介于0和1之間的概率值。

第一步，我使用樣本創(chuàng)建初始參數(shù)向量\（\ Theta_0 \）

在第二步中，我為估算制定了約束

請(qǐng)注意，參數(shù)的初始向量應(yīng)滿足約束條件

all(A%*%theta0 >= B)## [1] TRUE

最后，回想一下，通過(guò)構(gòu)建大多數(shù)優(yōu)化算法都可以搜索最小點(diǎn)。因此，我們需要將似然函數(shù)的輸出更改為負(fù)值。

## $par## [1] 1.7119528 -1.9981224 0.8345350 2.2183230 0.9365507 0.8487511## ## $value## [1] 174.7445## ## $counts## function gradient ## 1002 NA ## ## $convergence## [1] 0## ## $message## NULL## ## $outer.iterations## [1] 3## ## $barrier.value## [1] 6.538295e-05

為了檢查MLE值是否與真實(shí)參數(shù)一致，我們繪制估計(jì)值與真實(shí)值的關(guān)系圖：

plot(opt$par ~ theta_known,pch = 20,cex=2,ylab="MLE",xlab = "True")abline(a=0,b=1,lty=2)

總體而言，我們觀察到估計(jì)值非常一致，由于MLE的一致性屬性，這不足為奇。

估算

我將在下面演示如何使用r軟件復(fù)制人工估算的結(jié)果。

如果我們要忽略過(guò)程中的任何體制轉(zhuǎn)換，我們可以簡(jiǎn)單地將參數(shù)\（\ mu \）和\（\ sigma \）估計(jì)為

kable(mod_est, "html", booktabs = F,escape = F) %>% kable_styling(position = "center")

hat {mu}	hat {sigma}
0.6244574	2.198929

平均而言，我們應(yīng)該期望過(guò)程平均值約為1，即\（0.75 \ times 2 + 0.25 \ times（-1）= 1 \）。這是由總期望屬性定律得出的，其中我們知道\ [\ begin {equation} \ mathbb {E} [x] = \ sum _ {\ forall s} \ mathbb {E} \ left [x \ mid s \右] \ mathbb {P}（s）= \ sum _ {\ forall s} \ mu_ {s} \ pi_ {s} \ end {equation} \]

這樣

EX <- 0.75*2 + 0.25*-2EX## [1] 1

對(duì)于波動(dòng)率，適用相同的邏輯。

EX2 <- (2^2 + 1^2)*0.75 + ((-2)^2 + 2^2)*0.25VX <- EX2 - EX^2sqrt(VX)## [1] 2.179449

我們注意到，回歸估計(jì)值與波動(dòng)率的一致性高于均值。這主要是由于估算第一時(shí)刻與第二時(shí)刻的工作比較繁瑣。

上面的觀點(diǎn)是，估計(jì)值并未涵蓋數(shù)據(jù)的真實(shí)性質(zhì)。如果我們假設(shè)數(shù)據(jù)是固定的，那么我們錯(cuò)誤地估計(jì)過(guò)程的平均值為62％。但是，與此同時(shí)，我們通過(guò)構(gòu)造知道該過(guò)程表現(xiàn)出兩個(gè)平均結(jié)果-一個(gè)正面和一個(gè)負(fù)面。波動(dòng)性也是如此。

為了揭示這些模式，我們?cè)谙旅嫜菔救绾问褂蒙厦娴木€性模型部署狀態(tài)切換模型：

主要輸入是擬合模型， mod我們將其歸納為適應(yīng)切換狀態(tài)。第二個(gè) k是制度的數(shù)量。由于我們知道我們要處理兩個(gè)狀態(tài)，因此將其設(shè)置為2。但是，實(shí)際上，需要參考一種信息標(biāo)準(zhǔn)來(lái)確定最佳狀態(tài)數(shù)。根據(jù)定義，我們有兩個(gè)參數(shù)，均值\（\ mu_s \）和波動(dòng)率\（\ sigma_s \）。因此，我們添加一個(gè)true / false向量來(lái)指示正在切換的參數(shù)。在上面的命令中，我們?cè)试S兩個(gè)參數(shù)都切換。最后，我們可以指定估計(jì)過(guò)程是否正在使用并行計(jì)算進(jìn)行。

要了解模型的輸出，讓我們看一下

## Markov Switching Model## ## ## AIC BIC logLik## 352.2843 366.705 -174.1422## ## Coefficients:## (Intercept)(S) Std(S)## Model 1 1.711693 0.8346013## Model 2 -2.004137 2.2155742## ## Transition probabilities:## Regime 1 Regime 2## Regime 1 0.93767719 0.1510052## Regime 2 0.06232281 0.8489948

上面的輸出主要報(bào)告我們嘗試手動(dòng)估算的六個(gè)估算參數(shù)。首先，系數(shù)表報(bào)告了每個(gè)狀態(tài)的均值和波幅。模型1的平均值為1.71，波動(dòng)率接近1。模型2的平均值為-2，波動(dòng)率約為2。顯然，該模型針對(duì)數(shù)據(jù)確定了兩種具有不同均值和波動(dòng)率的不同狀態(tài)。其次，在輸出的底部，擬合的模型報(bào)告過(guò)渡概率。

plot(opt$par ~ theta_known,pch = 20,cex=2,ylab="MLE",xlab = "True")points(theta_mswm~theta_known,pch = 1,col = 2, cex = 2,lwd = 2)abline(a=0,b=1,lty=2)legend("topleft",c("Manual","MSwM"), pch = c(20,1), col = 1:2)

有趣的是，就每種狀態(tài)的過(guò)濾器而言，我們將從包中檢索到的狀態(tài)與手動(dòng)提取的狀態(tài)進(jìn)行比較。根據(jù)定義，可以使用mod.mswm 對(duì)象上的圖函數(shù) 來(lái)了解平滑概率以及確定的方案。

par(mar = 2*c(1,1,1,1),mfrow = c(2,1))plotProb(mod.mswm,2)

頂部的圖表示隨時(shí)間變化的過(guò)程\（x_t \），其中灰色陰影區(qū)域表示\（\ hat {\ xi} _ {T \ mid t，1}> 0.5 \）的時(shí)間段。換句話說(shuō)，灰色區(qū)域表示狀態(tài)1占優(yōu)勢(shì)的時(shí)間段。

plot(x~t_index,type ="l",col = 0,xlim=c(1,100))rect(xx-1,-10,xx,10,col = "lightgray",lty = 0)lines(x~t_index)points(x[s==2] ~ t_index[s==2],col = 1,pch = 20)

過(guò)濾器會(huì)在一個(gè)周期內(nèi)檢測(cè)到第二種狀態(tài)。發(fā)生這種情況是因?yàn)?nbsp; plotProb 在這種情況下，返回的是平滑概率，即在實(shí)現(xiàn)整個(gè)樣本\（T \）后處于每種狀態(tài)的概率，即\（\ hat {\ xi} _ {t \ mid T，1} \ ）。另一方面，來(lái)自手動(dòng)估計(jì)的一個(gè)對(duì)應(yīng)于推斷概率\（\ hat {\ xi} _ {t \ mid t，1} \）。

無(wú)論如何，由于我們知道狀態(tài)的真實(shí)值，因此可以確定我們是否處于真實(shí)狀態(tài)。我們?cè)谏厦娴膱D中使用黑點(diǎn)突出顯示狀態(tài)2。總的來(lái)說(shuō)，我們觀察到模型在檢測(cè)數(shù)據(jù)狀態(tài)方面表現(xiàn)非常出色。唯一的例外是第60天，其中推斷概率大于50％。要查看推理概率多長(zhǎng)時(shí)間正確一次，我們運(yùn)行以下命令

mean(Filter$Regime_1 == (s==1)*1)## [1] 0.96

結(jié)束語(yǔ)

此過(guò)程似乎運(yùn)行良好。然而，在實(shí)際數(shù)據(jù)實(shí)現(xiàn)方面仍然存在許多挑戰(zhàn)。首先，我們不具備有關(guān)數(shù)據(jù)生成過(guò)程的知識(shí)。其次，狀態(tài)不一定實(shí)現(xiàn)。因此，這兩個(gè)問(wèn)題可能會(huì)破壞狀態(tài)切換模式的可靠性。在應(yīng)用方面，通常部署此類模型以評(píng)估經(jīng)濟(jì)或市場(chǎng)狀況。從決策上來(lái)說(shuō)，這也可以為策略分配提供有趣的應(yīng)用。