預(yù)測包括，數(shù)值擬合，線性回歸，多元回歸，時(shí)間序列，神經(jīng)網(wǎng)絡(luò)等等

對于單變量的時(shí)間序列預(yù)測：模型有AR,MA,ARMA,ARIMA，綜合來說用ARIMA即可表示全部。

數(shù)據(jù)和代碼鏈接：數(shù)據(jù)和Jupyter文件

以預(yù)測美國未來10年GDP的變換情況為列：

目錄

第一步進(jìn)行數(shù)據(jù)導(dǎo)入

第二步進(jìn)行平穩(wěn)序列分析

?第三步進(jìn)行不平穩(wěn)序列的差分運(yùn)算

第四步進(jìn)行模型定階和模型選擇及擬合

?第五步進(jìn)行模型結(jié)果分析和模型檢驗(yàn)?

第六步進(jìn)行模型預(yù)測

PS：自動化AUTO-ARIMA的比較

ARIMA流程圖如下：

第一步進(jìn)行數(shù)據(jù)導(dǎo)入：

? ? 觀察要預(yù)測變量的數(shù)據(jù)情況，去除雜數(shù)據(jù)，得到數(shù)據(jù)本身的屬性如：大小，類型等

import pandas as pd,numpy as np
from matplotlib import pyplot as plt

#加載數(shù)據(jù)，sheet_name指定excel表的數(shù)據(jù)頁面，header指定指標(biāo)column屬性，loc去除雜數(shù)據(jù)，可選：parse_dates=['']，index_col='',use_cols=['']
df=pd.read_excel('./data/time1.xls',sheet_name='數(shù)據(jù)',header=1).loc[1:,:]

#DateFrame索引重置
df=df.set_index('DATE')#df.set_index('DATE',inplace=True)

#查看前5行
print(df.head(5))

#查看列索引
print(df.columns)#print(df.keys())

#查看表的維度
print(df.shape)

#查看行索引
print(df.index)

#np.array()  array1.reshape(,)  df.values.astype(int).tolist() np.vstack((a1,a2))  np.hstack((a1,a2))  round()  iloc

#時(shí)間索引拆分
# dates=pd.date_range(start='1991-01-01',end='2007-08-01',freq='MS')#日期取值和格式轉(zhuǎn)換，MS代表每月第一天
# years=[d.strftime('%Y-%m') for d in dates][0:200:25]
# years.append('2007-09')

結(jié)果如下：

第二步進(jìn)行平穩(wěn)序列分析：

? ? 我們研究分析的時(shí)間序列，即面板數(shù)據(jù)，只有是寬平穩(wěn)的才有研究意義，如果是非平穩(wěn)序列需要將其差分轉(zhuǎn)換為平穩(wěn)序列才能進(jìn)行分析，對于嚴(yán)平穩(wěn)序列，性質(zhì)不變化，即序列為白噪聲序列，這樣的序列沒有研究意義。
? ??
? ? 所以這里拿到GDP的時(shí)間序列數(shù)據(jù)后，先進(jìn)行原序列的白噪聲檢驗(yàn)，利用LB統(tǒng)計(jì)量，若p值小于顯著水平a=0.05，即認(rèn)為原序列為非白噪聲序列，有研究的意義。
? ??
? ? 接著畫時(shí)序圖，主觀上判斷GDP時(shí)間序列的平穩(wěn)情況，若呈現(xiàn)明顯的趨勢性，則為非平穩(wěn)序列；若沒有明顯趨勢，或者出現(xiàn)你不好判斷平穩(wěn)的情況，這時(shí)候你需要借助ADF單位根統(tǒng)計(jì)量來進(jìn)行輔助判斷它是否平穩(wěn)。對于ADF統(tǒng)計(jì)量，你可以比較統(tǒng)計(jì)量的值來判斷，也可以比較p值來判斷，如顯著水平a=0.05時(shí)候，如果你變量的ADF統(tǒng)計(jì)量小于a=0.05對于的ADF統(tǒng)計(jì)量，則說明該變量對應(yīng)的時(shí)間序列為平穩(wěn)序列；還一種更為直接的方法是，只需要判斷該變量的p值如果小于0.05，則也說明為平穩(wěn)序列。
? ?

