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

原文出處：拓端數(shù)據(jù)部落公眾號(hào)

馬爾可夫鏈?zhǔn)菑囊粋€(gè)“狀態(tài)”（一種情況或一組值）跳到另一個(gè)“狀態(tài)”的數(shù)學(xué)系統(tǒng)。本文介紹了馬爾可夫鏈和一種簡(jiǎn)單的狀態(tài)轉(zhuǎn)移模型，該模型構(gòu)成了隱馬爾可夫模型（HMM）的特例。從應(yīng)用的角度來看，這些模型在評(píng)估經(jīng)濟(jì)/市場(chǎng)狀態(tài)數(shù)據(jù)時(shí)（查看文末了解數(shù)據(jù)獲取方式）非常有用。這里的討論主要圍繞使用這些模型的科學(xué)性。

視頻：馬爾可夫鏈原理可視化解釋與R語言區(qū)制轉(zhuǎn)換Markov regime switching實(shí)例

馬爾可夫鏈原理可視化解釋與R語言區(qū)制轉(zhuǎn)換Markov?regime?switching實(shí)例

，時(shí)長(zhǎng)07:25

例如，如果您制作了嬰兒行為的馬爾可夫鏈模型，您可能會(huì)將“玩?！薄ⅰ俺燥垺?、“睡覺”和“哭泣”作為狀態(tài)，它們與其他行為一起可以形成“狀態(tài)空間”：所有可能狀態(tài)的列表。此外，在狀態(tài)空間之上，馬爾可夫鏈告訴您從一個(gè)狀態(tài)跳躍或“轉(zhuǎn)換”到任何其他狀態(tài)的概率——例如，正在玩耍的嬰兒在下一個(gè)狀態(tài)下入睡的可能性五分鐘不先哭。

一個(gè)簡(jiǎn)單的兩態(tài)馬爾可夫鏈如下所示。

在我們的狀態(tài)空間中有兩種狀態(tài)（A 和 B），有 4 種可能的轉(zhuǎn)換（不是 2 種，因?yàn)闋顟B(tài)可以轉(zhuǎn)換回自身）。如果我們?cè)凇癆”，我們可以過渡到“B”或留在“A”。如果我們?cè)凇癇”，我們可以過渡到“A”或留在“B”。在這兩個(gè)狀態(tài)圖中，從任何狀態(tài)轉(zhuǎn)換到任何其他狀態(tài)的概率為 0.5。

當(dāng)然，真正的建模者并不總是畫出馬爾可夫鏈圖。相反，他們使用“轉(zhuǎn)移矩陣”來計(jì)算轉(zhuǎn)移概率。狀態(tài)空間中的每個(gè)狀態(tài)都包含一次作為行和列，并且矩陣中的每個(gè)單元格都告訴您從其行狀態(tài)轉(zhuǎn)換到其列狀態(tài)的概率。因此，在矩陣中，單元格的作用與圖中箭頭的作用相同。

如果狀態(tài)空間添加一個(gè)狀態(tài)，我們添加一行和一列，每個(gè)現(xiàn)有的列和行添加一個(gè)單元格。這意味著當(dāng)我們向馬爾可夫鏈添加狀態(tài)時(shí)，單元格的數(shù)量呈二次增長(zhǎng)。因此，除非您想繪制叢林健身房馬爾可夫鏈圖，否則轉(zhuǎn)換矩陣很快就會(huì)派上用場(chǎng)。

馬爾可夫鏈的一種用途是在計(jì)算機(jī)模擬中包含真實(shí)世界的現(xiàn)象。例如，我們可能想要檢查新大壩的溢出頻率，這取決于連續(xù)下雨的天數(shù)。為了建立這個(gè)模型，我們從以下雨天?(R) 和晴天 (S) 開始：

模擬這種天氣的一種方法是只說“有一半的日子下雨。因此，在我們的模擬中，每天都有 50% 的幾率下雨。”

該規(guī)則將在模擬中生成以下序列：

你有沒有注意到上面的序列看起來和原來的不太一樣？第二個(gè)序列似乎跳來跳去，而第一個(gè)（真實(shí)數(shù)據(jù)）似乎具有“粘性”。在實(shí)際數(shù)據(jù)中，如果一天是晴天（S），那么第二天晴天的可能性也更大。

