python調用statsmodels模塊實現整合移動平均自回歸模型(ARIMA)——以預測股票收盤價為例

目錄

• 程序簡介
• 程序/數據集下載
• 代碼分析

程序簡介

調用statsmodels模塊對上證指數的收盤價進行ARIMA模型動態建模，ARIMA適合短期預測，因此輸入為15個數據，輸出為1個數據
程序輸入：原序列，需要往后預測的個數
程序輸出：預測序列，模型結構(白噪聲檢驗、單根檢驗、一階差分自相關圖、一階差分偏自相關圖)

差分整合移動平均自回歸模型(ARIMA)，ARIMA(p，d，q)中，AR是”自回歸”，p為自回歸項數；MA為”滑動平均”，q為滑動平均項數，d為使之成為平穩序列所做的差分次數（階數）。

程序/數據集下載

點擊進入下載地址

代碼分析

導入模塊、路徑

# -*- coding: utf-8 -*-
from Module.BuildModel import ARIMA
from sklearn.metrics import mean_absolute_error
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import os

#路徑目錄
baseDir = ''#當前目錄
staticDir = os.path.join(baseDir,'Static')#靜態文件目錄
resultDir = os.path.join(baseDir,'Result')#結果文件目錄

讀取數據、分成訓練集和測試集，查看前5行

#讀取數據
data = pd.read_csv(staticDir+'/000001.csv',encoding='gbk')
train = data['收盤價'].values[-20:-10]#訓練數據
test = data['收盤價'].values[-10:]#測試數據
data.head()

	日期	股票代碼	名稱	收盤價	最高價	最低價	開盤價	前收盤	漲跌額	漲跌幅	成交量	成交金額
0	2020-02-18	'000001	上證指數	2984.9716	2990.6003	2960.7751	2981.4097	2983.6224	1.3492	0.0452	311665913	3.74998562648e+11
1	2020-02-17	'000001	上證指數	2983.6224	2983.6371	2924.9913	2924.9913	2917.0077	66.6147	2.2837	313198007	3.67014340129e+11
2	2020-02-14	'000001	上證指數	2917.0077	2926.9427	2899.5739	2899.8659	2906.0735	10.9342	0.3763	250650627	3.08080368726e+11
3	2020-02-13	'000001	上證指數	2906.0735	2935.4060	2901.2425	2927.1443	2926.8991	-20.8256	-0.7115	274804844	3.34526327364e+11
4	2020-02-12	'000001	上證指數	2926.8991	2926.8991	2892.4240	2895.5561	2901.6744	25.2247	0.8693	248733429	2.97534420493e+11

調用ARIMA函數進行動態建模，函數自動繪制一階差分自相關圖和一階差分偏自相關圖，即每次預測1個值，持續多次，計算並打印MAE

其中ARIMA函數位置在Module/BuildModel.py，代碼會在下文給出

#ARIMA動態建模
yPre = []
for i in range(test.shape[0]):
    #只預測1個數
    result = ARIMA(train,1)
    yPre.append(result['predict']['value'][0])
    #更新訓練集
    train = train.tolist()[:-1]
    train.append(test[i])
    train = np.array(train).reshape(-1)
#計算MAE
MAE = mean_absolute_error(test,yPre)
print('誤差MAE',MAE)

這里給出上文用到的ARIMA函數，直接運行是直接驗證代碼有效性的另一個簡單實驗，程序會保留白噪聲檢驗、單根檢驗、一階差分自相關圖、一階差分偏自相圖

# -*- coding: utf-8 -*-
import os
import pandas as pd
import numpy as np
import statsmodels.api as sm
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.stattools import adfuller
from statsmodels.tsa import arima_model
from statsmodels.stats.diagnostic import acorr_ljungbox
import warnings
warnings.filterwarnings("ignore")


