時間序列模型簡介

平穩(wěn)序列
判斷樣本的平穩(wěn)性
時間序列模型
參數(shù)選擇
實驗代碼

1. 平穩(wěn)序列

時間序列是一列觀測值 $X_t$ 的集合, 其中每個觀測值是在時段 $t$ 觀測所得( $t$ 是自然數(shù) ). 給定時間序列 $\{X_t\}_{t=1}^n$ , 如果對任意的 $t=1,\ldots,n$ , 它滿足下列條件:
i. $\text{var}(X_t)<\infty$
ii. $\mathbb{E}(X_t) = \mu$
iii. $\text{cov}(X_r, X_s) = \text{cov}(X_{r+t}, X_{s+t}), \forall r,s = 1, ..., n$
我們把它叫做(弱)平穩(wěn)(weakly stationary)序列.(下文我們簡稱平穩(wěn)序列.)

通俗地講, 平穩(wěn)序列的期望, 方差, 協(xié)方差不隨時間變化. 例如, $X_t$ 服從同一個分布時, 它是平穩(wěn)的.

例1 下圖中的時間序列由 $X_t\sim N(2, 1)$ 生成. 從直觀上看, 這個序列是"平穩(wěn)的".

Fig1. 平穩(wěn)的時間序列

例2 下圖的中的時間序列由 $X_t= 2 + 0.9X_{t-1}+W_t$ 生成, 其中 $X_0=0$ , $W_t= N(0, 1)$ . 它起初有明顯地增長, 然后趨于平穩(wěn). 利用ADF檢驗(詳情見下文), 我們發(fā)現(xiàn)該序列是平穩(wěn)的(p-value < 0.01).

Fig2. 平穩(wěn)的時間序列

Remark 弱平穩(wěn)性的"弱"主要體現(xiàn)在時間序列在全局上是平穩(wěn)的, 即岸裙，時間序列局部是波動的，但整體上看是平穩(wěn)的, 或者隨著時間的變化其樣本的均值收斂.

2. 判斷樣本的平穩(wěn)性

我們用統(tǒng)計學(xué)中假設(shè)檢驗的方法來判斷樣本的平穩(wěn)性. 常用的是Augmented Dickey-Fuller(ADF)檢驗^[1].

Null Hypothesis $(H_0)$ : 樣本中存在unit root^[2]. 如果接受 $H_0$ , 則意味著樣本是非平穩(wěn)的.
Alternate Hypothesis $(H_1)$ : 樣本中不存在unit root. 如果拒絕 $H_0$ , 則意味著樣本是平穩(wěn)的.

在顯著水平 $\alpha=0.05$ 的條件下, 我們可以通過計算p-value來接受或者拒絕 $H_0$ :

$\text{p-value} > 0.05$ : 接受 $H_0$ .
$\text{p-value} \leq 0.05$ : 拒絕 $H_0$ .

Python3中statsmodels.tas.stattools中的adfuller函數(shù)^[3]實現(xiàn)了ADF檢驗. 使用方法如下所示.

from statsmodels.tsa.stattools import adfuller


def test_stationarity(data, alpha=0.05, print_detail=True):
    """ Test stationarity of time series data.
    :param data: time series data, formatted as list
    :param alpha: significance level.
    :param print_detail: if True print additional information. 
    """
    result = adfuller(data)
    is_stationary = True if result[1] <= alpha else False
    if print_detail:
        print('ADF statistic: %f' % result[0])
        print('p-value: %f' % result[1])
        print('critical values:')
        for key, value in result[4].items():
            print('\t%s: %.3f' % (key, value))
    return is_stationary

3. 時間序列模型

前面之所以介紹平穩(wěn)序列的概念及檢驗方法, 是因為它是很多基礎(chǔ)的時間序列模型的前提假設(shè). 在本節(jié)我們介紹一些常見的時間序列模型(更多內(nèi)容可以參考^[4], ^[5]).

AR(p)

AR代表自回歸(Autoregression). 假設(shè)時間序列 $\{X_t\}$ 是平穩(wěn)的, 它可以被表示成如下形式:
$X_t = \delta + \sum_{i=1}^p \phi_iX_{t-i} + W_t.$

$\delta$ 是常數(shù).
$W_t\overset{\text{iid}}{\sim} N(0, \sigma^2)$ , 它表示時段 $t$ 的誤差(隨機變量).
$p$ 代表自回歸階數(shù).