plt.rcParams['axes.unicode_minus'] = False  # 解決保存圖像是負(fù)號'-'顯示為方塊的問題
plt.rc('font',family='SimHei')
plt.style.use('ggplot')

df.plot(secondary_y=['CS','INV','P_GDP','GOV_NET'])#單個指標(biāo)時(shí)序圖 df['CS'].plot()
plt.xlabel('Date')
plt.ylabel('Value')
plt.title('Time Series Plot')
#plt.grid()
plt.show()

from statsmodels.tsa.stattools import adfuller# ADF檢驗(yàn)

for i in df.columns:
    data=df[i]
    print(f'{i}的單位根檢驗(yàn)')
    result = adfuller(data)#默認(rèn)情況下，regression參數(shù)為'c'，表示使用包含截距項(xiàng)的回歸模型。
    print('ADF Statistic: %f' % result[0])#ADF統(tǒng)計(jì)量
    print('p-value: %f' % result[1])#p值
    print('Critical Values:')#在置信水平下的臨界值
    for key, value in result[4].items():
        print('\t%s: %.3f' % (key, value))
    print()

代碼結(jié)果：

第三步進(jìn)行不平穩(wěn)序列的差分運(yùn)算：

? ? 一般時(shí)間序列的差分不會超過三階，對原始數(shù)據(jù)進(jìn)行差分運(yùn)算，重復(fù)第二步的時(shí)序圖和ADF單位根檢驗(yàn)操作，若發(fā)現(xiàn)差分后的序列為平穩(wěn)序列，則記下差分的階數(shù)，這里GDP的差分階數(shù)1，通過了平穩(wěn)序列檢驗(yàn)。再次對差分后的序列進(jìn)行白噪聲檢驗(yàn)，若通過為非白噪聲序列，則進(jìn)行下一步操作。
? ??

#差分時(shí)序圖
diff_data = df.diff(periods=1).dropna()# 創(chuàng)建一階差分的時(shí)間序列，加上dropna()后續(xù)不需要執(zhí)行[1:]
#print(diff_data)

diff_data.plot()
plt.xlabel('Date')
plt.ylabel('Value')
plt.title('Time Series Plot')
#plt.grid()
plt.show()

for i in diff_data.columns:
    data=diff_data[i]#Series索引選取，一階差分第一個數(shù)據(jù)為NA
    
    print(f'{i}的單位根檢驗(yàn)')
    result = adfuller(data)#結(jié)果對應(yīng)位置的數(shù)據(jù)需要自己判斷是什么含義
    print('ADF Statistic: %f' % result[0])#ADF統(tǒng)計(jì)量
    print('p-value: %f' % result[1])#p值
    print('Critical Values:')#在置信水平下的臨界值
    for key, value in result[4].items():
        print('\t%s: %.3f' % (key, value))
    print()

from statsmodels.stats.diagnostic import acorr_ljungbox

#這里為一階差分后的平穩(wěn)序列進(jìn)行白噪聲檢驗(yàn)，lags為1，否則lags為0，這里拿上述的GDP指標(biāo)進(jìn)行
lags = [1,4,8,16,32]

print('差分序列的白噪聲檢驗(yàn)結(jié)果為：'+'\n',acorr_ljungbox(diff_data['GDP'], lags))  # 返回統(tǒng)計(jì)量和p值,這里的lags對應(yīng)差分階數(shù)

print('原始數(shù)據(jù)序列的白噪聲檢驗(yàn)結(jié)果為：'+'\n',acorr_ljungbox(df['GDP'], lags))

代碼結(jié)果如下：

第四步進(jìn)行模型定階和模型選擇及擬合：

