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')