一、資料平穩性與差分法

1.1 平穩性

• 平穩性就是要求經由樣本時間序列所得到的擬合曲線在未來的一段期間內仍能順著現有的形態“慣性”地延續下去
• 平穩性要求序列的均值和方差不發生明顯變化

1.2 嚴平穩與弱平穩

嚴平穩：嚴平穩表示的分佈不隨時間的改變而改變。
如：白噪聲（正態），無論怎麼取，都是期望為0，方差為1

弱平穩：期望與相關係數（依賴性）不變
未來某時刻的t的值Xt就要依賴於它的過去資訊，所以需要依賴性

嚴平穩的條件只是理論上的存在，現實中用得比較多的是寬平穩的條件。

弱平穩也叫寬平穩或者二階平穩（均值和方差平穩），它應滿足：

• 常數均值
• 常數方差
• 常數自協方差

ARIMA 模型對時間序列的要求是平穩型。因此，當你得到一個非平穩的時間序列時，首先要做的即是做時間序列的差分，直到得到一個平穩時間序列。如果你對時間序列做d次差分才能得到一個平穩序列，那麼可以使用ARIMA(p,d,q)模型，其中d是差分次數。

1.3 差分法

時間序列在t與t-1時刻的差值

從統計學的角度來看，資料差分是指將非平穩資料轉換為平穩的過程，去除其非恆定的趨勢。“差分消除了時間序列水平的變化，消除了趨勢和季節性，從而穩定了時間序列的平均值。”

1.4 差分實驗

資料來自：https://github.com/SimiY/pydata-sf-2016-arima-tutorial

import sys
import os
 
import pandas as pd
import numpy as np
 
import matplotlib.pylab as plt
import seaborn as sns
 
pd.set_option('display.float_format', lambda x: '%.5f' % x) # pandas
np.set_printoptions(precision=5, suppress=True) # numpy
 
pd.set_option('display.max_columns', 100)
pd.set_option('display.max_rows', 100)
 
# seaborn plotting style
sns.set(style='ticks', context='poster')

#美國消費者信心指數
Sentiment = 'data/sentiment.csv'
Sentiment = pd.read_csv(Sentiment, index_col=0, parse_dates=[0])
Sentiment.head()
 

# Select the series from 2005 - 2016
sentiment_short = Sentiment.loc['2005':'2016']

sentiment_short.plot(figsize=(12,8))
plt.legend(bbox_to_anchor=(1.25, 0.5))
plt.title("Consumer Sentiment")
sns.despine()

sentiment_short['diff_1'] = sentiment_short['UMCSENT'].diff(1)#求差分值，一階差分。     1指的是1個時間間隔，可更改。
 
sentiment_short['diff_2'] = sentiment_short['diff_1'].diff(1)#再求差分，二階差分。
 
sentiment_short.plot(subplots=True, figsize=(18, 12))
 

二、ARIMA模型拆解

基本模型：自迴歸移動平均模型(ARMA(p，q))是時間序列中最為重要的模型之一。它主要由兩部分組成： AR代表p階自迴歸過程，MA代表q階移動平均過程。

2.1 AR - 自迴歸

自迴歸模型的限制:

• 自迴歸模型是用自身的資料來進行預測
• 必須具有平穩性
• 必須具有自相關性，如果自相關係數(φi)小於0.5，則不宜採用
• 自迴歸只適用於預測與自身前期相關的現象

2.2 MA - 移動平均

2.3 I - 表示綜合

ARIMA中的“I”指的是它的綜合方面。當應用差分步驟時，資料是“綜合”的，以消除非平穩性。表示原始觀測值的差異，以允許時間序列變得平穩，即資料值被資料值和以前的值之間的差異替換。

I表示差分項，用d表示，等於1為1階差分，2為2階差分，等於0不用做，一般做1階差分就夠了。

2.4 最終模型

ARIMA(p，d，q)模型全稱為差分自迴歸移動平均模型(Autoregressive Integrated Moving Average Model,簡記ARIMA)

• AR是自迴歸，p為自迴歸項；
• MA為移動平均
• q為移動平均項數，d為時間序列成為平穩時所做的差分次數

原理：將非平穩時間序列轉化為平穩時間序列然後將因變數。僅對它的滯後值以及隨機誤差項的現值和滯後值進行迴歸所建立的模型。

滯後值指的是階數

三、自相關和偏自相關

3.1 自相關函式ACF

• 自相關函式ACF(autocorrelation function)
• 對一個時間序列，現在值與其過去值的相關性。如果相關性為正，則說明現有趨勢將繼續保持。
• 自相關函式反映了同一序列在不同時序的取值之間的相關性

公式：

Pk的取值範圍為[-1,1]

3.2 偏自相關函式(PACF)

• 偏自相關函式(PACF)(partial autocorrelation function)

• 對於一個平穩AR(p)模型，求出滯後k自相關係數p(k)時實際上得到並不是x(t)與x(t-k)之間單純的相關關係