? ? 模型定階先采用自相關(guān)圖和偏自相關(guān)圖進(jìn)行初步判斷：
? ? ? ? ? ?自相關(guān)系數(shù) ? ? ? ? ? ? ? ?偏自相關(guān)系數(shù) ? ? ? ? ? ? 差分
? ? AR ? ? ? ? ? ? 拖尾 ? ? ? ? ? ? ? ? ? ? ? ? ?p階截尾 ? ? ? ? ? ? ? ? ? ? ?0
? ? MA ? ? ? ? ? ?q階截尾 ? ? ? ? ? ? ? ? ? ?拖尾 ? ? ? ? ? ? ? ? ? ? ? ? ? ?0
? ? ARMA ? ? ? 拖尾 ? ? ? ? ? ? ? ? ? ? ? ? ?拖尾 ? ? ? ? ? ? ? ? ? ? ? ? ? ?0
? ? ARIMA ? ? ?拖尾 ? ? ? ? ? ? ? ? ? ? ? ? ?拖尾 ? ? ? ? ? ? ? ? ? ? ? ? ? ?d
? ? 總之，上述的模型都可以用ARIMA(p,d,q)來實(shí)現(xiàn)，AR,MA,ARMA都是ARIMA的特殊情況。
? ??
? ? 對于拖尾和截尾的經(jīng)驗(yàn)判斷：
? ? 拖尾：負(fù)指數(shù)單調(diào)收斂到0，或者呈現(xiàn)余弦式衰減
? ? 截尾：迅速衰減到0，并且在0附近波動
? ? 截尾比拖尾趨近于0更加迅速，截尾在后期不會再有明顯的增加
? ??
? ? ? ?

????????在利用圖表進(jìn)行初步判斷后，依舊無法拿準(zhǔn)模型的階數(shù)，可以利用AIC和BIC準(zhǔn)則進(jìn)行輔助判斷，或者叫做優(yōu)化。比如我判斷模型的階數(shù)可能為2，1，4，或者3，1，3，這時(shí)候我拿不準(zhǔn)，你可以先做ARIMA(2,1,4)，如果發(fā)現(xiàn)后續(xù)的模型顯著性檢驗(yàn)和參數(shù)檢驗(yàn)都通過了，那么你可以直接利用ARIMA(2,1,4)進(jìn)行預(yù)測并得到結(jié)論，但你或許也去嘗試附近的幾個階數(shù)，發(fā)ARIMA(3,1,3),ARIMA(2,1,3)等等也通過了，但是ARIMA(3,1,3)的AIC和BIC的值是最小的，即該模型的擬合精度更高，是一個相對最優(yōu)的模型，這一步便是模型的優(yōu)化。

? ??
? ? 這邊的定階和優(yōu)化在python里的實(shí)現(xiàn)是：從q為0到（數(shù)據(jù)長度/10），從p為0到（數(shù)據(jù)長度/10），一般階數(shù)不會超過（數(shù)據(jù)長度/10）；然后進(jìn)行q*p次模型擬合，利用模型擬合的結(jié)果比較AIC和BIC的值，BIC值越小越小的那個模型就是相對最優(yōu)模型，再用該模型的階數(shù)作為最終模型擬合的階數(shù)。
? ??

#生成自相關(guān)圖，偏自相關(guān)圖
from statsmodels.graphics.tsaplots import plot_acf,plot_pacf

fig, ax = plt.subplots(figsize=(5, 4))
plot_acf(diff_data['GDP'],ax=ax)#可以換成df的數(shù)據(jù)，這里用一階差分?jǐn)?shù)據(jù)得到平穩(wěn)序列
    
ax.set_title("Autocorrelation Plot")
lags=list(range(24))
ax.set_xticks(lags)#這里把x坐標(biāo)刻度變更精細(xì)，加網(wǎng)格圖更方便,xticklabels替換標(biāo)簽
ax.set_xlabel("Lag")
ax.set_ylabel("Autocorrelation")
ax.grid(alpha=0.5)
plt.legend(['GDP'])#這里區(qū)別直接放入df.columns[i]，如果是多字符如‘CS’，這樣會被認(rèn)為是一個序列，拆成C和S的圖例
plt.show()
    