def ARIMA(series,n):
    '''
    只討論一階差分的ARIMA模型，預測，數字索引從1開始
    series:時間序列
    n:需要往后預測的個數
    '''
    series = np.array(series)
    series = pd.Series(series.reshape(-1))
    currentDir = os.getcwd()#當前工作路徑
    #一階差分數據
    fd = series.diff(1)[1:]
    plot_acf(fd).savefig(currentDir+'/一階差分自相關圖.png')
    plot_pacf(fd).savefig(currentDir+'/一階差分偏自相關圖.png')
    #一階差分單位根檢驗
    unitP = adfuller(fd)[1]
    if unitP>0.05:
        unitAssess = '單位根檢驗中p值為%.2f，大於0.05，該一階差分序列可能為非平穩序列'%(unitP)
        #print('單位根檢驗中p值為%.2f，大於0.05，認為該一階差分序列判斷為非平穩序列'%(unitP))
    else:
        unitAssess = '單位根檢驗中p值為%.2f，小於0.05，認為該一階差分序列為平穩序列'%(unitP)
        #print('單位根檢驗中p值為%.2f，小於0.05，認為該一階差分序列判斷為平穩序列'%(unitP))
    #白噪聲檢驗
    noiseP = acorr_ljungbox(fd, lags=1)[-1]
    if noiseP<=0.05:
        noiseAssess = '白噪聲檢驗中p值為%.2f，小於0.05，認為該一階差分序列為非白噪聲'%noiseP
        #print('白噪聲檢驗中p值為%.2f，小於0.05，認為該一階差分序列為非白噪聲'%noiseP)
    else:
        noiseAssess = '白噪聲檢驗中p值%.2f，大於0.05，該一階差分序列可能為白噪聲'%noiseP
        #print('白噪聲檢驗中%.2f，大於0.05，認為該一階差分序列為白噪聲'%noiseP)
    #BIC准則確定p、q值
    pMax = int(series.shape[0]/10)# 一般階數不超過length/10
    qMax = pMax# 一般階數不超過length/10
    bics = list()
    for p in range(pMax + 1):
        tmp = list()
        for q in range(qMax + 1):
            try:
                tmp.append(arima_model.ARIMA(series, (p, 1, q)).fit().bic)
            except Exception as e:
                #print(str(e))
                tmp.append(1e+10)#加入一個很大的數
        bics.append(tmp)
    bics = pd.DataFrame(bics)
    p, q = bics.stack().idxmin()
    #print('BIC准則下確定p,q為%s,%s'%(p,q))
    #建模
    model = arima_model.ARIMA(series,order=(p, 1, q)).fit()
    predict = model.forecast(n)[0]
    return {
            'model':{'value':model,'desc':'模型'},
            'unitP':{'value':unitP,'desc':unitAssess},
            'noiseP':{'value':noiseP[0],'desc':noiseAssess},
            'p':{'value':p,'desc':'AR模型階數'},
            'q':{'value':q,'desc':'MA模型階數'},
            'params':{'value':model.params,'desc':'模型系數'},
            'predict':{'value':predict,'desc':'往后預測%d個的序列'%(n)}
            }
    
if __name__ == "__main__":
    data = data = np.array([1.2,2.2,3.1,4.5,5.6,6.7,7.1,8.2,9.6,10.6,11,12.4,13.5,14.7,15.2])
    x = data[0:10]#輸入數據
    y = data[10:]#需要預測的數據
    result = ARIMA(x,len(y))#預測結果,一階差分偏自相關圖,一階差分自相關圖
    predict = result['predict']['value']
    predict = np.round(predict,1)
    print('真實值:',y)
    print('預測值:',predict)
    print(result)

打印建立的股票模型假設檢驗結論，可以看到序列可能為白噪聲，但是仍然可以強行建模

#打印模型結構
print(result['noiseP']['desc'])
print(result['unitP']['desc'])

白噪聲檢驗中p值0.98，大於0.05，該一階差分序列可能為白噪聲
單位根檢驗中p值為0.99，大於0.05，該一階差分序列可能為非平穩序列

觀測值預測值對比可視化

#可視化
plt.clf()#清理歷史繪圖
#用來正常顯示中文標簽
plt.rcParams['font.sans-serif']=['SimHei'] 
#用來正常顯示負號
plt.rcParams['axes.unicode_minus']=False
plt.plot(range(test.shape[0]),yPre,label="預測值")
plt.plot(range(test.shape[0]),test,label="觀測值")
plt.legend()
plt.title('ARIMA預測效果，MAE：%2f'%MAE)
plt.savefig(resultDir+'/ARIMA預測效果.png',dpi=100,bbox_inches='tight')