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

最近我們被客戶要求撰寫關(guān)于狀態(tài)空間模型的研究報(bào)告，包括一些圖形和統(tǒng)計(jì)輸出。

狀態(tài)空間建模是一種高效、靈活的方法，用于對大量的時(shí)間序列和其他數(shù)據(jù)進(jìn)行統(tǒng)計(jì)推斷

摘要

本文介紹了狀態(tài)空間建模，其觀測值來自指數(shù)族，即高斯、泊松、二項(xiàng)、負(fù)二項(xiàng)和伽馬分布。在介紹了高斯和非高斯?fàn)顟B(tài)空間模型的基本理論后，提供了一個(gè)泊松時(shí)間序列預(yù)測的說明性例子。最后，介紹了與擬合非高斯時(shí)間序列建模的其他方法的比較。

緒論

狀態(tài)空間模型為幾種類型的時(shí)間序列和其他數(shù)據(jù)的建模提供了一個(gè)統(tǒng)一的框架。結(jié)構(gòu)性時(shí)間序列、自回歸綜合移動(dòng)平均模型（ARIMA）、簡單回歸、廣義線性混合模型和三次樣條平滑模型只是一些可以表示為狀態(tài)空間模型的統(tǒng)計(jì)模型的例子。最簡單的一類狀態(tài)空間模型是線性高斯?fàn)顟B(tài)空間模型（也被稱為動(dòng)態(tài)線性模型），經(jīng)常被用于許多科學(xué)領(lǐng)域。

高斯?fàn)顟B(tài)空間模型

本節(jié)將介紹有關(guān)高斯?fàn)顟B(tài)空間模型理論的關(guān)鍵概念。由于卡爾曼濾波（Kalman filtering）背后的算法主要是基于Durbin和Koopman（2012）以及同一作者的相關(guān)文章。對于具有連續(xù)狀態(tài)和離散時(shí)間間隔的線性高斯?fàn)顟B(tài)空間模型t=1, . . . ，n，我們有

其中t～N（0，Ht），ηt～N（0，Qt）和α1～N（a1，P1）相互獨(dú)立。我們假設(shè)yt是一個(gè)p×1，αt+1是一個(gè)m×1，ηt是一個(gè)k×1的向量。α = (α > 1 , . . . , α> n ) >，同樣y = (y > 1 , . . , y> n ) >。

狀態(tài)空間建模的主要目標(biāo)是在給定觀測值y的情況下獲得潛狀態(tài)α的知識。這可以通過兩個(gè)遞歸算法實(shí)現(xiàn)，即卡爾曼濾波和平滑算法。從卡爾曼濾波算法中，我們可以得到先行一步的預(yù)測結(jié)果和預(yù)測誤差

和相關(guān)的協(xié)方差矩陣

利用卡爾曼濾波的結(jié)果，我們建立了狀態(tài)平滑方程，在時(shí)間上向后運(yùn)行，產(chǎn)生了

對于干擾項(xiàng)t和ηt，對于信號θt = Ztαt，也可以計(jì)算類似的平滑估計(jì)。

高斯?fàn)顟B(tài)空間模型的例子

現(xiàn)在通過例子來說明。我們的時(shí)間序列包括1969-2007年40-49歲年齡組每年每10萬人中酒精相關(guān)的死亡人數(shù)（圖1）。數(shù)據(jù)取自統(tǒng)計(jì)局。對于觀測值 y1, ... . , yn，我們假設(shè)在所有t = 1, . . . , n，其中μt是一個(gè)隨機(jī)游走的漂移過程

ηt～N(0, σ2 η)。假設(shè)我們沒有關(guān)于初始狀態(tài)μ1或斜率ν的先驗(yàn)信息。這個(gè)模型可以用狀態(tài)空間的形式來寫，定義為

在KFAS中，這個(gè)模型可以用以下代碼來寫。為了說明問題，我們手動(dòng)定義所有的系統(tǒng)矩陣，而不采用默認(rèn)值。

