使用ARIMA进行时间序列预测（Python）

栏目: Python · 发布时间: 7年前

内容简介：前几篇文章介绍了使用时间序列进行预测的大致思路及流程，今天又找到一篇使用ARIMA进行时间序列的文章，处理的流程与先前的有些差异，简单的翻译出来供参考。首先是加载需要使用到的模块：这里比较特别的是指定

前几篇文章介绍了使用时间序列进行预测的大致思路及流程，今天又找到一篇使用ARIMA进行时间序列的文章，处理的流程与先前的有些差异，简单的翻译出来供参考。

导入数据

首先是加载需要使用到的模块：

import itertools
import pandas as pd
import numpy as np
import statsmodels.api as sm
import matplotlib.pyplot as plt
plt.style.use('fivethirtyeight')

这里比较特别的是指定 matplotlib 的样式 。

使用的数据为“美国夏威夷莫纳罗亚火山天文台空气样本中大气二氧化碳数据”，该数据手机了1958年3月至2001年12月的二氧化碳样本数据。重要的是statsmodels中自带此部分数据，可通过如下方式获取：

data = sm.datasets.co2.load_pandas().data

返回如下错误：

---------------------------------------------------------------------------
TypeError                                 Traceback (most recent call last)
<ipython-input-2-2db2ec2e4244> in <module>()
----> 1 data = sm.datasets.co2.load_pandas().data
 
D:\ProgramData\Anaconda3\lib\site-packages\statsmodels\datasets\co2\data.py in load_pandas()
     64     index = pd.DatetimeIndex(start=data.data['date'][0].decode('utf-8'),
     65                              periods=len(data.data), format='%Y%m%d',
---> 66                              freq='W-SAT')
     67     dataset = pd.DataFrame(data.data['co2'], index=index, columns=['co2'])
     68     #NOTE: this is how I got the missing values in co2.csv
 
TypeError: __new__() got an unexpected keyword argument 'format'

经查询，发现此为statsmodels的一个Bug，该版本 已在Master分支中修复，但未在released的版本中修复 。可以采取的方案是使用源代码方式进行模块的重新安装。这里只是为了做演示，所以直接拷贝了statsmodels模块文件夹下的co2.csv文件进行手动加载。 

data = pd.read_csv('co2.csv', parse_dates=['date'], index_col='date')
data = data['co2'].resample('MS').mean()
data = data.fillna(data.bfill())
print(data.head())
data.plot(figsize=(12, 6))
plt.show()

使用ARIMA进行时间序列预测（Python）

这里涉及到的内容有：

从绘制的图形看，整体的数据呈上升趋势，并且存在一定的季节性波动。

为ARIMA模型选择最优参数

