使用ARIMA预测欧元汇率

这篇文章主要是用于复习ARIMA，针对一些可能存在的问题进行进一步的理清。对预测的整个流程再做一次清晰的回顾。

导入数据

数据集: https://pan.baidu.com/s/1x5ug4x72CyFGc9TaTwuTKw 提取码: n76d

导入及观察数据：

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
data = pd.read_csv('BOE-XUDLERD.csv', parse_dates=['Date'],index_col='Date') # 读取欧元汇率数据集
ts = data['Value'] # 将DataFrame切片为一维Series，Series包含日期索引
ts.plot(figsize=(12,6))

ts_week =ts.resample('W').mean() # 时间重取样
ts_week.plot(figsize=(12,6))

ADF检验

在ARMA/ARIMA这样的自回归模型中，模型对时间序列数据的平稳是有要求的，因此，需要对数据或者数据的n阶差分进行平稳检验，而一种常见的方法就是ADF检验，即单位根检验。

平稳随机过程

在数学中，平稳随机过程（Stationary random process）或者严平稳随机过程（Strictly-sense stationary random process），又称狭义平稳过程，是在固定时间和位置的概率分布与所有时间和位置的概率分布相同的随机过程：即随机过程的统计特性不随时间的推移而变化。这样，数学期望和方差这些参数也不随时间和位置变化。平稳在理论上有严平稳和宽平稳两种，在实际应用上宽平稳使用较多。宽平稳的数学定义为：对于时间序列 yt，若对任意的t,k,m，满足：

则称时间序列 yt 是宽平稳的。

平稳是自回归模型ARMA的必要条件，因此对于时间序列，首先要保证应用自回归的n阶差分序列是平稳的。

单位根检验

单位根检验是指检验序列中是否存在单位根，因为存在单位根就是非平稳时间序列了。单位根就是指单位根过程，可以证明，序列中存在单位根过程就不平稳，会使回归分析中存在伪回归。而迪基-福勒检验（Dickey-Fuller test）和扩展迪基-福勒检验（Augmented Dickey-Fuller test可以测试一个自回归模型是否存在单位根（unit root）。

ADF Test

在 Python 中对时间序列的建模通常使用statsmodel库，在 statsmodels.tsa.stattools.adfuller 中可进行adf校验，一般传入一个1d 的array like的data就行，包括list，numpy array 和 pandas series都可以作为输入，其他参数可以保留默认。

单位根检验(test_stationarity，ADF检验)，用于检验序列是否是平稳的，测试结果由测试统计量和一些置信区间的临界值组成。如果“测试统计量”少于“临界值”，我们可以拒绝无效假设，并认为序列是稳定的。当p-value<0.05，且TestStatistic显著小于Critical Value (5%)时，数列稳定。

检测数据

from statsmodels.tsa.stattools import adfuller
def test_stationarity(timeseries, window=12):
    rolmean = timeseries.rolling(window=window, center=False).mean()
    rolstd = timeseries.rolling(window=window, center=False).std()
    orig = plt.plot(timeseries, color='blue', label='Original')# 设置原始图，移动平均图和标准差图的式样
    mean = plt.plot(rolmean, color='red', label='Rolling Mean')
    std = plt.plot(rolstd, color='black', label='Rolling Std')
    plt.legend(loc='best')  # 使用自动最佳的图例显示位置
    plt.title('Rolling Mean & Standard Deviation')
    plt.show() #供肉眼观察是否平稳
    print('ADF检验结果：')
    dftest = adfuller(timeseries, autolag='AIC')  # 使用减小AIC的办法估算ADF测试所需的滞后数
    # 将ADF测试结果、显著性概率、所用的滞后数和所用的观测数打印出来
    dfoutput = pd.Series(dftest[0:4], index=['Test Statistic', 'p-value', 'Num Lags Used', 'Num Observations Used'])
    for key, value in dftest[4].items():
        dfoutput['Critical Value (%s)' % key] = value
    print(dfoutput)
test_stationarity(ts_week)

从上面的检测结果可以看到并没有满足平稳性要求。

白噪声检测

