程序简介
调用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\')