fig, ax = plt.subplots(figsize=(5, 4))
lags=list(range(24))
ax.set_xticks(lags)#這里把x坐標(biāo)刻度變更精細(xì)，加網(wǎng)格圖更方便,xticklabels替換標(biāo)簽
plot_pacf(diff_data['GDP'], ax=ax,method='ywm')#ywm替換默認(rèn)的yw，去除警告
ax.set_xlabel('Lags')
ax.set_ylabel('Partial Autocorrelation')
ax.grid(alpha=0.5)
plt.legend(['GDP'])
plt.show()

import warnings
warnings.filterwarnings("ignore")
import statsmodels.api as sm


#diff_data['GDP'].values.astype(float),這里發(fā)現(xiàn)Serise的dtype為object，模型用的應(yīng)該為float或者int類型，需要注意原數(shù)據(jù)的數(shù)據(jù)類型是否一致

Min=float('inf')
for i in range(0,6):#AIC，BIC最小找到p，q階數(shù)來定階，從0開始定階是否可行？？
    for j in range(0,6):
        result=sm.tsa.ARIMA(df['GDP'].values.astype(float),order=(i,1,j)).fit()
        print([i,j,result.aic,result.bic])
        if result.bic<Min:
            Min=result.bic
            best_pq=[i,j,result.aic,result.bic]
print(f'最優(yōu)定階為{best_pq}')

代碼結(jié)果：

第五步進(jìn)行模型結(jié)果分析和模型檢驗(yàn)?

? ? 模型檢驗(yàn)分為模型參數(shù)檢驗(yàn)和模型顯著性檢驗(yàn)。

? ? 模型顯著性檢驗(yàn)：即殘差的白噪聲檢驗(yàn)，若殘差為白噪聲序列，即原序列信息提取充分，看LB統(tǒng)計(jì)量，python中模型結(jié)果的LB統(tǒng)計(jì)量即為殘差的LB統(tǒng)計(jì)量，若p值小于0.05，則為非白噪聲序列，若p值大于0.05，則說明殘差為白噪聲序列，這是我們想要的結(jié)果。這里你可以單獨(dú)拿出模型的殘差值，自己繪制殘差的時(shí)序圖，QQ圖，正態(tài)分布圖，或者自己進(jìn)行白噪聲檢驗(yàn)輔助判斷。白噪聲序列即服從正態(tài)分布，時(shí)序圖是平穩(wěn)波動,QQ圖上的數(shù)值點(diǎn)在對角線附近。
? ??
? ? 模型的參數(shù)檢驗(yàn)：檢驗(yàn)每一個未知參數(shù)是否顯著為0，檢驗(yàn)?zāi)Ｐ褪欠褡罹?。參?shù)不顯著非零，則可以從擬合模型從剔除，看t統(tǒng)計(jì)量。
? ??
? ? 模型結(jié)果部分解釋：
? ? const：常數(shù)項(xiàng) ? ? ??
? ? ar.L1：自回歸項(xiàng)系數(shù)
? ? ma.L1：移動平均項(xiàng)系數(shù)
? ? sigma2：方差
? ? 各個參數(shù)下的P>|z|這一行，若小于a=0.05則拒絕原假設(shè)，認(rèn)為參數(shù)顯著非0，即不需要簡化模型。
? ? Ljung-Box：LB統(tǒng)計(jì)量,這里需要注意的是LB的P值需要>0.95,即判斷殘差序列為白噪聲序列，小于0.05為非白噪聲序列。
? ? Jarque-Bera：JB統(tǒng)計(jì)量 ??
? ? Heteroskedasticity檢驗(yàn)的結(jié)果表明方差穩(wěn)定情況

Result=sm.tsa.ARIMA(df['GDP'].values.astype(float),order=(best_pq[0],1,best_pq[1])).fit()
print(Result.summary())  #顯示模型的所有信息

print(len(Result.resid))
#print(Result.resid)這里觀察到殘差的第一項(xiàng)為原數(shù)據(jù)的1239.5，即差分?jǐn)?shù)據(jù)不管第一項(xiàng)，這里需要調(diào)整殘差的觀測

#這里就可以觀察到原始模型的結(jié)果LB統(tǒng)計(jì)量和這里的白噪聲檢驗(yàn)是一致的，p>0.05,即認(rèn)為殘差為白噪聲序列，原序列信息提取充分。
lags = [1,4,8,16,32]
print('差分序列的白噪聲檢驗(yàn)結(jié)果為：'+'\n',acorr_ljungbox(Result.resid[1:], lags))