• x(t)同時還會受到中間k-1個隨機變數x(t-1)、x(t-2)、……、x(t-k+1)的影響，而這k-1個隨機變數又都和x(t-k)具有相關關係。所以自相關係數p(k)裡實際摻雜了其他變數對x(t)與x(t-k)的影響

• 剔除了中間k-1個隨機變數x(t-1)、x(t-2)、……、x(t-k+1)的干擾之後x(t-k)對x(t)影響的相關程度。

• ACF還包含了其他變數的影響，而偏自相關係數PACF是嚴格這兩個變數之間的相關性

https://www.cnblogs.com/tianqizhi/p/9277376.html

https://blog.csdn.net/fengdu78/article/details/121347188

3.3 pdq階數確定

截尾：落在置信區間內（95%的點都符合該規則），可簡單理解為從某階之後直接就變為0。

3.4 ACF、PACF求解

del sentiment_short['diff_2']
del sentiment_short['diff_1']
sentiment_short.head()
print (type(sentiment_short))

<class 'pandas.core.frame.DataFrame'>

from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
fig = plt.figure(figsize=(12,8))
#acf 自相關
ax1 = fig.add_subplot(211)
plot_acf(sentiment_short,lags=20,title='acf', ax=ax1)  
ax1.xaxis.set_ticks_position('bottom')
fig.tight_layout();

#pacf 偏自相關
ax2 = fig.add_subplot(212)
plot_pacf(sentiment_short,lags=20,title='pacf', ax=ax2)  
ax2.xaxis.set_ticks_position('bottom')
fig.tight_layout();

# 散點圖也可以表示
 
lags=9
ncols=3
nrows=int(np.ceil(lags/ncols))
 
fig, axes = plt.subplots(ncols=ncols, nrows=nrows, figsize=(4*ncols, 4*nrows))
 
for ax, lag in zip(axes.flat, np.arange(1,lags+1, 1)):
    lag_str = 't-{}'.format(lag)
    X = (pd.concat([sentiment_short, sentiment_short.shift(-lag)], axis=1,
                   keys=['y'] + [lag_str]).dropna())
 
    X.plot(ax=ax, kind='scatter', y='y', x=lag_str);
    corr = X.corr().values[0][1]
    ax.set_ylabel('Original')
    ax.set_title('Lag: {} (corr={:.2f})'.format(lag_str, corr));
    ax.set_aspect('equal');
    sns.despine();
 
fig.tight_layout();

# 更直觀一些
#模板，使用時直接改自己的資料就行，用以下四個圖進行評估和分析就可以
def tsplot(y, lags=None, title='', figsize=(20, 8)):
    
    fig = plt.figure(figsize=figsize)
    layout = (2, 2)
    ts_ax   = plt.subplot2grid(layout, (0, 0))
    hist_ax = plt.subplot2grid(layout, (0, 1))
    acf_ax  = plt.subplot2grid(layout, (1, 0))
    pacf_ax = plt.subplot2grid(layout, (1, 1))
     
    y.plot(ax=ts_ax)
    ts_ax.set_title(title)
    y.plot(ax=hist_ax, kind='hist', bins=25)
    hist_ax.set_title('Histogram')
    plot_acf(y, lags=lags, ax=acf_ax)
    plot_pacf(y, lags=lags, ax=pacf_ax)
    [ax.set_xlim(0) for ax in [acf_ax, pacf_ax]]
    sns.despine()
    plt.tight_layout()
    return ts_ax, acf_ax, pacf_ax
 
tsplot(sentiment_short, title='Consumer Sentiment', lags=36);

四、AIC和BIC

AIC和BIC兩個值均是越低越好

AIC表達含義是引數和最終結果之間的權衡

AIC和BIC就是在保持精度的情況下，使K值越小越好

五、ARIMA完整流程

%load_ext autoreload
%autoreload 2
%matplotlib inline
%config InlineBackend.figure_format='retina'
 
from __future__ import absolute_import, division, print_function
import sys
import os
 
import pandas as pd
import numpy as np
 
import matplotlib.pylab as plt
import seaborn as sns
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.statespace.sarimax import SARIMAX
 
pd.set_option('display.float_format', lambda x: '%.5f' % x) # pandas
np.set_printoptions(precision=5, suppress=True) # numpy
 
pd.set_option('display.max_columns', 100)
pd.set_option('display.max_rows', 100)
 
# seaborn plotting style
sns.set(style='ticks', context='poster')

讀入資料

filename_ts = 'data/series2.csv'
ts_df = pd.read_csv(filename_ts, index_col=0, parse_dates=[0])
n_sample = ts_df.shape[0]
print(ts_df.shape)
print(ts_df.head())

(250, 1)
value
1998-06-01 -0.59883
1998-07-01 -0.80018
1998-08-01 2.29898
1998-09-01 1.15039
1998-10-01 -1.19258