R>?Zt?<-?matrix(c(1,?0),?1,?2)R>?model_gaussian?<-Model(deaths?/?population?~?-1?+custom(Z?=?Zt)

第一個(gè)參數(shù)是定義觀測值的公式（左側(cè)~）和狀態(tài)方程的結(jié)構(gòu)（右側(cè)）。這里死亡人數(shù)/人口是一個(gè)單變量時(shí)間序列，狀態(tài)方程是用矩陣來定義的，為了保持模型的可識別性，截距項(xiàng)用-1省略。觀測水平方差通過參數(shù)H定義，NA值代表未知方差參數(shù)σ 2和σ 2 η。估計(jì)之后，進(jìn)行過濾和平滑遞歸。

KF(fit_gaussian)

在這種情況下，最大似然估計(jì)值對于σ 2是9.5，對于σ 2 η是4.3。從卡爾曼濾波算法中，我們得到了對狀態(tài)的一步超前預(yù)測，at = (μt , νt) 。請注意，即使斜率項(xiàng)ν在我們的模型中被定義為時(shí)間不變量（νt = ν），它也是由卡爾曼濾波算法遞歸估計(jì)的。因此，在每個(gè)時(shí)間點(diǎn)t，當(dāng)新的觀測值yt可用時(shí)，ν的估計(jì)值被更新，以考慮到y(tǒng)t所提供的新信息。在卡爾曼濾波結(jié)束時(shí)，an+1給出了我們對所有數(shù)據(jù)下恒定斜率項(xiàng)的最終估計(jì)。這里斜率項(xiàng)被估計(jì)為0.84，標(biāo)準(zhǔn)誤差為0.34。對于μt，卡爾曼濾波給出了一步的預(yù)測，但是由于狀態(tài)是時(shí)變的，如果我們對t=1, ..., n的μt估計(jì)值感興趣，我們還需要運(yùn)行平滑算法。n的估計(jì)值。

圖1顯示了帶有一步預(yù)測（紅色）和平滑化（藍(lán)色）的隨機(jī)行走過程μt的估計(jì)值的觀察結(jié)果。請注意典型的模型；在時(shí)間t，卡爾曼濾波器計(jì)算一步向前預(yù)測誤差vt = yt - μt，并使用它和先前的預(yù)測來修正下一個(gè)時(shí)間點(diǎn)的預(yù)測。在這里，這在系列的開始階段最容易看到，我們的預(yù)測似乎落后于觀測值一個(gè)時(shí)間步長。另一方面，平滑算法同時(shí)考慮了每個(gè)時(shí)間點(diǎn)的過去和未來的數(shù)值，從而產(chǎn)生了更平滑的潛過程的估計(jì)。

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

卡爾曼濾波器：用R語言中的KFAS建模時(shí)間序列

左右滑動(dòng)查看更多

01

02

03

04

非高斯?fàn)顟B(tài)空間模型的例子

與酒精有關(guān)的死亡也可以自然地被建模為泊松過程?，F(xiàn)在我們的觀測值yt是第t年與酒精有關(guān)的死亡的實(shí)際計(jì)數(shù)，而變化的人口規(guī)模則由暴露項(xiàng)ut來考慮。狀態(tài)方程保持不變，但觀察方程現(xiàn)在的形式是p(yt |μt) = Poisson(ute μt)。

R>?Model(deaths?~?-1?++?distribution?=?"poisson")

與之前的高斯模型相比，我們現(xiàn)在需要用參數(shù)distribution（默認(rèn)為 "高斯"）來定義觀測數(shù)據(jù)的分布。我們還通過參數(shù)u來定義暴露項(xiàng)，并使用a1和P1的默認(rèn)值。在這個(gè)模型中，只有一個(gè)未知參數(shù)，即σ 2 η。這個(gè)參數(shù)被估計(jì)為0.0053，但是高斯模型和泊松模型之間σ 2 η的實(shí)際值不能直接比較，因?yàn)椴煌Ｐ蛯Ζ蘴的解釋不同。泊松模型的斜率項(xiàng)估計(jì)為0.022，標(biāo)準(zhǔn)誤差為1.4×10-4，對應(yīng)于死亡人數(shù)每年增加2.3%。

圖2顯示了以高斯過程（藍(lán)色）和泊松過程（紅色）為模型（每10萬人的死亡人數(shù)）的平滑估計(jì)。

任意的狀態(tài)空間模型