对于纯随机序列，又称白噪声序列，序列的各项数值之间没有任何相关关系，序列在进行完全无序的随机波动，可以终止对该序列的分析。白噪声序列是没有信息可提取的平稳序列。Ljung-Box test是对randomness的检验，statsmodel库中 statsmodels.stats.diagnostic.acorr_ljungbox ：

#生成一个伪随机白噪声用于测试acorr_ljungbox的可靠性
from statsmodels.stats.diagnostic import acorr_ljungbox  #导入白噪声检验函数
from random import gauss
whitenoise = pd.Series([gauss(0.0, 1.0) for i in range(1000)])#创建一个高斯分布的白噪声
print(u'白噪声检验结果为：', acorr_ljungbox(whitenoise,lags=1))#检验结果：平稳度，p-value。p-value>0.05为白噪声

白噪声检验原假设：时间序列是白噪声。p-value值小于0.05，拒绝原假设，所以一阶差分后的序列为平稳非白噪声序列(95%的把握)。

对ts_week进行白噪声检测：

# 是否为白噪声测试,p-value<0.05非白噪声。平稳且非白噪声函数，可用于时间序列建模
print(u'白噪声检验结果为：', acorr_ljungbox(ts_week, lags=1))

白噪声检验结果为： (array([2281.72237192]), array([0.]))，可以看到其p-value值为0，所以并不是白噪音序列。

季节、趋势分解

statsmodels.tsa.seasonal.seasonal_decompose 可以将时间序列的数据分解为趋势(trend),季节性(seasonality)和残差(residual)三部分。

from statsmodels.tsa.seasonal import seasonal_decompose  #导入季节性分解函数，将数列分解为趋势、季节性和残差三部分
def decompose_plot(series):
    decomposition = seasonal_decompose(series)  # 季节性分解函数
    trend = decomposition.trend
    seasonal = decomposition.seasonal
    residual = decomposition.resid
    fig = decomposition.plot()
    fig.set_size_inches(12, 8)
    fig.tight_layout()
decompose_plot(ts_week)

ACF和PACF图

自相关函数（ACF）：这是时间序列和它自身滞后版本之间的相关性的测试。比如在自相关函数可以比较时间的瞬间‘t1’…’t2’以及序列的瞬间‘t1-5’…’t2-5’ (t1-5和t2 是结束点)。

部分自相关函数(PACF):这是时间序列和它自身滞后版本之间的相关性测试，但是是在预测（已经通过比较干预得到解释）的变量后。如：滞后值为5，它将检查相关性，但是会删除从滞后值1到4得到的结果。

根据序列的自相关函数和偏自相关函数的特征可以初步判断模型类型，如下表：

截尾就是在某阶之后，系数都为 0（或靠近0） 。 拖尾就是有一个衰减的趋势，但是不都为 0 。

相关方法：

• graphics.tsaplots.plot_acf
• graphics.tsaplots.plot_pacf

# 自相关性和偏自相关性(acf_pacf)，通过截尾或拖尾的lag值，初步确认p,q。也可以用来检验序列是否平稳
from statsmodels.graphics.tsaplots import plot_pacf, plot_acf  #导入自相关和偏自相关的绘图函数
from matplotlib.ticker import MaxNLocator # 导入自动查找到最佳的最大刻度函数
def acf_pacf(series, lags=40, title=None, figsize=(12, 6)):
    # 求自相关函数
    fig = plt.figure(figsize=figsize) # figure指设置图形的特征。figsize为设置图形大小，单位为inch
    ax1 = fig.add_subplot(211)  # 子图2行1列的第一张,大于10后用逗号分隔，如(3,4,10)
    ax1.set_xlabel('lags')  # 横坐标为滞后值
    ax1.set_ylabel('AutoCorrelation')  # 纵坐标为自相关系数
    ax1.xaxis.set_major_locator(MaxNLocator(integer=True))  # 设置主坐标轴为自动设置最佳的整数型坐标标签
    plot_acf(series.dropna(), lags=lags, ax=ax1)  # 没有title参数，需要删除title
    # 求偏自相关函数
    ax2 = fig.add_subplot(212)
    ax2.set_xlabel('lags')
    ax2.set_ylabel('Partial AutoCorrelation')
    ax2.xaxis.set_major_locator(MaxNLocator(integer=True))
    plot_pacf(series.dropna(), lags=lags, ax=ax2)  # 没有title参数，需要删除title
    plt.tight_layout()  #设置为紧凑布局
    plt.show()