在先前的时间序列预测完整教程中，详细介绍了ARIMA模型及，对应的(p,d,q）参数，当时介绍的方案是采用ACF和PACF图来确定具体的参数值。这里采用了另外一种方法“网格搜索”，通过迭代不同参数的组合来获取最佳参数。 

# Define the p, d and q parameters to take any value between 0 and 2
p = d = q = range(0, 2)
 
# Generate all different combinations of p, q and q triplets
pdq = list(itertools.product(p, d, q))
print(pdq)
 
# Generate all different combinations of seasonal p, q and q triplets
seasonal_pdq = [(x[0], x[1], x[2], 12) for x in list(itertools.product(p, d, q))]
print(seasonal_pdq)

接下来我们就可以通过不同的参数组合来多ARIMA进行评估。这里会使用到SARIMAX()，在评估和比较不同参数的模型时，我们可以使用 AIC 值来评估。该值将在模型中直接返回。在这里AIC的值越小，则模型越优。 

import warnings
warnings.filterwarnings("ignore") # specify to ignore warning messages
 
for param in pdq:
    for param_seasonal in seasonal_pdq:
        try:
            model = sm.tsa.statespace.SARIMAX(data, order=param, seasonal_order=param_seasonal, enforce_stationarity=False, enforce_invertibility=False)
            results = model.fit()
            print('ARIMA{}x{} - AIC:{}'.format(param, param_seasonal, results.aic))
        except:
            continue

由于某些参数组合可能会产生数值错误，所以在这里显式地禁用了警告消息，上述代码执行后获得的结果：

ARIMA(0, 0, 0)x(0, 0, 0, 12) - AIC:7612.583429881011
ARIMA(0, 0, 0)x(0, 0, 1, 12) - AIC:6787.343624031752
ARIMA(0, 0, 0)x(0, 1, 0, 12) - AIC:1854.828234141261
ARIMA(0, 0, 0)x(0, 1, 1, 12) - AIC:1596.711172763898
ARIMA(0, 0, 0)x(1, 0, 0, 12) - AIC:1058.9388921320026
ARIMA(0, 0, 0)x(1, 0, 1, 12) - AIC:1056.2878557033246
ARIMA(0, 0, 0)x(1, 1, 0, 12) - AIC:1361.6578978072075
ARIMA(0, 0, 0)x(1, 1, 1, 12) - AIC:1044.7647912998186
ARIMA(0, 0, 1)x(0, 0, 0, 12) - AIC:6881.048755126321
ARIMA(0, 0, 1)x(0, 0, 1, 12) - AIC:6072.662327947137
ARIMA(0, 0, 1)x(0, 1, 0, 12) - AIC:1379.1941067330986
ARIMA(0, 0, 1)x(0, 1, 1, 12) - AIC:1241.4174716851521
ARIMA(0, 0, 1)x(1, 0, 0, 12) - AIC:1084.9453872064087
ARIMA(0, 0, 1)x(1, 0, 1, 12) - AIC:780.4316223926196
ARIMA(0, 0, 1)x(1, 1, 0, 12) - AIC:1119.5957893630987
ARIMA(0, 0, 1)x(1, 1, 1, 12) - AIC:807.0912988896297
ARIMA(0, 1, 0)x(0, 0, 0, 12) - AIC:1675.8086923024293
ARIMA(0, 1, 0)x(0, 0, 1, 12) - AIC:1240.2211199194085
ARIMA(0, 1, 0)x(0, 1, 0, 12) - AIC:633.4425586468699
ARIMA(0, 1, 0)x(0, 1, 1, 12) - AIC:337.793854858494
ARIMA(0, 1, 0)x(1, 0, 0, 12) - AIC:619.9501759055394
ARIMA(0, 1, 0)x(1, 0, 1, 12) - AIC:376.9283760898942
ARIMA(0, 1, 0)x(1, 1, 0, 12) - AIC:478.3296906672489
ARIMA(0, 1, 0)x(1, 1, 1, 12) - AIC:323.07764996953097
ARIMA(0, 1, 1)x(0, 0, 0, 12) - AIC:1371.187260233786
ARIMA(0, 1, 1)x(0, 0, 1, 12) - AIC:1101.8410734303056
ARIMA(0, 1, 1)x(0, 1, 0, 12) - AIC:587.9479709744935
ARIMA(0, 1, 1)x(0, 1, 1, 12) - AIC:302.49490008406354
ARIMA(0, 1, 1)x(1, 0, 0, 12) - AIC:584.4333533493711
ARIMA(0, 1, 1)x(1, 0, 1, 12) - AIC:337.1999050779197
ARIMA(0, 1, 1)x(1, 1, 0, 12) - AIC:433.0863608210692
ARIMA(0, 1, 1)x(1, 1, 1, 12) - AIC:281.5190053393452
ARIMA(1, 0, 0)x(0, 0, 0, 12) - AIC:1676.8881767362054
ARIMA(1, 0, 0)x(0, 0, 1, 12) - AIC:1241.9353188506311
ARIMA(1, 0, 0)x(0, 1, 0, 12) - AIC:624.2602350563734
ARIMA(1, 0, 0)x(0, 1, 1, 12) - AIC:341.2896611730041
ARIMA(1, 0, 0)x(1, 0, 0, 12) - AIC:579.3896901160358
ARIMA(1, 0, 0)x(1, 0, 1, 12) - AIC:370.591911050293
ARIMA(1, 0, 0)x(1, 1, 0, 12) - AIC:476.0500429124218
ARIMA(1, 0, 0)x(1, 1, 1, 12) - AIC:329.5844988945166
ARIMA(1, 0, 1)x(0, 0, 0, 12) - AIC:1372.6085881684037
ARIMA(1, 0, 1)x(0, 0, 1, 12) - AIC:1199.4888202249872
ARIMA(1, 0, 1)x(0, 1, 0, 12) - AIC:586.4485732044088
ARIMA(1, 0, 1)x(0, 1, 1, 12) - AIC:305.62738298712236
ARIMA(1, 0, 1)x(1, 0, 0, 12) - AIC:586.5467067203238
ARIMA(1, 0, 1)x(1, 0, 1, 12) - AIC:400.3606781886424
ARIMA(1, 0, 1)x(1, 1, 0, 12) - AIC:433.54694646294087
ARIMA(1, 0, 1)x(1, 1, 1, 12) - AIC:284.35964168784733
ARIMA(1, 1, 0)x(0, 0, 0, 12) - AIC:1324.311112732457
ARIMA(1, 1, 0)x(0, 0, 1, 12) - AIC:1060.935191442204
ARIMA(1, 1, 0)x(0, 1, 0, 12) - AIC:600.7412682874252
ARIMA(1, 1, 0)x(0, 1, 1, 12) - AIC:312.1329632683155
ARIMA(1, 1, 0)x(1, 0, 0, 12) - AIC:593.6637754853627
ARIMA(1, 1, 0)x(1, 0, 1, 12) - AIC:349.209146507186
ARIMA(1, 1, 0)x(1, 1, 0, 12) - AIC:440.1375884358338
ARIMA(1, 1, 0)x(1, 1, 1, 12) - AIC:293.56145595783397
ARIMA(1, 1, 1)x(0, 0, 0, 12) - AIC:1262.6545542448305
ARIMA(1, 1, 1)x(0, 0, 1, 12) - AIC:1052.0636724058634
ARIMA(1, 1, 1)x(0, 1, 0, 12) - AIC:581.3099935252751
ARIMA(1, 1, 1)x(0, 1, 1, 12) - AIC:295.9374058139739
ARIMA(1, 1, 1)x(1, 0, 0, 12) - AIC:576.8647112079228
ARIMA(1, 1, 1)x(1, 0, 1, 12) - AIC:327.90491219439946
ARIMA(1, 1, 1)x(1, 1, 0, 12) - AIC:444.12436865154666
ARIMA(1, 1, 1)x(1, 1, 1, 12) - AIC:277.7802201071108

从上面返回的结果中，我们发现ARIMA(1, 1, 1)x(1, 1, 1, 12)返回的AIC为277.78，值最小。因此，我们应该认为这是我们所考虑的所有模型中的最佳选择。

拟合ARIMA模型

通过“网格搜索”我们找到了最佳拟合模型的参数，接下里我们将最佳参数值输入到一个新的 SARIMAX 模型：

model = sm.tsa.statespace.SARIMAX(data, order=(1, 1, 1), seasonal_order=(1, 1, 1, 12), enforce_stationarity=False, enforce_invertibility=False)
results = model.fit()
print(results.summary().tables[1])

输出内容为：

==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
ar.L1          0.3182      0.092      3.442      0.001       0.137       0.499
ma.L1         -0.6254      0.077     -8.163      0.000      -0.776      -0.475
ar.S.L12       0.0010      0.001      1.732      0.083      -0.000       0.002
ma.S.L12      -0.8769      0.026    -33.812      0.000      -0.928      -0.826
sigma2         0.0972      0.004     22.632      0.000       0.089       0.106
==============================================================================

从 SARIMAX 输出结果得到的 summary 属性返回大量信息，但我们将把注意力集中在coef列上。在这里每列的P值都币0.05小或接近0.05，因此我们可以保留模型中的所有权重。在拟合季节 ARIMA 模型(以及其他相关的模型)时，重要的是要运行模型诊断，以确保模型所做的假设是正确的。 plot_diagnostics()允许我们快速生成模型诊断并调查任何异常行为。

results.plot_diagnostics(figsize=(12, 12))
plt.show()

使用ARIMA进行时间序列预测（Python）

这里我们主要要确保残差稳定，且平均分布为0。上面的模型诊断表明：

  • 右上角的红色KDE线和黄色N(0,1)线非常的接近（N(0,1)表示，平均值为0，标准差为1），表明残差是正态分布。
  • 左下角对的 QQ分位图 表明，从N(0,1)的标准正态分布抽取的样本，残差(蓝点)的有序分布服从线性趋势。 同样，这是一个强烈的迹象表明，残差是正态分布。
  • 随着时间的推移（左上图）残差不会显示任何明显的季节性，似乎是白噪声。这通过右下角的自相关（即相关图）证实了这一点，它表明时间序列的残差与其自身的滞后版本有很低的相关性。

这些观察结果使我们得出结论，我们的模型产生了令人满意的结果，可以帮助我们预测未来的数据。

备注：虽然我们得到了令人满意的拟合参数，但是ARIMA模型的一些参数还是可以改变，以改进我们的模型拟合。 例如，我们的网格搜索只考虑受限制的参数组合集，所以如果我们扩大网格搜索，我们可能会找到更好的模型。

验证预测

我们已经获得了我们时间序列的模型，现在可以用来产生预测。 我们首先将预测值与时间序列的实际值进行比较，这将有助于我们了解我们的预测的准确性。 get_prediction()和conf_int()属性允许我们获得时间序列预测的值和相关的置信区间。

pred = results.get_prediction(start=pd.to_datetime('1998-01-01'), dynamic=False)
pred_ci = pred.conf_int()

代码中设置从1998年1月开始预测，通过dynamic=False设置后续的每一个都是使用1998年1月之前的数据预测生产的。

我们可以绘制二氧化碳时间序列的实际值和预测值，以评估我们做得如何。

ax = data['1990':].plot(label='Observed',figsize=(12, 6))
pred.predicted_mean.plot(ax=ax, label='One-step ahead Forecast', alpha=.7)
 
ax.fill_between(pred_ci.index,
                pred_ci.iloc[:, 0],
                pred_ci.iloc[:, 1], color='k', alpha=.2)
 
ax.set_xlabel('Date')
ax.set_ylabel('CO2 Levels')
plt.legend()
 
plt.show()

使用ARIMA进行时间序列预测（Python）

总体而言，我们的预测值与真实值保持一致，呈现总体增长趋势。量化我们的预测的准确性也是有用的。 我们将使用MSE（均方误差）来评估预测的准确性。

data_forecasted = pred.predicted_mean
data_truth = data['1998-01-01':]
 
# Compute the mean square error
mse = ((data_forecasted - data_truth) ** 2).mean()
print('The Mean Squared Error of our forecasts is {}'.format(round(mse, 2)))

这里预测的MSE值为0.07，是接近0的非常低的值。0的MSE是估计器将以完美的精度预测参数的观测值，这将是一个理想的场景但通常不可能。在这种情况下，我们只使用时间序列中的信息到某一点，之后，使用先前预测时间点的值生成预测。然而，使用动态预测可以获得更好地表达我们的真实预测能力。 在下面的代码块中，我们指定从1998年1月起开始计算动态预测和置信区间。

pred_dynamic = results.get_prediction(start=pd.to_datetime('1998-01-01'), dynamic=True, full_results=True)
pred_dynamic_ci = pred_dynamic.conf_int()
 
ax = data['1990':].plot(label='Observed', figsize=(12, 6))
pred_dynamic.predicted_mean.plot(label='Dynamic Forecast', ax=ax)
 
ax.fill_between(pred_dynamic_ci.index,
                pred_dynamic_ci.iloc[:, 0],
                pred_dynamic_ci.iloc[:, 1], color='k', alpha=.25)
 
ax.fill_betweenx(ax.get_ylim(), pd.to_datetime('1998-01-01'), data.index[-1],
                 alpha=.1, zorder=-1)
 
ax.set_xlabel('Date')
ax.set_ylabel('CO2 Levels')
 
plt.legend()
plt.show()

使用ARIMA进行时间序列预测（Python）

我们看到即使使用动态预测，总体预测也是准确的。 所有预测值（红线）与实际值（蓝线）相当吻合，并且在我们预测的置信区间内。我们再次通过计算MSE量化我们预测的预测的准确性：

# Extract the predicted and true values of our time series
data_forecasted = pred_dynamic.predicted_mean
data_truth = data['1998-01-01':]
 
# Compute the mean square error
mse = ((data_forecasted - data_truth) ** 2).mean()
print('The Mean Squared Error of our forecasts is {}'.format(round(mse, 2)))

从动态预测获得的预测值产生的MSE 为1.01。要比非动态要高，这个也在预料之中，原因是可供参考的历史数据较少。

采用动态或非动态的预测都证实了时间序列模型的有效。然而，关于时间序列预测的大部分兴趣是能够及时预测未来价值观。

预测及可视化

在本教程的最后一步，我们将介绍如何利用季节性ARIMA时间序列模型来预测未来的价值。 我们的时间序列对象的get_forecast()属性可以计算预先指定数量的“步数”的预测值。

# Get forecast 500 steps ahead in future
pred_uc = results.get_forecast(steps=500)
 
# Get confidence intervals of forecasts
pred_ci = pred_uc.conf_int()

我们可以使用此代码的输出绘制其未来值的时间序列和预测。

ax = data.plot(label='Observed', figsize=(12, 6))
pred_uc.predicted_mean.plot(ax=ax, label='Forecast')
ax.fill_between(pred_ci.index,
                pred_ci.iloc[:, 0],
                pred_ci.iloc[:, 1], color='k', alpha=.25)
ax.set_xlabel('Date')
ax.set_ylabel('CO2 Levels')
 
plt.legend()
plt.show()

使用ARIMA进行时间序列预测（Python）

我们所作的预测和相关的置信区间现在都可以用来进一步理解时间序列，并预测将会发生什么。 我们的预测显示，预计时间序列将以稳定的速度持续增长。

总结

上面的代码中主要用到了statsmodels，对于如何使用好statsmodels还需要进一步学习。另外SARIMAX与ARIMA的差异也没有讲清楚。

原文链接：： https://www.digitalocean.com/community/tutorials/a-guide-to-time-series-forecasting-with-arima-in-python-3


以上就是本文的全部内容，希望对大家的学习有所帮助，也希望大家多多支持 码农网

查看所有标签

猜你喜欢:

本站部分资源来源于网络，本站转载出于传递更多信息之目的，版权归原作者或者来源机构所有，如转载稿涉及版权问题，请联系我们。

码农书籍
Agile Web Development with Rails 4

Agile Web Development with Rails 4

Sam Ruby、Dave Thomas、David Heinemeier Hansson / Pragmatic Bookshelf / 2013-10-11 / USD 43.95

Ruby on Rails helps you produce high-quality, beautiful-looking web applications quickly. You concentrate on creating the application, and Rails takes care of the details. Tens of thousands of deve......一起来看看 《Agile Web Development with Rails 4》 这本书的介绍吧!

码农工具
JSON 在线解析
JSON 在线解析

在线 JSON 格式化工具

图片转BASE64编码
图片转BASE64编码

在线图片转Base64编码工具

html转js在线工具
html转js在线工具

html转js在线工具

  • New
  • 文章
  • 话题
  • 教程
    • 关注 码农网 公众号
    版权所有，保留一切权利！© 2018-2026 码农网 粤ICP备17054400号-3