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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

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

基本案例

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

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

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

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％的概率進入狀態(tài)2。

P <- matrix(c(p11,1-p22,1-p11,p22),2,2)
P\[1,\]

## \[1\] 0.95 0.05

因為它先前的狀態(tài)，模擬s_t 是遞歸的。因此，我們需要構(gòu)造一個循環(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的“實現(xiàn)”多于狀態(tài)2。實際上，這可以由穩(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個周期的小樣本，因此我們觀察到穩(wěn)定概率為69％，接近但不完全等于75％。

結(jié)果

給定模擬的馬爾可夫過程，結(jié)果的模擬非常簡單。一個簡單的技巧是模擬的T周期和的 T  周期。然后，給定 s\_t  的模擬，我們針對每個狀態(tài)創(chuàng)建結(jié)果變量 x\_t  。

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

---

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

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

左右滑動查看更多

01

02

03

04

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

估計參數(shù)

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

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

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

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

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

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

## \[1\] TRUE

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

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

手動估算

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

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

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

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

all(A%*%theta0 >= B)

## \[1\] TRUE

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

## $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ù)一致，我們繪制估計值與真實值的關(guān)系圖：

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

總體而言，我們觀察到估計值非常一致。

估算 

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

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

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

對于波動率，適用相同的邏輯。

## \[1\] 2.179449

我們注意到，回歸估計值與波動率的一致性高于均值。

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

為了揭示這些模式，我們在下面演示如何使用上面的線性模型建立區(qū)制轉(zhuǎn)移模型：

主要輸入是擬合模型，  mod我們將其歸納為擬合轉(zhuǎn)移狀態(tài)。第二個  k是區(qū)制的數(shù)量。由于我們知道我們要處理兩個狀態(tài)，因此將其設(shè)置為2。但是，實際上，需要參考一種信息標準來確定最佳狀態(tài)數(shù)。根據(jù)定義，我們有兩個參數(shù)，均值 mu\_s 和波動率 sigma\_s 。因此，我們添加一個true / false向量來指示正在轉(zhuǎn)移的參數(shù)。在上面的命令中，我們允許兩個參數(shù)都轉(zhuǎ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

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

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

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

頂部的圖表示隨時間變化的過程 x_t ，其中灰色陰影區(qū)域表示的時間段。換句話說，灰色區(qū)域表示狀態(tài)1占優(yōu)勢的時間段。

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

過濾器會在一個周期內(nèi)檢測到第二種狀態(tài)。發(fā)生這種情況是因為在這種情況下，返回的是平滑概率，即在實現(xiàn)整個樣本 T  后處于每種狀態(tài)的概率，即。另一方面，來自手動估計的推斷概率。

無論如何，由于我們知道狀態(tài)的真實值，因此可以確定我們是否處于真實狀態(tài)。我們在上面的圖中使用黑點突出顯示狀態(tài)2?？偟膩碚f，我們觀察到模型在檢測數(shù)據(jù)狀態(tài)方面表現(xiàn)非常好。唯一的例外是第60天，其中推斷概率大于50％。要查看推斷概率多長時間正確一次，我們運行以下命令

mean(Filter$Regime_1 == (s==1)*1)

## \[1\] 0.96

結(jié)束語

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