我們可以使用兩態(tài)馬爾可夫鏈來縮小這種“粘性”。當(dāng)馬爾可夫鏈處于“R”狀態(tài)時(shí)，它有 0.9 的概率留在原地，有 0.1 的機(jī)會(huì)離開“S”狀態(tài)。同樣，“S”狀態(tài)有 0.9 的概率保持原狀，并且有 0.1 的機(jī)會(huì)轉(zhuǎn)換到“R”狀態(tài)。

在氣象學(xué)家、生態(tài)學(xué)家、計(jì)算機(jī)科學(xué)家、金融工程師和其他需要對(duì)大現(xiàn)象建模的人群中，馬爾可夫鏈可以變得非常龐大和強(qiáng)大。例如，谷歌用來確定搜索結(jié)果順序的算法，稱為PageRank，是一種馬爾可夫鏈。

馬爾可夫區(qū)制轉(zhuǎn)移模型Markov regime switching

本文簡(jiǎn)要介紹了一種簡(jiǎn)單的狀態(tài)轉(zhuǎn)移模型，該模型構(gòu)成了隱馬爾可夫模型（HMM）的特例。這些模型擬合時(shí)間序列數(shù)據(jù)中的非平穩(wěn)性。

基本案例

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

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

為了演示，我們準(zhǔn)備一些數(shù)據(jù)并嘗試進(jìn)行反向推測(cè)。每個(gè)狀態(tài)都具有不同的均值和波動(dòng)率。

theta_v?<-?data.frame(t(c(2.00,-2.00,1.00,2.00,0.95,0.85))) kable(theta_v,?"html",?booktabs?=?F,escape?=?F)?%>%? ????????kable_styling(position?=?"center")

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

馬爾可夫過程

為了模擬過程x\_t ，我們從模擬馬爾可夫過程s\_t 開始。為了模擬T 期間的過程，首先，我們需要構(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。

P?<-?matrix(c(p11,1-p22,1-p11,p22),2,2) P\[1,\]##?\[1\]?0.95?0.05

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

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

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

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

因此，有15％的概率處于1狀態(tài)，而有25％的概率處于狀態(tài)2。這應(yīng)該反映在模擬過程中??s，從而

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

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

結(jié)果

給定模擬的馬爾可夫過程，結(jié)果的模擬非常簡(jiǎn)單。一個(gè)簡(jiǎn)單的技巧是模擬

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

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

在本節(jié)中，我將使用R軟件手動(dòng)（從頭開始）和非手動(dòng)進(jìn)行統(tǒng)計(jì)分解。在前者中，我將演示如何構(gòu)造似然函數(shù)，然后使用約束優(yōu)化問題來估計(jì)參數(shù)。

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

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

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

顯然，這兩種狀態(tài)的每次過濾器的總和應(yīng)為1?？雌饋?，我們可以處于狀態(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)

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

手動(dòng)估算

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

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

在第二步中，我為估算設(shè)置了約束

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

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

最后，回想一下，通過構(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ì)值非常一致。

估算?

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

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

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

平均而言，我們應(yīng)該期望過程平均值約為1，即

。這是由期望定律得出的，其中我們知道

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

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

##?\[1\]?2.179449

我們注意到，回歸估計(jì)值與波動(dòng)率的一致性高于均值。

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

為了揭示這些模式，我們?cè)谙旅嫜菔救绾问褂蒙厦娴木€性模型建立區(qū)制轉(zhuǎn)移模型：

主要輸入是擬合模型，??mod我們將其歸納為擬合轉(zhuǎn)移狀態(tài)。第二個(gè)??k是區(qū)制的數(shù)量。由于我們知道我們要處理兩個(gè)狀態(tài)，因此將其設(shè)置為2。但是，實(shí)際上，需要參考一種信息標(biāo)準(zhǔn)來確定最佳狀態(tài)數(shù)。根據(jù)定義，我們有兩個(gè)參數(shù)，均值 mu\_s 和波動(dòng)率 sigma\_s 。因此，我們添加一個(gè)true / false向量來指示正在轉(zhuǎn)移的參數(shù)。在上面的命令中，我們?cè)试S兩個(gè)參數(shù)都轉(zhuǎn)移。最后，我們可以指定估計(jì)過程是否正在使用并行計(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)的均值和波動(dòng)。模型1的平均值為1.71，波動(dòng)率接近1。模型2的平均值為-2，波動(dòng)率約為2。顯然，該模型針對(duì)數(shù)據(jù)確定了兩種具有不同均值和波動(dòng)率的不同狀態(tài)。其次，在輸出的底部，擬合的模型報(bào)告了轉(zhuǎn)移概率。

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