# 自相关和偏自相关初步确认p,q
acf_pacf(ts_week, lags=40, title='Weekly Euro Dollar Exchange Rate')

数据处理及校验

对数据进行一阶差分去除趋势后进行数据校验：

ts_week_diff = ts_week - ts_week.shift(1)
test_stationarity(ts_week_diff.dropna())

print(u'白噪声检验结果为：',acorr_ljungbox(ts_week_diff.dropna(), lags=1))

白噪声检验结果为： (array([155.18954176]), array([1.27283835e-35]))

acf_pacf(ts_week_diff.dropna(), lags=20, title='Weekly Euro Dollar ExchangeRate')# 自相关和偏自相关初步确认p,q

一阶差分分析结果：平稳，非白噪声，自相关和偏自相关lag=1时截尾，可以做时间序列建模

对数据进行log后一阶差分分析

ts_week_log =np.log(ts_week)
ts_week_log_diff = ts_week_log.diff(1)

log后是否更平稳，比较Test Statistic和p-value，越小越平稳(负数直接比较，不是比较绝对值)

print(u'白噪声检验结果为：',acorr_ljungbox(ts_week_log_diff.dropna(), lags=1))

白噪声检验结果为： (array([142.03403004]), array([9.55961654e-33]))

acf_pacf(ts_week_log_diff.dropna(),lags=20, title='WeeklyEuro Dollar Exchange Rate')

确定ARIMA参数值

从上面的ACF和PACF可大致确定ARIMA的参数，p=偏自相关截尾的lag。如果偏自相关截尾后，自相关仍然拖尾，考虑增加p；q=自相关截尾的lag。如果自相关截尾后，偏自相关仍然拖尾，考虑增加q。以上仅是理论上的说法，具体最合适的p,q，还是要以AIC、BIC或预测结果残差最小来循环确定。ARIMA模型是通过优化AIC, HQIC, BIC等得出的, 一般来说AIC较小的模型拟合度较好,但拟合度较好不能说明预测能力好。

from statsmodels.tsa.arima_model import ARIMA  # 导入ARIMA模型
d = 1 # 一阶差分后平稳，所以d=1。p,q参数使用循环暴力调参确定最佳值
pdq_result = pd.DataFrame()
for p in range(0, 5):
    for q in range(0, 5):
        try:
            model = ARIMA(ts_week_log,order=(p, d, q))  # ARIMA的参数: order = p, d, q
            results_ARIMA = model.fit()
            pdq_result = pdq_result.append(pd.DataFrame(
                {'p': [p], 'd': [d],'q': [q], 'AIC': [results_ARIMA.aic], 'BIC': [results_ARIMA.bic],
                 'HQIC':[results_ARIMA.hqic]}))
            #print(results_ARIMA.summary())  #summary2()也可以使用
        except:
            pass
print(pdq_result)

   p  d  q           AIC           BIC          HQIC
0  0  1  0 -13939.315506 -13927.843765 -13935.132034
0  0  1  1 -14083.366353 -14066.158742 -14077.091145
0  0  1  2 -14083.421767 -14060.478286 -14075.054823
0  0  1  3 -14081.931844 -14053.252492 -14071.473164
0  0  1  4 -14080.862357 -14046.447135 -14068.311941
0  1  1  0 -14083.792730 -14066.585119 -14077.517522
0  1  1  1 -14083.699762 -14060.756281 -14075.332818
0  1  1  2 -14082.602763 -14053.923412 -14072.144084
0  1  1  4 -14078.946543 -14038.795451 -14064.304391
0  2  1  0 -14083.401327 -14060.457845 -14075.034383
0  2  1  1 -14081.986100 -14053.306749 -14071.527420
0  2  1  2 -14081.512486 -14047.097264 -14068.962070
0  3  1  0 -14082.610717 -14053.931365 -14072.152037
0  3  1  1 -14080.785531 -14046.370309 -14068.235115
0  3  1  2 -14087.953053 -14047.801961 -14073.310901
0  3  1  3 -14086.313032 -14040.426069 -14069.579144
0  4  1  0 -14080.937459 -14046.522237 -14068.387043
0  4  1  1 -14078.958559 -14038.807467 -14064.316408
0  4  1  2 -14086.257204 -14040.370241 -14069.523316
0  4  1  4 -14082.963600 -14025.604897 -14062.046240

