ARIMA初步研究

import pandas as pd
import numpy asnp
import matplotlib.pyplot as plt
from statsmodels.tsa.stattools import adfuller #导入ADF检验函数
from statsmodels.tsa.seasonalimport seasonal_decompose  #导入季节性分解函数，将数列分解为趋势、季节性和残差三部分
from statsmodels.stats.diagnosticimport acorr_ljungbox  #导入白噪声检验函数
from statsmodels.graphics.tsaplotsimport plot_pacf, plot_acf  #导入自相关和偏自相关的绘图函数
from matplotlib.tickerimport MaxNLocator #导入自动查找到最佳的最大刻度函数
from statsmodels.tsa.arima_modelimport ARIMA  #导入ARIMA模型
from sklearn.metrics import mean_squared_error


# 单位根检验(test_stationarity，ADF检验)，用于检验序列是否是平稳的
# 季节性分解函数(seasonal_decompose)，通过分解后的趋势、季节性确认序列是否是平稳的
# 白噪声检验函数
# 自相关性和偏自相关性(acf_pacf)，通过截尾或拖尾的lag值，初步确认p,q。也可以用来检验序列是否平稳

# ADF检验：这是一种检查数据稳定性的统计测试。无效假设：时间序列是不稳定的。
# 测试结果由测试统计量和一些置信区间的临界值组成。
# 如果“测试统计量”少于“临界值”，我们可以拒绝无效假设，并认为序列是稳定的。
# 当p-value<0.05，且TestStatistic显著小于Critical Value (5%)时，数列稳定
# 主要看p-value，显著小于的判断不精确
def test_stationarity(timeseries, window=12):
    rolmean = timeseries.rolling(window=window,center=False).mean()
    rolstd = timeseries.rolling(window=window,center=False).std()
    # 旧版方法，即将被移除
    # rolmean =pd.rolling_mean(timeseries, window=window)
    # rolstd = pd.rolling_std(timeseries,window=window)
    # 设置原始图，移动平均图和标准差图的式样
    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 indftest[4].items():
        dfoutput['Critical Value (%s)' % key] = value
    print(dfoutput)


# 按照滑动均值将维经过指数转化的时间序列分为趋势(trend),季节性(seasonality)和残差(residual)三部分
def decompose_plot(series, title=''):
    decomposition =seasonal_decompose(series)  # 季节性分解函数
    trend = decomposition.trend  #分解出的趋势，包含NaN值，原因不明
    seasonal = decomposition.seasonal  #分解出的季节性
    residual = decomposition.resid  #分解出的残差，包含NaN值，原因不明
    fig = decomposition.plot()
    fig.set_size_inches(12, 6)
    fig.suptitle(title)
    fig.tight_layout()
    fig2 = acf_pacf(residual, title='Residuals', figsize=(12, 6)) # 分析残差的自相关，偏自相关
    #test_stationarity(residual.dropna())  #Dropna后才能做稳定性检验
    # 原数据的残差一般是不稳定的，经过差分后的数据残差可能是平稳的


# 定义一个画自相关，偏自相关的函数
# series 输入的时间序列
# lags 自相关和偏自相关函数的滞后取值范围
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()


# #生成一个伪随机白噪声用于测试acorr_ljungbox的可靠性
# from random import gauss
# from random import seed
# from pandas import Series
# # seed random number generator
# # seed(1) #指定生成伪随机数的种子，指定后每次生成的随机数相同
# # create white noise series
# whitenoise = Series([gauss(0.0, 1.0) for i in range(1000)])#创建一个高斯分布的白噪声
# acf_pacf(whitenoise, lags=40, title='White Noise Series')
# print(u'白噪声检验结果为：', acorr_ljungbox(whitenoise,lags=1))#检验结果：平稳度，p-value。p-value>0.05为白噪声

# 欧元汇率预测
df_fx_data =pd.read_csv('D:\\BOE-XUDLERD.csv') # 读取数据集
df_fx_data['Date'] = pd.to_datetime(df_fx_data['Date']) # Date列转换为日期型
indexed_df =df_fx_data.set_index('Date') # 将日期设置为索引
ts = indexed_df['Value'] # 将DataFrame切片为一维Series，Series包含日期索引
# plt.plot(ts)#绘制折线图
# plt.show()
ts_week =ts.resample('W').mean() #时间重取样
# plt.plot(ts_week)
# plt.show()
ts_week = ts_week[1900:]
# # 原数据分析结果：不平稳，非白噪声，自相关和偏自相关不截尾，无法做时间序列建模
# # 平稳性测试
# test_stationarity(ts_week)
# decompose_plot(ts_week, title='ts_week')
# # 是否为白噪声测试,p-value<0.05非白噪声。平稳且非白噪声函数，可用于时间序列建模
# print(u'白噪声检验结果为：', acorr_ljungbox(ts_week, lags=1))
# # 自相关和偏自相关初步确认p,q
# acf_pacf(ts_week, lags=20, title='Weekly Euro Dollar Exchange Rate')