MA(q)

MA代表移動平均(Moving Average). 假設(shè)時間序列 $\{X_t\}$ 是平穩(wěn)的, 它可以被表示成如下形式:
$X_t = \delta + \sum_{i=1}^q \theta_iW_{t-i} + W_t.$

$\delta$ 是常數(shù).
$W_t\overset{\text{iid}}{\sim} N(0, \sigma^2)$ , 它表示時段 $t$ 的誤差(隨機變量).
$q$ 代表移動平均階數(shù).

ARMA(p,q)

ARMA模型是AR和MA的組合. 假設(shè)同上. 它可以被表示為如下形式:
$X_t = \delta + \sum_{i=1}^p \phi_iX_{t-i} + \sum_{i=1}^q \theta_iW_{t-i} + W_t.$

$p$ 是自回歸階數(shù)
$q$ 是移動平均階數(shù)

ARIMA(p,d,q)

ARIMA模型是ARMA模型的推廣, 全稱是Autoregressive Integrated Moving Average. 當(dāng)時間序列 $\{X_t\}$ 不滿足平穩(wěn)性時, 我們通常使用差分的技巧把序列變得平穩(wěn), 然后再應(yīng)用ARMA模型.

參數(shù) $d$ 代表差分的階數(shù). 下面是差分的計算公式( $\Delta$ 為差分算子):

一階差分 $\Delta X_t = X_t - X_{t-1}$
二階差分 $\Delta^{(2)} X_t = \Delta(X_t-X_{t-1}) = X_t - 2X_{t-1}+X_{t-2} = (1-\Delta)^2X_t$
$d$ 階差分 $\Delta^{(d)} X_t = (1-\Delta)^dX_t$

例3 下圖是原始的時間序列. 通過觀察, 它的均值有明顯的上升趨勢且不收斂, 因此不是平穩(wěn)序列(ADF檢驗的p-value為0.94).

對該序列進行一階差分后, 我們得到如下平穩(wěn)的時間序列(p-value為0.00).

ARIMA(p,d,q) $\times$ (P,D,Q,s)

該記號代表季節(jié)性(或周期性)ARIMA模型, 詳細的表達式可以參考^[4](4.1 Seasonal ARIMA models), 其中

$p,d,q$ 的意義同上.
$P$ 代表周期性自回歸階數(shù)(前 $P$ 個周期對應(yīng)觀測值的自回歸).
$D$ 代表周期性差分階數(shù).
$Q$ 代表周期性移動平均階數(shù)(前 $Q$ 個周期對應(yīng)的移動平均).
$s$ 代表一個周期的長度.

我們可以把它看成兩階段模型: 第一階段在全局使用ARIMA(p,d,q); 第二階段通過指定周期長度 $s$ , 再利用ARIMA(P,Q,D)模型考慮周期之間的關(guān)系.

例4 考慮如下周期性的平穩(wěn)時間序列( $s=18$ ).

season.png

對序列進行周期性差分: $Y_t = X_t - X_{t-18}$ 得到新的時間序列 $\{Y_t\}$ 如下圖所示(紅色部分)

deseason.png

通過使用周期性差分, 我們可以把原有時間序列的周期性移除. 同理, 通過采用周期性的自回歸和移動平均系數(shù), 我們可以把周期之間的依賴關(guān)系考慮進模型.

例5 考慮周期s=18的數(shù)據(jù)(藍色曲線). 用 $\text{ARIMA}(1,0, 0)$ 和 $\text{ARIMA}(1, 0, 0)\times(0, 1, 0, 18)$ 分別進行預(yù)測的結(jié)果如下.

不考慮周期性的ARIMA模型的預(yù)測結(jié)果(灰色曲線)逐漸收斂到時間序列的均值. 由于序列是平穩(wěn)的, 這樣的預(yù)測結(jié)果符合我們的期望. 考慮到該時間序列有比較強的周期性, 且通過觀察發(fā)現(xiàn)周期 $s=18$ . 在本例中, 我們僅使用周期差分, 最終得到了如圖所示(紅色曲線)的周期性預(yù)測結(jié)果.

ARCH(p)

