時(shí)間序列中p,q值選擇

1.模型識(shí)別：

對(duì)平穩(wěn)時(shí)間序列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.模型定階：

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

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

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

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

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

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

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

import statsmodels

print statsmodels.__file__
1
2
3

根據(jù)輸出找到源代碼所在位置，用來(lái)畫(huà)偏自相關(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

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

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

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

代碼如下：

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

#用來(lái)檢查時(shí)間序列穩(wěn)定性的，代碼中選用的臨界值為5%，p-value選用的為0.1，這個(gè)可以根據(jù)實(shí)際進(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值

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

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)用的流程

1. 根據(jù)時(shí)間序列的散點(diǎn)圖、自相關(guān)函數(shù)和偏自相關(guān)函數(shù)圖識(shí)別其平穩(wěn)性。

2. 對(duì)非平穩(wěn)的時(shí)間序列數(shù)據(jù)進(jìn)行平穩(wěn)化處理。直到處理后的自相關(guān)函數(shù)和偏自相關(guān)函數(shù)的數(shù)值非顯著非零

3. 根據(jù)所識(shí)別出來(lái)的特征建立相應(yīng)的時(shí)間序列模型。
   平穩(wěn)化處理后，若偏自相關(guān)函數(shù)是截尾的，而自相關(guān)函數(shù)是拖尾的，則建立AR模型；
   若偏自相關(guān)函數(shù)和自相關(guān)函數(shù)均是拖尾的，則序列適合ARMA模型。

4. 參數(shù)估計(jì)，檢驗(yàn)是否具有統(tǒng)計(jì)意義。

5. 假設(shè)檢驗(yàn)，判斷（診斷）殘差序列是否為白噪聲序列。

6. 利用已通過(guò)檢驗(yàn)的模型進(jìn)行預(yù)測(cè)。