# # 一阶差分分析结果：平稳，非白噪声，自相关和偏自相关lag=1时截尾，可以做时间序列建模
# ts_week_diff = ts_week - ts_week.shift(1)
# test_stationarity(ts_week_diff.dropna())
# decompose_plot(ts_week_diff.dropna(), title='ts_week_diff')
# print(u'白噪声检验结果为：',acorr_ljungbox(ts_week_diff.dropna(), lags=1))
# # 自相关和偏自相关初步确认p,q
# acf_pacf(ts_week_diff.dropna(), lags=20, title='Weekly Euro Dollar ExchangeRate')

# log后一阶差分分析结果：平稳，非白噪声，自相关和偏自相关lag=1时截尾，可以做时间序列建模
# log后是否更平稳，比较Test Statistic和p-value，越小越平稳(负数直接比较，不是比较绝对值)
ts_week_log =np.log(ts_week)
ts_week_log_diff = ts_week_log.diff(1)#ts_week_log - ts_week_log.shift(1)
test_stationarity(ts_week_log_diff.dropna())
decompose_plot(ts_week_log_diff.dropna(), title='ts_week_log_diff')
print(u'白噪声检验结果为：',acorr_ljungbox(ts_week_log_diff.dropna(),lags=1))
# 自相关和偏自相关初步确认p,q
acf_pacf(ts_week_log_diff.dropna(),lags=20,title='WeeklyEuro Dollar Exchange Rate')
# 这里log后比log前稍微平稳一些

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

# # 一阶差分后平稳，所以d=1。p,q参数使用循环暴力调参确定最佳值
# d = 1
# pdq_result = pd.DataFrame()
# for p in np.arange(0, 7):
#     for q in np.arange(0, 7):
#         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:
#             print('NG')
# print(pdq_result)
# # 比较拟合的数据和原始数据
model =ARIMA(ts_week_log, order=(2,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比较
# 拟合的数据和原始数据ts_week_log的残差
residuals=pd.DataFrame(results_ARIMA.resid)
# # 计算拟合后残差的平方和(RSS)
# plt.title('RSS: %.f'% sum(results_ARIMA.fittedvalues-ts_week_log_diff)**2)
plt.show()
#计算拟合后残差的自相关和偏自相关
acf_pacf(residuals,lags=10,title='Residuals')#等于acf_pacf(results_ARIMA.fittedvalues-ts_week_log_diff,lags=10,title='Residuals')
# 因为ARIMA模型会使拟合和原始数据的残差变为白噪声，所以残差的核密度(概率密度)为正态分布
residuals.plot(kind='kde')#概率密度图kde
#打印残差的特征
print ('\n\nResiduals Summary:\n', residuals.describe())

# 将拟合的值逆向转换并可视化
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()

# 预测结果
pdq = (2, 1,2)
test_point = 50  #测试点数
train_size = int(len(ts_week_log) - test_point) #测试集大小
train, test =ts_week_log[:train_size], ts_week_log[train_size:]  #切分数据集
model2 = ARIMA(train,order=pdq).fit()
predictions2 = model2.forecast(test_point) # 连续预测N个值
# predictions2=model2.predict(start=train.size,end=train.size+50, dynamic=True)
# forecast返回值为有3个元素的元组(tuple)，每个元素都是一个array，说明：forecast : array， stderr : array，conf_int : array2D

# predictions没有Index，无法在图上显示出来，必须先追加Index
# predictions_series1 = pd.Series(predictions1, index=test.index)
predictions_series2= pd.Series(predictions2[0],index=test.index)
#a=np.exp(predictions_series1)
b =np.exp(predictions_series2)
plt.figure(figsize=(20,8))
plt.title('SpotExchange Rate, Euro into USD')
plt.plot(ts_week[-(test_point+100):], 'o',label='observed')
#plt.plot(np.exp(predictions_series1), 'g', label='one timepredict steps=50',color='red')  #一次性预测50步的方法
plt.plot(np.exp(predictions_series2),'g',label='rollingone-step out-of-sample forecast',color='green') # 每次预测下一步的方法
plt.legend(loc='upperleft')
plt.show()
# fig, ax =plt.subplots(figsize=(10, 6))
# ax.set(title='Spot Exchange Rate, Euro into USD', xlabel='Date', ylabel='Eurointo USD')
# ax.plot(ts_week[-test_point:], 'o', label='observed')
# ax.plot(np.exp(predictions_series), 'g', label='rolling one-stepout-of-sample forecast')#须将预测做指数转换
# legend = ax.legend(loc='upper left')
print('end')
# 指数加权移动平均法是很受欢迎的方法，所有的权重被指定给先前的值连同衰减系数。这可以通过pandas实现：
# AIC负数时比较方法，例如-3好过-1