劃分訓練集和測試集

# Create a training sample and testing sample before analyzing the series
 
n_train=int(0.95*n_sample)+1
n_forecast=n_sample-n_train
#ts_df
ts_train = ts_df.iloc[:n_train]['value']
ts_test = ts_df.iloc[n_train:]['value']
print(ts_train.shape)
print(ts_test.shape)
print("Training Series:", "\n", ts_train.tail(), "\n")
print("Testing Series:", "\n", ts_test.head())

(238,)
(12,)
Training Series:
2017-11-01 1.08165
2017-12-01 0.40939
2018-01-01 2.17708
2018-02-01 2.21133
2018-03-01 -0.39728
Name: value, dtype: float64

Testing Series:
2018-04-01 0.82064
2018-05-01 0.66053
2018-06-01 0.78946
2018-07-01 -0.02444
2018-08-01 -0.39888
Name: value, dtype: float64

def tsplot(y, lags=None,figsize=(14, 8)):
     
    fig = plt.figure(figsize=figsize)
    layout = (2, 2)
    ts_ax   = plt.subplot2grid(layout, (0, 0))
    hist_ax = plt.subplot2grid(layout, (0, 1))
    acf_ax  = plt.subplot2grid(layout, (1, 0))
    pacf_ax = plt.subplot2grid(layout, (1, 1))
    y.plot(ax=ts_ax)
    ts_ax.set_title('A Given Training Series')
    y.plot(ax=hist_ax, kind='hist', bins=25)
    hist_ax.set_title('Histogram')
    plot_acf(y, lags=lags, ax=acf_ax)
    plot_pacf(y, lags=lags, ax=pacf_ax)
    [ax.set_xlim(0) for ax in [acf_ax, pacf_ax]]
    sns.despine()
    
    fig.tight_layout()
    
    return ts_ax, acf_ax, pacf_ax
    
tsplot(ts_train, lags=20);

模型訓練

#Model Estimation
 
# Fit the model
arima200 = SARIMAX(ts_train, order=(2,0,0))#order裡邊的三個引數p,d,q
model_results = arima200.fit()#fit模型

import itertools
#當多組值都不符合時，遍歷多組值，得出最好的值
p_min = 0
d_min = 0
q_min = 0
p_max = 4
d_max = 0
q_max = 4
 
# Initialize a DataFrame to store the results
results_bic = pd.DataFrame(index=['AR{}'.format(i) for i in range(p_min,p_max+1)],
                           columns=['MA{}'.format(i) for i in range(q_min,q_max+1)])
 
for p,d,q in itertools.product(range(p_min,p_max+1),
                               range(d_min,d_max+1),
                               range(q_min,q_max+1)):
    if p==0 and d==0 and q==0:
        results_bic.loc['AR{}'.format(p), 'MA{}'.format(q)] = np.nan
        continue
     
    try:
        model = SARIMAX(ts_train, order=(p, d, q),
                               #enforce_stationarity=False,
                               #enforce_invertibility=False,
                              )
        results = model.fit()
        results_bic.loc['AR{}'.format(p), 'MA{}'.format(q)] = results.bic
    except:
        continue
results_bic = results_bic[results_bic.columns].astype(float)

fig, ax = plt.subplots(figsize=(10, 8))
ax = sns.heatmap(results_bic,
                 mask=results_bic.isnull(),
                 ax=ax,
                 annot=True,
                 fmt='.2f',
                 );
ax.set_title('BIC');

# Alternative model selection method, limited to only searching AR and MA parameters
from statsmodels.tsa.stattools import arma_order_select_ic
 
train_results = arma_order_select_ic(ts_train, ic=['aic', 'bic'], trend='nc', max_ar=4, max_ma=4)
 
print('AIC', train_results.aic_min_order)
print('BIC', train_results.bic_min_order)

AIC (1, 3)
BIC (2, 0)

結果：得出兩個不同的標準，比較尷尬，還需要進行篩選

模型殘差檢驗：
ARIMA模型的殘差是否是平均值為0且方差為常數的正態分佈
QQ圖：線性即正態分佈

#殘差分析 正態分佈 QQ圖線性
model_results.plot_diagnostics(figsize=(20, 10));#statsmodels庫

Q-Q圖：越像直線，則是正態分佈；越不是直線，離正態分佈越遠。

時間序列建模基本步驟：

1. 獲取被觀測系統時間序列資料；
2. 對資料繪圖，觀測是否為平穩時間序列；對於非平穩時間序列要先進行d階差分運算，化為平穩時間序列；
3. 經過第二步處理，已經得到平穩時間序列。要對平穩時間序列分別求得其自相關係數ACF 和偏自相關係數PACF ，通過對自相關圖和偏自相關圖的分析，得到最佳的階層 p 和階數 q
4. 由以上得到的 ，得到ARIMA模型。然後開始對得到的模型進行模型檢驗。