通過結(jié)合前面的方法，可以相對容易地構(gòu)建大量的模型。對于這樣做還不夠的情況，可以通過直接定義系統(tǒng)矩陣來構(gòu)建任意狀態(tài)空間模型。作為一個(gè)例子，我們修改了之前的泊松模型，增加了一個(gè)額外的白噪聲項(xiàng)，試圖捕捉數(shù)據(jù)的可能的過度離散。現(xiàn)在我們的泊松強(qiáng)度模型是ut exp(μt + t)，即

其中ηt～N(0, σ2 η)如前，t～N(0, σ2)。這個(gè)模型可以用狀態(tài)空間的形式來寫，定義為

?Model(deaths?~?trend(2,?Q?=?list(NA,?0))?+distribution?=?"poisson")

由于模型包含P1中的未知參數(shù)，我們需要提供一個(gè)特定的模型更新函數(shù)。

R>?update?<-?function(pars,?model)?{+?model[?"custom"]?<-?exp(pars)+?}

?fit(model_poisson,method?=?"BFGS")

從圖3中我們看到，高斯結(jié)構(gòu)時(shí)間序列模型和帶有額外白噪聲的泊松結(jié)構(gòu)時(shí)間序列模型對平滑趨勢μt的估計(jì)幾乎沒有區(qū)別。這是由于泊松過程的強(qiáng)度相對較高。

例子

我現(xiàn)在用一個(gè)比前面的例子更完整的例子來說明KFAS的使用。數(shù)據(jù)還是由酒精有關(guān)的死亡組成，但現(xiàn)在有四個(gè)年齡組，即30-39歲、40-49歲、50-59歲和60-69歲，被一起作為一個(gè)多變量泊松模型來建模。

1969-2012年的死亡人數(shù)和相應(yīng)年齡組的年人口規(guī)模都有，但作為說明，我們只使用2007年之前的數(shù)據(jù)，并對2008-2013年進(jìn)行預(yù)測。圖4顯示了所有年齡組的每10萬人的死亡人數(shù)。

ts.plot(window(data[,?1:4]?/?data[,?5:8],?end?=?2007)

這里我選擇了之前使用的泊松模型的一個(gè)多變量擴(kuò)展。

這里μt是帶有漂移成分的隨機(jī)游走，νt是一個(gè)恒定的斜率，t是一個(gè)額外的白噪聲成分，用于捕捉序列的額外變化。我對水平和噪聲成分的協(xié)方差結(jié)構(gòu)不做限制。模型（4）可以用KFAS構(gòu)建如下。

R>?Model(Pred[,?1:4]?~+?trend(2,?Q?=?list(matrix(NA,?4,?4))?distribution?=?"poisson"

更新函數(shù)為

R>?updatefn?<-?function(pars,?model,?...){+?model[?etas?]?<-?crossprod(Q)+?crossprod(Q)+?model+?}

我們可以先不通過模擬來估計(jì)模型參數(shù)，然后用這些估計(jì)值作為初始值再次運(yùn)行重要性抽樣的估計(jì)程序。在這種情況下，從重要性抽樣步驟得到的結(jié)果實(shí)際上與從初始步驟得到的結(jié)果相同。

fit(model,?update,+?method?=?"BFGS")R>?fit?<-?fit(model,?updatefn?=?updatefn,?inits?=optimpar)

使用擬合模型的提取方法，我們可以檢查估計(jì)的協(xié)方差和相關(guān)矩陣。

R>?varcordel["Q",???"level"]R>?varcordel["Q",??"custom"]

狀態(tài)空間模型的參數(shù)估計(jì)通常工作量很大，因?yàn)樗迫幻姘鄠€(gè)最大值，從而使優(yōu)化問題高度依賴于初始值。通常情況下，未知參數(shù)與未觀察到的潛在狀態(tài)有關(guān)，如本例中的協(xié)方差矩陣，幾乎沒有先驗(yàn)知識。

因此，要猜出好的初始值是很有挑戰(zhàn)性的，特別是在更復(fù)雜的環(huán)境中。因此，在可以合理地確定找到適當(dāng)?shù)淖顑?yōu)值之前，建議使用多種初始值配置，可能有幾種不同類型的優(yōu)化方法。這里我們使用觀察到的系列的協(xié)方差矩陣作為協(xié)方差結(jié)構(gòu)的初始值。在非高斯模型的情況下，另一個(gè)問題是，似然計(jì)算是基于迭代程序的，該程序使用一些終止條件（如對數(shù)似然的相對變化）停止，因此對數(shù)似然函數(shù)實(shí)際上包含一些噪聲。這反過來又會影響B(tài)FGS等方法的梯度計(jì)算，在理論上可以得到不可靠的結(jié)果。因此，有時(shí)建議使用無導(dǎo)數(shù)的方法，如Nelder-Mead。另一方面，BFGS通常比Nelder-Mead快得多，因此我更愿意先嘗試BFGS，至少在初步分析中。我們可以計(jì)算出狀態(tài)的平滑估計(jì)。