ARCH的全稱是Autoregressive Conditionally Heteroscedasticity, 它可以用來考慮樣本的方差隨著時間變化(或震蕩)的時間序列. 設(shè)時間序列 $\{X_t\}$ 是平穩(wěn)的, $\text{ARCH}(p)$ 模型可以被表示成如下形式:

$X_t = \sigma_tW_t,\quad W_t \overset{\text{iid}}{\sim} N(0,1)$
其中
$\sigma_t^2 = \alpha_0 + \sum_{i=1}^p \alpha_iX_{t-i}^2.$

$p$ 代表 $X_t^2$ 的自回歸階數(shù).

GARCH(p,q)

GARCH即Generalized ARCH, 是ARCH模型的推廣^[6]. 設(shè)時間序列 $\{X_t\}$ 是平穩(wěn)的, $\text{GARCH}(p,q)$ 模型可以被表示成如下形式:
$X_t = \sigma_tW_t,\quad W_t \overset{\text{iid}}{\sim} N(0,1)$
其中
$\sigma_t^2 = \alpha_0 + \sum_{i=1}^p \alpha_iX_{t-i}^2 + \sum_{i=1}^q \beta_i\sigma^2_{t-i}.$

$p$ 代表 $X_t^2$ 的自回歸階數(shù).
$q$ 代表 $\sigma_t^2$ 的移動平均階數(shù).

Remark ARCH/GARCH隨機過程產(chǎn)生的數(shù)據(jù)是什么樣的? 前面提到它們允許樣本的方差隨時間變化, 但是由于 $\{X_t\}$ 必須滿足平穩(wěn)性(前提假設(shè)), 因此樣本的方差從局部看是變化(震蕩)的, 但從整體看應(yīng)該是"平穩(wěn)的"序列. 例如下圖是一個 $\text{GARCH}(1,1)$ 過程生成的時間序列( $\alpha_0=5, \alpha_1=\beta_1=0.5$ ).

VAR(p)

VAR即Vector Autoregression, 它是多變量的自回歸模型. 類似地, 我們有 $\text{VARMA}(p, q)$ , 它是 $\text{ARMA}(p,q)$ 的向量版本. 需要注意的是, VARMA模型處理的時間序列可以有趨勢. 我們不做詳細的展開, 感興趣的讀者可以參考^[4]章節(jié)11.2: Vector Autoregressive models VAR(p) models.

4. 參數(shù)選擇

給定時間序列的觀測樣本, 選定預(yù)測模型之后如何確定模型的參數(shù)? 本節(jié)我們介紹兩種常用的方法: 1. 畫出ACF/PACF圖, 然后觀察出 $p,q$ 的值; 2. 通過計算相關(guān)的統(tǒng)計指標, 自動化地選擇參數(shù).

4.1 觀察ACF/PACF

ACF

ACF的全稱是Autocorrelation Function. 對變量 $h=1,2, ...$ , ACF的值代表 $X_t$ 與 $X_{t-h}$ 之間的相關(guān)性.
$\text{ACF}(h) = \frac{\text{cov}(X_t, X_{t-h})}{\text{var}(X_t)}$

PACF

PACF的全稱是Partial Autocorrelation Function. 對變量 $h=1,2,...$ , PACF的值代表已知 $X_{t-1}, ... X_{t-h+1}$ 的條件下, $X_t$ 與 $X_{t-h}$ 之間的相關(guān)性.
$\text{PACF}(h) = \frac{\text{cov}(X_t,X_{t-h}|x_{t-1}, ... x_{t-h+1})}{\sqrt{\text{var}(X_t|x_{t-1}, ... x_{t-h+1})\text{var}(X_{t-h}|x_{t-1}, ... x_{t-h+1})}}$

例6 設(shè) $W_t\sim N(0, 1)$ . 考慮下面三個模型生成的時間序列, 并計算相應(yīng)的ACF/PACF.

$\text{AR}(1): X_t=2 + 0.5X_{t-1}+W_t$
$\text{MA}(1): X_t=2 + 0.5W_{t-1}+W_t$
$\text{ARMA}(1,1): X_t=2 + 0.5X_{t-1}+ 0.5W_{t-1}+W_t$

指導(dǎo)原則(參考^[4])

