時間序列ARMA中p，q選擇

imelee 2018-09-17

展開全文

時間序列中p,q值選擇

1.模型識別：

對平穩(wěn)時間序列Yn，求得其自相關(guān)函數(shù)（ACF）和偏自相關(guān)函數(shù)（PACF）序列。
若PACF序列滿足在p步截尾，且ACF序列被負(fù)指數(shù)函數(shù)控制收斂到0，則Yn為AR(p)序列。
若ACF序列滿足在q步截尾，且PACF序列被負(fù)指數(shù)函數(shù)控制收斂到0，則Yn為MA(q)序列。
若ACF序列和PACF序列滿足皆不截尾，但都被負(fù)指數(shù)函數(shù)控制收斂到0，則Yn為ARMA序列。

2.模型定階：

對于有N個觀察值的序列，求得相應(yīng)于AR(p)、MA(q) 和 ARMA(p,q)三種模型的殘差方差，出現(xiàn)模型最小殘差方差時的模型階數(shù)就是各個模型的最佳階數(shù)。

3.模型參數(shù)估計：

AR模型的參數(shù)可以根據(jù)ACF序列構(gòu)成的矩陣及其矩陣之間的轉(zhuǎn)化關(guān)系求得。
MA模型的參數(shù)采用線性迭代法即可求出。
ARMA模型參數(shù)估計方法是按上述求解AR模型和MA模型參數(shù)的方法分別對AR和MA模型進(jìn)行參數(shù)估計，即可得到ARMA模型的參數(shù)。

4.模型估計函數(shù)：

根據(jù)對應(yīng)的模型以及估計參數(shù)等帶入估計函數(shù)計算出估計值。

python中使用Statsmodel庫進(jìn)行p，q值選擇

最近工作中用到時間序列模型中的ARMA（ARIMA），需要自動選擇p，q值，但是我查到的資料都是根據(jù)自相關(guān)圖和偏自相關(guān)圖來觀察拖尾和截尾，以此來選擇p，q值，剛開始時一籌莫展，后來靈機(jī)一動，為何不看下導(dǎo)入的statsmodel庫中對應(yīng)畫圖的函數(shù)調(diào)用時用的源代碼呢：

import statsmodels

print statsmodels.__file__1
2
3

根據(jù)輸出找到源代碼所在位置，用來畫偏自相關(guān)圖的代碼部分為：

acf_x, confint = pacf(x, nlags=nlags, alpha=alpha, method=method)

    if use_vlines:
        ax.vlines(lags, [0], acf_x, **kwargs)
        ax.axhline(**kwargs)

    # center the confidence interval TODO: do in acf?
    confint = confint - confint.mean(1)[:,None]
    kwargs.setdefault('marker', 'o')
    kwargs.setdefault('markersize', 5)
    kwargs.setdefault('linestyle', 'None')
    ax.margins(.05)
    ax.plot(lags, acf_x, **kwargs)
    ax.fill_between(lags, confint[:,0], confint[:,1], alpha=.25)1
2
3
4
5
6
7
8
9
10
11
12
13
14

從畫圖部分可以看到置信區(qū)間上界為confint[:,0]，下界線為 confint[:,1]，又有：

confint = confint - confint.mean(1)[:,None]1

因此可以寫個循環(huán)，當(dāng)出現(xiàn)截尾時返回當(dāng)前p,q值，也就可以確定選擇AR或MA模型了，如果兩者都拖尾，則需要用赤池信息準(zhǔn)則，或者貝葉斯信息準(zhǔn)則或別的準(zhǔn)則來進(jìn)一步判斷，那就是另一回事了，除了運(yùn)行時間稍微有點長外并沒有什么困難的了。

代碼如下：

1、判斷時間序列穩(wěn)定性

#用來檢查時間序列穩(wěn)定性的，代碼中選用的臨界值為5%，p-value選用的為0.1，這個可以根據(jù)實際進(jìn)行修改。

from statsmodels.tsa.arima_model import ARIMA
from statsmodels.tsa.stattools import adfuller, acf, pacf


def testStationarity(timeSer):

    stationarity = False

    dftest = adfuller(timeSer)
    dfoutput = Series(dftest[:4], index=[
                      'Test Statistic', 'p-value', 'lags', 'nobs'])

    for key, value in dftest[4].items():
        dfoutput['Critical values (%s)' % key] = value

    if dfoutput['Test Statistic'] < dfoutput['Critical values (5%)']:
        if dfoutput['p-value'] < 0.1:
            stationarity = True

    return stationarity
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23

2、選擇合適的p，q值

（由于自動判斷，因此收斂速度狀況并沒有做判斷，這個是因為我使用的數(shù)據(jù)在滿足第一個在置信區(qū)間內(nèi)的值后即是局部最優(yōu)，實際中因數(shù)據(jù)的不同代碼需要部分修改，這里只是提供了一個思路。）

def p_q_choice(timeSer, nlags=40, alpha=.05):

    kwargs = {'nlags': nlags, 'alpha': alpha}
    acf_x, confint = acf(timeSer, **kwargs)
    acf_px, confint2 = pacf(timeSer, **kwargs)

    confint = confint - confint.mean(1)[:, None]
    confint2 = confint2 - confint2.mean(1)[:, None]

    for key1, x, y, z in zip(range(nlags), acf_x, confint[:,0], confint[:,1]):
        if x > y and x < z:
            q = key1
            break

    for key2, x, y, z in zip(range(nlags), acf_px, confint2[:,0], confint[:,1]):
        if x > y and x < z:
            p = key2
            break

    return p, q
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21

附錄：
ARIMA模型運(yùn)用的流程

根據(jù)時間序列的散點圖、自相關(guān)函數(shù)和偏自相關(guān)函數(shù)圖識別其平穩(wěn)性。
對非平穩(wěn)的時間序列數(shù)據(jù)進(jìn)行平穩(wěn)化處理。直到處理后的自相關(guān)函數(shù)和偏自相關(guān)函數(shù)的數(shù)值非顯著非零
根據(jù)所識別出來的特征建立相應(yīng)的時間序列模型。
平穩(wěn)化處理后，若偏自相關(guān)函數(shù)是截尾的，而自相關(guān)函數(shù)是拖尾的，則建立AR模型；
若偏自相關(guān)函數(shù)和自相關(guān)函數(shù)均是拖尾的，則序列適合ARMA模型。
參數(shù)估計，檢驗是否具有統(tǒng)計意義。
假設(shè)檢驗，判斷（診斷）殘差序列是否為白噪聲序列。
利用已通過檢驗的模型進(jìn)行預(yù)測。