R>?out?<-?KF(model,)

我們看到殘差之間偶爾有滯后的交叉相關(guān)，但總體上我們可以對我們的模型相對滿意?，F(xiàn)在我們可以用我們估計(jì)的模型預(yù)測2008-2013年每個(gè)年齡組與酒精有關(guān)的死亡強(qiáng)度e θt。由于我們的模型是時(shí)間變化的（u變化），我們需要通過newdata參數(shù)為未來的觀察樣本提供模型。

predict(model,+?newdata?++?interval?=?"confidence",)

for?(i?in?1:4)?ts.plot(data[,?i]/data[,?4?+?i],?trend[,?i],?pred[[i]]

圖7顯示了觀察到的死亡人數(shù)，1969-2007年的平滑趨勢，以及2008-2013年的預(yù)測，還有95%的預(yù)測區(qū)間。當(dāng)我們將我們的預(yù)測與真實(shí)的觀察結(jié)果進(jìn)行比較時(shí)，我們看到在現(xiàn)實(shí)中，最年長的年齡組（60-69歲）的死亡人數(shù)略有增加，而在預(yù)測期間，另一個(gè)年齡組的死亡人數(shù)大幅下降。部分原因是在此期間酒精消費(fèi)總量幾乎單調(diào)下降，而這又可能是由于2008年、2009年和2012年酒精稅的增加造成的。

討論

狀態(tài)空間模型提供了解決一大類統(tǒng)計(jì)問題的工具。在這里，我介紹了一個(gè)用于線性狀態(tài)空間建模的方法。

點(diǎn)擊文末?“閱讀原文”

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

本文選自《R語言狀態(tài)空間模型和卡爾曼濾波預(yù)測酒精死亡人數(shù)時(shí)間序列》。

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

matlab實(shí)現(xiàn)擴(kuò)展卡爾曼濾波(EKF)進(jìn)行故障檢測
卡爾曼濾波器：用R語言中的KFAS建模時(shí)間序列
狀態(tài)空間模型：卡爾曼濾波器KFAS建模時(shí)間序列
R語言用LOESS(局部加權(quán)回歸)季節(jié)趨勢分解（STL）進(jìn)行時(shí)間序列異常檢測
使用R語言隨機(jī)波動(dòng)模型SV處理時(shí)間序列中的隨機(jī)波動(dòng)率
PYTHON用時(shí)變馬爾可夫區(qū)制轉(zhuǎn)換（MRS）自回歸模型分析經(jīng)濟(jì)時(shí)間序列
R語言有限混合模型(FMM,finite mixture model)EM算法聚類分析間歇泉噴發(fā)時(shí)間
長短期記憶網(wǎng)絡(luò)LSTM在時(shí)間序列預(yù)測和文本分類中的應(yīng)用
Python隨機(jī)波動(dòng)率(SV)模型對標(biāo)普500指數(shù)時(shí)間序列波動(dòng)性預(yù)測
R語言中ARMA，ARIMA（Box-Jenkins），SARIMA和ARIMAX模型用于預(yù)測時(shí)間序列數(shù)據(jù)
R語言使用ARIMAX預(yù)測失業(yè)率經(jīng)濟(jì)時(shí)間序列數(shù)據(jù)
R語言用ARIMA模型，ARIMAX模型預(yù)測冰淇淋消費(fèi)時(shí)間序列數(shù)據(jù)
R語言經(jīng)濟(jì)學(xué)：動(dòng)態(tài)模型平均(DMA)、動(dòng)態(tài)模型選擇(DMS)預(yù)測原油時(shí)間序列價(jià)格