模型	ACF	PACF	說明
AR(p)	逐漸趨近于0或像正弦曲線一樣收斂到0	前 $p$ 個值非常顯著, 其余的值不顯著	$p$ 值主要參考PACF
MA(q)	前 $q$ 個值非常顯著, 其余的值不顯著	逐漸趨近于0或像正弦曲線一樣收斂到0	$q$ 值主要參考ACF
ARMA(p,q)	逐漸趨近于0或像正弦曲線一樣收斂到0	逐漸趨近于0或像正弦曲線一樣收斂到0	$p,q$ 值靠猜

4.2 自動化地決定參數(shù)

基本思想是通過計算一些指標, 并選擇參數(shù)使得相關(guān)的指標值盡可能小. 下面我們介紹一些常用的指標.

為方便描述, 我們先定義一些記號.

$n$ = 樣本的大小
$k$ = 模型中需要擬合的參數(shù)數(shù)量(例如正態(tài)分布有的數(shù)量是2: $\mu$ 和 $\sigma$ )
$L_{\max}$ = 通過最大似然估計得到的最大Likelihood

AIC(Akaike Information Criterion)^[7]

$\text{AIC} = 2k - 2\ln(L_{\max})$

AICc^[8]

(AIC的改良版, 解決小樣本過擬合的問題)
$\text{AICc} = \text{AIC} + \frac{2k^2 + 2k}{n-k-1}$

BIC(Bayesian Information Criterion)^[9]

(也稱為Schwartz Criterion, SBC, SBIC)

$\text{BIC} = \ln(n)k - 2\ln(L_{\max})$

HQIC(Hannan–Quinn Information Criterion)^[10]

$\text{HQIC} = 2k\ln(\ln(n)) -2\ln(L_{\max})$

Remark 建議在實際中綜合考慮這些指標.

5. 實驗代碼

Python3 code on Github

參考文獻

Wikipedia. Augmented Dickey-Fuller test. ?
Dickey, D. A.; Fuller, W. A. Distribution of the Estimators for Autoregressive Time Series with a Unit Root. Journal of the American Statistical Association. 74(366): 427–431, 1979. ?
https://www.statsmodels.org/dev/generated/statsmodels.tsa.stattools.adfuller.html ?
Lecture notes of Applied Time Series Analysis (SATA 510). The Pennsylvania State University. ? ? ? ?
Jan Grandell. Time series analysis (lecture notes). ?
Bollerslev, T. Generalized autoregressive conditional heteroscedasticity. Journal of Econometrics, 31, 307–327, 1986. ?
Hirotugu Akaike. A new look at the statistical model identification, IEEE Transactions on Automatic Control, 19 (6): 716–723, 1974. ?
Wikipedia. https://en.wikipedia.org/wiki/Akaike_information_criterion#AICc. ?
Schwartz E.S. The Stochastic Behavior of Commodity Prices: Implications for Valuation and Hedging'. J Finance 52(3) Papers and Proceedings Fifty-Seventh Annual Meeting, American Finance Association, New Orleans, Louisiana, 923-973, 1997. ?
Hannan, E. J. and B. G. Quinn. The Determination of the order of an autoregression, Journal of the Royal Statistical Society, Series B, 41: 190–195, 1979. ?

最后編輯于：2018.10.20 10:46:09

?著作權(quán)歸作者所有,轉(zhuǎn)載或內(nèi)容合作請聯(lián)系作者