从上面的结果看，(p,d,q)比较适合选择(3,1,2)。所以用(3,1,2)参数比较拟合的数据和原始数据：

model =ARIMA(ts_week_log, order=(3, 1, 2))  # ARIMA的参数: order = p, d, q
results_ARIMA =model.fit()
plt.plot(ts_week_log_diff)
plt.plot(results_ARIMA.fittedvalues, color='red')  # fittedvalues返回的是d次差分后的序列, 因此和ts_week_log_diff比较
residuals=pd.DataFrame(results_ARIMA.resid) # 拟合的数据和原始数据ts_week_log的残差
plt.show()

计算拟合后残差的自相关和偏自相关

acf_pacf(residuals,lags=10,title='Residuals') #等于acf_pacf(results_ARIMA.fittedvalues-ts_week_log_diff,lags=10,title='Residuals')

上图中我们看到完全的没有拖尾了。因为ARIMA模型会使拟合和原始数据的残差变为白噪声，所以残差的核密度(概率密度)为正态分布：

print('Residuals Summary:\n', residuals.describe())
residuals.plot(kind='kde') #概率密度图kde

将拟合的值逆向转化并可视化：

predictions_ARIMA_diff= pd.Series(results_ARIMA.fittedvalues, copy=True)  # 复制Log和差分后的预测值
predictions_ARIMA_diff_cumsum= predictions_ARIMA_diff.cumsum()  # 一阶差分的逆向转换为累计和
predictions_ARIMA_log= pd.Series(ts_week_log.iloc[0], index=ts_week_log.index)  #将一个新序列的每个元素都填充为ts_week_log的第一个元素
predictions_ARIMA_log= predictions_ARIMA_log.add(predictions_ARIMA_diff_cumsum,fill_value=0)  # 将累计和追加到第一个元素上，还原log序列。其中的空值填充为0
predictions_ARIMA= np.exp(predictions_ARIMA_log)  # 进行指数运算，还原为原始值
plt.plot(ts_week)
plt.plot(predictions_ARIMA)
plt.title('RMSE:%.4f' % np.sqrt(sum((predictions_ARIMA - ts_week) ** 2) / len(ts_week)))
plt.show()

预测结果

相关文档： statsmodels.tsa.arima_model.ARIMA

from sklearn.metrics import mean_squared_error
from statsmodels.tsa.arima_model import ARIMA
 
pdq = (3, 1, 2)
test_point = 15 # 测试点数
size = len(ts_week_log) - 15
train, test = ts_week_log[0:size], ts_week_log[size:len(ts_week_log)] #切分数据集
 
model = ARIMA(train, order=pdq)
model_fit = model.fit()
predictions = model_fit.predict(start=size,end=size+test_point-1, dynamic=True).tolist()
 
error = mean_squared_error(test, predictions)
print('Printing Mean Squared Error of Predictions...')
 
predictions_series = pd.Series(predictions, index = test.index) # predictions没有Index，无法在图上显示出来，必须先追加Index
 
fig, ax = plt.subplots()
ax.set(title='Spot Exchange Rate, Euro into USD', xlabel='Date', ylabel='Euro into USD')
ax.plot(ts_week[-50:], 'o', label='observed')
ax.plot(np.exp(euro_predictions_series), 'g', label='rolling one-step out-of-sample forecast')
legend = ax.legend(loc='upper left')
legend.get_frame().set_facecolor('w')

参考链接： https://dataplatform.cloud.ibm.com/exchange/public/entry/view/815137c868b916821dec777bdc23013c