1. 程式人生 > >時間序列相關演算法與分析步驟

時間序列相關演算法與分析步驟

1.純隨機序列(白噪聲序列),這時候可以停止分析,因為就像預測下一次硬幣哪一面朝上一樣毫無規律。

2.平穩非白噪聲序列,它們的均值和方差是常數,對於這類序列,有成熟的模型來擬合這個序列在未來的發展狀況,如AR,MA,ARMA等(具體模型演算法及實現在後面)

3.非平穩序列,一般做法是把他們轉化為平穩的序列,在按照平穩序列的演算法進行擬合。如果經過差分後平穩,則應使用ARIMA模型進行擬合。

注:本文模型採用的資料為某餐廳一個多月內的銷量資料,包含兩個特徵:時間和銷量

Q1:序列的平穩性用什麼來衡量呢?

方法1:

根據時序圖和自相關圖的特徵做出主觀的判斷,如下圖:
時序圖:
這裡寫圖片描述
自相關圖:
這裡寫圖片描述
從上圖可以基本看出,自相關係數的絕對值長期都保持了較大的值,所以可以判斷上述時間序列存在自相關性。

平穩的序列自相關圖和偏自相關圖不是拖尾就是截尾。

截尾就是在某階之後,係數都為 0 。
拖尾就是有一個衰減的趨勢,但是不都為 0 。

從自相關圖來看,呈現三角對稱形式,不存在截尾或拖尾,屬於單調序列的典型表現形式,原始資料屬於不平穩序列。

注:

  • 如果自相關是拖尾,偏相關截尾,則用 AR 演算法

  • 如果自相關截尾,偏相關拖尾,則用 MA 演算法

  • 如果自相關和偏相關都是拖尾,則用 ARMA 演算法, ARIMA 是 ARMA 演算法的擴充套件版,用法類似 。

相關係數的計算方法:
這裡寫圖片描述
VAR表示方差

方法2:

根據單位根檢驗

如果存在單位根,則此序列為隨機非平穩序列

Q2:平穩序列應該怎麼分析呢?

目前最常用的擬合平穩序列的模型為ARMA(Autoregressive moving average)模型,全稱是自迴歸移動平均模型,他又可以分為AR模型,MA模型和ARMA模型三大類。

1.自迴歸AR(p)模型

這裡寫圖片描述
自迴歸模型描述的是當前值與歷史值之間的關係。

2.移動平均MA(q)模型

這裡寫圖片描述
移動平均模型描述的是自迴歸部分的誤差累計。

3.ARMA(p,q)模型

ARMA(p,q)模型中包含了p個自迴歸項和q個移動平均項,ARMA(p,q)模型可以表示為:
這裡寫圖片描述

當q=0時,是AR(p)模型
當p=0時,是MA(q)模型

一般分析步驟:
這裡寫圖片描述

Q3:非平穩序列怎麼分析呢?

從上面的模型中可以看出,如果是非平穩序列,我們需要先把它轉為平穩序列之後再進行分析。

一般我們使用ARIMA(Autoregressive Integrated Moving Average model)進行分析

ARIMA(p,d,q)中,AR是”自迴歸”,p為自迴歸項數;MA為”滑動平均”,q為滑動平均項數,d為使之成為平穩序列所做的差分次數(階數)

“差分”一詞雖未出現在ARIMA的英文名稱中,卻是關鍵步驟。

Q4:舉個栗子看下唄!

讀取資料

#-*- coding: utf-8 -*-
#arima時序模型

import pandas as pd

#引數初始化
discfile = '../data/arima_data.xls'
forecastnum = 5

#讀取資料,指定日期列為指標,Pandas自動將“日期”列識別為Datetime格式
data = pd.read_excel(discfile, index_col = u'日期')

自相關檢測


#時序圖
import matplotlib.pyplot as plt
plt.rcParams['font.sans-serif'] = ['SimHei'] #用來正常顯示中文標籤
plt.rcParams['axes.unicode_minus'] = False #用來正常顯示負號
data.plot()
plt.show()

#自相關圖
from statsmodels.graphics.tsaplots import plot_acf
plot_acf(data).show()

#平穩性檢測
from statsmodels.tsa.stattools import adfuller as ADF
print(u'原始序列的ADF檢驗結果為:', ADF(data[u'銷量']))
#返回值依次為adf、pvalue、usedlag、nobs、critical values、icbest、regresults、

自相關圖
這裡寫圖片描述
可以看出自相關係數的絕對值長期都保持很大,多以基本判斷存在自相關性。

ADF檢測結果p值顯著大於0.05(p=0.9983),最終判斷為非平穩序列

一階差分後繼續檢測

#差分後的結果
D_data = data.diff().dropna()
D_data.columns = [u'銷量差分']
D_data.plot() #時序圖
plt.show()
plot_acf(D_data).show() #自相關圖
from statsmodels.graphics.tsaplots import plot_pacf
plot_pacf(D_data).show() #偏自相關圖
print(u'差分序列的ADF檢驗結果為:', ADF(D_data[u'銷量差分'])) #平穩性檢測

#白噪聲檢驗
from statsmodels.stats.diagnostic import acorr_ljungbox
print(u'差分序列的白噪聲檢驗結果為:', acorr_ljungbox(D_data, lags=1)) #返回統計量和p值

這裡寫圖片描述
上圖是差分後的銷量結果

這裡寫圖片描述
自相關圖顯示出1階截尾的性質
這裡寫圖片描述
偏自相關圖顯示出1階拖尾的性質

從ADF的結果(p=0.0226)和自相關圖以及偏自相關圖中可以看出一階差分後的序列是平穩的非白噪聲序列。

給ARIMA模型定階
從一階差分後的序列是平穩的非白噪聲序列可以看出ARIMA模型中的d=1

定階方法:
1.人為判斷:自相關圖顯示出從第1階之後的截尾性質,偏自相關圖從第1階之後顯示出拖尾的性質,所以人為判斷使用MA(1)模型,即ARMA(0,1,1)
2.相對最優模型識別,當p和q均小於等於3的所有組合的BIC資訊量,取其中BIC資訊量達到最小的模型階數。

#定階
pmax = int(len(D_data)/10) #一般階數不超過length/10
qmax = int(len(D_data)/10) #一般階數不超過length/10
bic_matrix = [] #bic矩陣
for p in range(pmax+1):
  tmp = []
  for q in range(qmax+1):
    try: #存在部分報錯,所以用try來跳過報錯。
      tmp.append(ARIMA(data, (p,1,q)).fit().bic)
    except:
      tmp.append(None)
  bic_matrix.append(tmp)

bic_matrix = pd.DataFrame(bic_matrix) #從中可以找出最小值

p,q = bic_matrix.stack().idxmin() #先用stack展平,然後用idxmin找出最小值位置。
print(u'BIC最小的p值和q值為:%
  • 17

BIC矩陣
取其中BIC資訊量達到最小的模型階數。
這裡寫圖片描述
確定p=0,q=1

擬合模型

model = ARIMA(data, (p,1,q)).fit() #建立ARIMA(0, 1, 1)模型
model.summary2() #給出一份模型報告
model.forecast(5) #作為期5天的預測,返回預測結果、標準誤差、置信區間。

最終得到模型的預測結果

資料和完整程式碼可以通過在留言中留下郵箱獲取哦~