par(mar?=?2*c(1,1,1,1),mfrow?=?c(2,1))

頂部的圖表示隨時(shí)間變化的過程 x_t ，其中灰色陰影區(qū)域表示

的時(shí)間段。換句話說，灰色區(qū)域表示狀態(tài)1占優(yōu)勢(shì)的時(shí)間段。

plot(x~t_index,type?="l",col?=?0,xlim=c(1,100))

過濾器會(huì)在一個(gè)周期內(nèi)檢測(cè)到第二種狀態(tài)。發(fā)生這種情況是因?yàn)樵谶@種情況下，返回的是平滑概率，即在實(shí)現(xiàn)整個(gè)樣本 T ?后處于每種狀態(tài)的概率，即

。另一方面，來自手動(dòng)估計(jì)的推斷概率。

無論如何，由于我們知道狀態(tài)的真實(shí)值，因此可以確定我們是否處于真實(shí)狀態(tài)。我們?cè)谏厦娴膱D中使用黑點(diǎn)突出顯示狀態(tài)2。總的來說，我們觀察到模型在檢測(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é)束語

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

數(shù)據(jù)獲取

在下面公眾號(hào)后臺(tái)回復(fù)“MRS數(shù)據(jù)”，可獲取完整數(shù)據(jù)。

本文摘選《R語言馬爾可夫區(qū)制轉(zhuǎn)移模型Markov regime switching》，點(diǎn)擊“閱讀原文”獲取全文完整資料。

點(diǎn)擊標(biāo)題查閱往期內(nèi)容

時(shí)變馬爾可夫區(qū)制轉(zhuǎn)換MRS自回歸模型分析經(jīng)濟(jì)時(shí)間序列

馬爾可夫轉(zhuǎn)換模型研究交通傷亡人數(shù)事故時(shí)間序列預(yù)測(cè)

如何實(shí)現(xiàn)馬爾可夫鏈蒙特卡羅MCMC模型、Metropolis算法？

Matlab用BUGS馬爾可夫區(qū)制轉(zhuǎn)換Markov switching隨機(jī)波動(dòng)率模型、序列蒙特卡羅SMC、M H采樣分析時(shí)間序列

matlab用馬爾可夫鏈蒙特卡羅 (MCMC) 的Logistic邏輯回歸模型分析汽車實(shí)驗(yàn)數(shù)據(jù)

R語言BUGS序列蒙特卡羅SMC、馬爾可夫轉(zhuǎn)換隨機(jī)波動(dòng)率SV模型、粒子濾波、Metropolis Hasting采樣時(shí)間序列分析

stata馬爾可夫Markov區(qū)制轉(zhuǎn)移模型分析基金利率

PYTHON用時(shí)變馬爾可夫區(qū)制轉(zhuǎn)換（MRS）自回歸模型分析經(jīng)濟(jì)時(shí)間序列

R語言使用馬爾可夫鏈對(duì)營銷中的渠道歸因建模

matlab實(shí)現(xiàn)MCMC的馬爾可夫轉(zhuǎn)換ARMA - GARCH模型估計(jì)

R語言隱馬爾可夫模型HMM識(shí)別不斷變化的股票市場(chǎng)條件

R語言中的隱馬爾可夫HMM模型實(shí)例

用機(jī)器學(xué)習(xí)識(shí)別不斷變化的股市狀況—隱馬爾科夫模型(HMM)

Matlab馬爾可夫鏈蒙特卡羅法（MCMC）估計(jì)隨機(jī)波動(dòng)率（SV，Stochastic Volatility） 模型

MATLAB中的馬爾可夫區(qū)制轉(zhuǎn)移(Markov regime switching)模型

Matlab馬爾可夫區(qū)制轉(zhuǎn)換動(dòng)態(tài)回歸模型估計(jì)GDP增長(zhǎng)率

R語言馬爾可夫區(qū)制轉(zhuǎn)移模型Markov regime switching

stata馬爾可夫Markov區(qū)制轉(zhuǎn)移模型分析基金利率

R語言如何做馬爾可夫轉(zhuǎn)換模型markov switching model

R語言隱馬爾可夫模型HMM識(shí)別股市變化分析報(bào)告

R語言中實(shí)現(xiàn)馬爾可夫鏈蒙特卡羅MCMC模型