人面猴
序言：七十年代末，一起剝皮案震驚了整個濱河市，隨后出現(xiàn)的幾起案子映砖，更是在濱河造成了極大的恐慌，老刑警劉巖弯淘，帶你破解...
沈念sama閱讀 217,185評論 6贊 503
死咒
序言：濱河連續(xù)發(fā)生了三起死亡事件颜武，死亡現(xiàn)場離奇詭異，居然都是意外死亡拙寡，警方通過查閱死者的電腦和手機授滓，發(fā)現(xiàn)死者居然都...
沈念sama閱讀 92,652評論 3贊 393
救了他兩次的神仙讓他今天三更去死
文/潘曉璐我一進店門，熙熙樓的掌柜王于貴愁眉苦臉地迎上來肆糕，“玉大人般堆，你說我怎么就攤上這事〕峡校” “怎么了郁妈？”我有些...
開封第一講書人閱讀 163,524評論 0贊 353
道士緝兇錄：失蹤的賣姜人
文/不壞的土叔我叫張陵，是天一觀的道長绍申。經(jīng)常有香客問我噩咪，道長顾彰，這世上最難降的妖魔是什么？我笑而不...
開封第一講書人閱讀 58,339評論 1贊 293
?港島之戀（遺憾婚禮）
正文為了忘掉前任胃碾，我火速辦了婚禮涨享，結(jié)果婚禮上，老公的妹妹穿的比我還像新娘仆百。我一直安慰自己厕隧，他們只是感情好，可當(dāng)我...
茶點故事閱讀 67,387評論 6贊 391
惡毒庶女頂嫁案：這布局不是一般人想出來的
文/花漫我一把揭開白布俄周。她就那樣靜靜地躺著吁讨，像睡著了一般。火紅的嫁衣襯著肌膚如雪峦朗。梳的紋絲不亂的頭發(fā)上建丧，一...
開封第一講書人閱讀 51,287評論 1贊 301
城市分裂傳說
那天，我揣著相機與錄音波势，去河邊找鬼翎朱。笑死，一個胖子當(dāng)著我的面吹牛尺铣，可吹牛的內(nèi)容都是我干的拴曲。我是一名探鬼主播，決...
沈念sama閱讀 40,130評論 3贊 418
雙鴛鴦連環(huán)套：你想象不到人心有多黑
文/蒼蘭香墨我猛地睜開眼凛忿，長吁一口氣：“原來是場噩夢啊……” “哼澈灼！你這毒婦竟也來了？” 一聲冷哼從身側(cè)響起店溢，我...
開封第一講書人閱讀 38,985評論 0贊 275
萬榮殺人案實錄
序言：老撾萬榮一對情侶失蹤叁熔，失蹤者是張志新（化名）和其女友劉穎，沒想到半個月后逞怨，有當(dāng)?shù)厝嗽跇淞掷锇l(fā)現(xiàn)了一具尸體者疤，經(jīng)...
沈念sama閱讀 45,420評論 1贊 313
?護林員之死
正文獨居荒郊野嶺守林人離奇死亡，尸身上長有42處帶血的膿包…… 初始之章·張勛以下內(nèi)容為張勛視角年9月15日...
茶點故事閱讀 37,617評論 3贊 334
?白月光啟示錄
正文我和宋清朗相戀三年叠赦，在試婚紗的時候發(fā)現(xiàn)自己被綠了驹马。大學(xué)時的朋友給我發(fā)了我未婚夫和他白月光在一起吃飯的照片。...
茶點故事閱讀 39,779評論 1贊 348
活死人
序言：一個原本活蹦亂跳的男人離奇死亡除秀，死狀恐怖糯累，靈堂內(nèi)的尸體忽然破棺而出，到底是詐尸還是另有隱情册踩，我是刑警寧澤泳姐，帶...
沈念sama閱讀 35,477評論 5贊 345
?日本核電站爆炸內(nèi)幕
正文年R本政府宣布，位于F島的核電站暂吉，受9級特大地震影響胖秒，放射性物質(zhì)發(fā)生泄漏缎患。R本人自食惡果不足惜，卻給世界環(huán)境...
茶點故事閱讀 41,088評論 3贊 328
男人毒藥：我在死后第九天來索命
文/蒙蒙一阎肝、第九天我趴在偏房一處隱蔽的房頂上張望挤渔。院中可真熱鬧，春花似錦风题、人聲如沸判导。這莊子的主人今日做“春日...
開封第一講書人閱讀 31,716評論 0贊 22
一樁弒父案沛硅，背后竟有這般陰謀
文/蒼蘭香墨我抬頭看了看天上的太陽眼刃。三九已至，卻和暖如春摇肌，著一層夾襖步出監(jiān)牢的瞬間擂红，已是汗流浹背。一陣腳步聲響...
開封第一講書人閱讀 32,857評論 1贊 269
情欲美人皮
我被黑心中介騙來泰國打工朦蕴，沒想到剛下飛機就差點兒被人妖公主榨干…… 1. 我叫王不留篮条，地道東北人弟头。一個月前我還...
沈念sama閱讀 47,876評論 2贊 370
代替公主和親
正文我出身青樓吩抓，卻偏偏與公主長得像，于是被迫代替她去往敵國和親赴恨。傳聞我的和親對象是個殘疾皇子疹娶，可洞房花燭夜當(dāng)晚...
茶點故事閱讀 44,700評論 2贊 354