## 查看模型的擬合殘差分布
fig = plt.figure(figsize=(12,5))
ax = fig.add_subplot(1,2,1)
plt.plot(Result.resid[1:])
plt.title("ARIMA(2,1,1)殘差曲線")

## 檢查殘差是否符合正太分布
ax = fig.add_subplot(1,2,2)
sm.qqplot(Result.resid[1:], line='q', ax=ax)
plt.title("ARIMA(2,1,0)殘差Q-Q圖")
plt.tight_layout()
plt.show()

fig = plt.figure(figsize=(12,5))
Residual=pd.DataFrame(Result.resid[1:])
Residual.plot(kind='kde', title='密度')
plt.legend('')
plt.show()

代碼結(jié)果為：

??

第六步進(jìn)行模型預(yù)測：

? ? 模型預(yù)測即利用上述的相對最優(yōu)模型進(jìn)行原始數(shù)據(jù)變量后續(xù)時(shí)間對應(yīng)的值的預(yù)測。
? ? 注意在python中，predict函數(shù)如果是對差分?jǐn)?shù)據(jù)進(jìn)行預(yù)測，注意start的起始和end的中止情況，比如一階差分的數(shù)據(jù)是從第二個觀測值
? ? 開始進(jìn)行，對應(yīng)的殘差也是如此，在后續(xù)繪圖的時(shí)候需要注意到差分模型繪圖的情況。

#預(yù)測，繪制原序列和預(yù)測序列值對比圖

Predict=Result.predict(start=1, end=len(df['GDP'])-1+1+10); #不加參數(shù)默認(rèn)0到n-1，要加預(yù)測個數(shù)在end后面N-1+預(yù)測n即可
#如果是一階差分的序列預(yù)測，第一個數(shù)據(jù)已經(jīng)差分消去了，應(yīng)該start從第二個觀測數(shù)據(jù)開始，即n=1；如果是0階，則不需要按默認(rèn)0到n-1

print(list(zip(range(193,203),Predict[-10:])))#打印預(yù)測值

plt.figure()
plt.plot(range(193),df['GDP'].values)#'o-k'
plt.plot(range(193+10),Predict)#'P--'
plt.legend(('原始觀測值','預(yù)測值'))
plt.xticks(list(range(0,203,10)),rotation=90)
plt.show()

plt.figure()
plt.plot(range(193),df['GDP'].values)#'o-k'
plt.plot(range(192,193+10),Predict[-11:])#'P--'#接著原數(shù)據(jù)最后一個，進(jìn)行擬合預(yù)測表示
plt.legend(('原始觀測值','預(yù)測值'))
plt.xticks(list(range(0,203,10)),rotation=90)
plt.show()

代碼結(jié)果：

PS：自動化AUTO-ARIMA的比較

import pmdarima as pm
# ## 自動搜索合適的參數(shù)
model = pm.auto_arima(df['GDP'].values,
                      start_p=1, start_q=1, # p,q的開始值
                      max_p=12, max_q=12, # 最大的p和q
                      d = 0,            # 尋找ARMA模型參數(shù)
                      m=1,              # 序列的周期
                      seasonal=False,   # 沒有季節(jié)性趨勢
                      trace=True,error_action='ignore',  
                      suppress_warnings=True, stepwise=True)
 
print(model.summary())

這一塊個人沒有深入研究，auto的方法有優(yōu)點(diǎn)也有缺點(diǎn)，?也算提供一個寫代碼的思路。文章來源地址http://www.zghlxwxcb.cn/news/detail-533015.html

到了這里，關(guān)于一文搞懂Python時(shí)間序列預(yù)測（步驟，模板，python代碼）的文章就介紹完了。如果您還想了解更多內(nèi)容，請?jiān)谟疑辖撬阉鱐OY模板網(wǎng)以前的文章或繼續(xù)瀏覽下面的相關(guān)文章，希望大家以后多多支持TOY模板網(wǎng)！