| 理论加实践，终于把时间序列预测ARIMA模型讲明白了 | 您所在的位置：网站首页 › 时间序列的预测程序简述 › 理论加实践，终于把时间序列预测ARIMA模型讲明白了 |
理论加实践，终于把时间序列预测ARIMA模型讲明白了
|
上篇我们一起学习了一些关于时间序列预测的知识。而本文将通过一段时间内电力负荷波动的数据集来实战演示完整的ARIMA模型的建模及参数选择过程，其中包括数据准备、随机性、稳定性检验。本文旨在实践中学习，在实战过程中穿插理论知识梳理和学习，相信大家一定有所收获。 本文主要内容
时间序列建模基本步骤
获取被观测系统时间序列数据。 对数据绘图，观测是否为平稳时间序列；对于非平稳时间序列要先进行 阶差分运算，化为平稳时间序列。 经过第二步处理，已经得到平稳时间序列。要对平稳时间序列分别求得其自相关系数ACF 和偏自相关系数PACF ，通过对自相关图和偏自相关图的分析，得到最佳的阶层 和阶数 。 由以上得到的 、、，得到ARIMA模型。然后开始对得到的模型进行模型检验。
ARIMA模型介绍
ARIMA 模型[1]是一种流行且广泛使用的时间序列预测统计方法。 ARIMA 是代表autoRegressive I integrated Moving a average[2]自回归综合移动平均线的首字母缩写词，它是一类在时间序列数据中捕获一组不同标准时间结构的模型。预测方程中平稳序列的滞后称为“自回归”项，预测误差的滞后称为“移动平均”项，需要差分才能使其平稳的时间序列被称为平稳序列的“综合”版本。随机游走和随机趋势模型、自回归模型和指数平滑模型都是 ARIMA 模型的特例。 ARIMA 模型可以被视为一个“过滤器”，它试图将信号与噪声分开，然后将信号外推到未来以获得预测。ARIMA模型特别适合于拟合显示非平稳性的数据。 一般概念为了能够使用ARIMA，你需要了解一些概念。 平稳性从统计学的角度来看，平稳性是指数据的分布在时间上平移时不发生变化。因此，非平稳数据显示了由于趋势而产生的波动，必须对其进行转换才能进行分析。例如，季节性会导致数据的波动，并可以通过“季节性差异”过程消除。 差分从统计学的角度来看，数据差分是指将非平稳数据转换为平稳的过程，去除其非恒定的趋势。“差分消除了时间序列水平的变化，消除了趋势和季节性，从而稳定了时间序列的平均值。” 季节性差分应用于季节性时间序列以去除季节性成分。 ARIMA模型拆解剖析ARIMA的各个部分，以便更好地理解它如何帮助我们时间序列建模，并对其进行预测。 AR - 自回归自回归模型，顾名思义，就是及时地“回顾”过去，分析数据中先前的值，并对它们做出假设。这些先前的值称为“滞后”。一个例子是显示每月铅笔销售的数据。每个月的销售总额将被认为是数据集中的一个“进化变量”。这个模型是作为“利益的演化变量根据其自身的滞后值(即先验值)进行回归”而建立的。 I - 表示综合与类似的“ARMA”模型相反，ARIMA中的“I”指的是它的综合方面。当应用差分步骤时，数据是“综合”的，以消除非平稳性。表示原始观测值的差异，以允许时间序列变得平稳，即数据值被数据值和以前的值之间的差异替换。 MA - 移动平均线该模型的移动平均方面，是将观测值与应用于滞后观测值的移动平均模型的残差之间的相关性合并。 ARIMA用于使模型尽可能地符合时间序列数据的特殊形式。 ARIMA模型建立 一般步骤① 首先需要对观测值序列进行平稳性检测，如果不平稳，则对其进行差分运算直到差分后的数据平稳；② 在数据平稳后则对其进行白噪声检验，白噪声是指零均值常方差的随机平稳序列；③ 如果是平稳非白噪声序列就计算ACF（自相关系数）、PACF（偏自相关系数），进行ARMA等模型识别；④ 对已识别好的模型，确定模型参数，最后应用预测并进行误差分析。 一般地，对于给定的时间序列 ，平稳序列的建模过程可以用下图中的流程图表示。 
ARIMA实战剖析
导入statmodelsPython库已使用ARIMA模型。 import os import warnings import matplotlib.pyplot as plt import numpy as np import pandas as pd import datetime as dt import math from pandas.plotting import autocorrelation_plot from statsmodels.tsa.statespace.sarimax import SARIMAX from sklearn.preprocessing import MinMaxScaler from common.utils import load_data, mape from IPython.display import Image from statsmodels.graphics.tsaplots import plot_acf, plot_pacf from statsmodels.tsa.stattools import adfuller # adf检验库 from statsmodels.stats.diagnostic import acorr_ljungbox # 随机性检验库 from statsmodels.tsa.arima_model import ARMA %matplotlib inline plt.rcParams['figure.figsize'] = (12,6) pd.options.display.float_format = '{:,.2f}'.format np.set_printoptions(precision=2) warnings.filterwarnings("ignore") # specify to ignore warning messages 导入数据 energy = pd.read_csv('./data/energy.csv') energy.head(10)
绘制从2012年1月到2014年12月的所有可用能源数据。看到这些数据，并不陌生，因为在之前的文章中已经展示了部分数据。 energy.plot(y='load', subplots=True, figsize=(15, 8), fontsize=12) plt.xlabel('timestamp', fontsize=12) plt.ylabel('load', fontsize=12) plt.show()
划分训练和测试数据集
现在已经加载了数据，可以将其划分为训练集和测试集。要在训练集上训练模型。通常，在模型完成训练后，将使用测试集评估它的准确性。需要确保测试集涵盖了来自训练集的较晚时间段，以确保模型不会从未来时间段获取信息。 从2014年9月1日到10月31日，分配两个月的时间给训练集。测试集将包括2014年11月1日至12月31日两个月的时间段: train_start_dt = '2014-11-01 00:00:00' test_start_dt = '2014-12-30 00:00:00'由于这一数据反映的是每日能源消费，因此存在强烈的季节性模式，但当前消费与最近几天的消费规律最为相似。 可视化差异为了更加直观地看出训练集和测试集的差异，我们在同一张图中用不同颜色区分两个测试集，蓝色为训练集、橙色为测试集。 energy[(energy.index = train_start_dt)][['load']].rename(columns={'load':'train'}) \ .join(energy[test_start_dt:][['load']].rename(columns={'load':'test'}), how='outer') \ .plot(y=['train', 'test'], figsize=(15, 8), fontsize=12) plt.xlabel('timestamp', fontsize=12) plt.ylabel('load', fontsize=12) plt.show()
使用一个相对较小的时间窗口来训练数据就足够了。 准备训练数据现在需要通过对数据进行筛选和归一化来为模型训练准备数据。筛选需要的时间段和列的数据，并且对其进行归一化，其作用的是将数据投影在0-1之间。 ① 过滤原始数据集，只包括前面提到的每个set的时间段，只包括所需的列'load'加上日期索引。 train = energy.copy()[(energy.index >= train_start_dt) & (energy.index = test_start_dt][['load']] print('Training data shape: ', train.shape) print('Test data shape: ', test.shape) Training data shape: (1416, 1) Test data shape: (48, 1)② 使用MinMaxScaler()对训练数据进行 (0, 1) 标准化。 scaler = MinMaxScaler() train['load'] = scaler.fit_transform(train) train.head(10)
③ 原始数据和标准化数据进行可视化比较。 energy[(energy.index >= train_start_dt) & (energy.index 1): eval_df['APE'] = (eval_df['prediction'] - eval_df['actual']).abs() / eval_df['actual'] print(eval_df.groupby('h')['APE'].mean())② 计算一步的MAPE: print('One step forecast MAPE: ', (mape(eval_df[eval_df['h'] == 't+1']['prediction'], eval_df[eval_df['h'] == 't+1']['actual']) )*100, '%') One step forecast MAPE: 0.5570581332313952 %③ 打印多步预测MAPE: print('Multi-step forecast MAPE: ', mape(eval_df['prediction'], eval_df['actual'])*100, '%') Multi-step forecast MAPE: 1.1460048657704118 %结果值较低是很好的：考虑到一个MAPE为10的预测会下降10%。 ④ 为更加容易直观地看到这种精度测量，把他们可视化出来。 if(HORIZON == 1): ## Plotting single step forecast eval_df.plot(x='timestamp', y=['actual', 'prediction'], style=['r', 'b'], figsize=(15, 8)) else: ## Plotting multi step forecast plot_df = eval_df[(eval_df.h=='t+1')][['timestamp', 'actual']] for t in range(1, HORIZON+1): plot_df['t+'+str(t)] = eval_df[(eval_df.h=='t+'+str(t))]['prediction'].values fig = plt.figure(figsize=(15, 8)) ax = plt.plot(plot_df['timestamp'], plot_df['actual'], color='red', linewidth=4.0) ax = fig.add_subplot(111) for t in range(1, HORIZON+1): x = plot_df['timestamp'][(t-1):] y = plot_df['t+'+str(t)][0:len(x)] ax.plot(x, y, color='blue', linewidth=4*math.pow(.9,t), alpha=math.pow(0.8,t)) ax.legend(loc='best') plt.xlabel('timestamp', fontsize=12) plt.ylabel('load', fontsize=12) plt.show()
综上所述，这个过程的步骤如下： 模型识别。使用绘图和汇总统计来识别趋势、季节性和自回归元素，以了解所需的差异量和滞后大小。 参数估计。使用拟合程序找到回归模型的系数。 模型检查。使用残差的绘图和统计检验来确定模型未捕获的时间结构的数量和类型。 重复该过程，直到在样本内或样本外观察（例如训练或测试数据集）上达到理想的拟合水平。
网格搜索选择超参数
将网格搜索定义为一个函数evaluate_arima_model()，该函数以时间序列数据集作为输入，以及元组(p,d,q)作为参数用于评估模型。 数据集分为两部分：初始训练数据集为 66％，测试数据集为剩余的 34％。 迭代测试集的每个时间步。一次迭代就可以训练一个模型，然后使用该模型对新数据进行预测。每次迭代都进行预测并存储在列表中。最后用测试集将所有预测值与预期值列表进行比较，并计算并返回均方误差分数。 # evaluate an ARIMA model for a given order (p,d,q) def evaluate_arima_model(X, arima_order): # prepare training dataset train_size = int(len(X) * 0.66) train, test = X[0:train_size], X[train_size:] history = [x for x in train] # make predictions predictions = list() for t in range(len(test)): model = ARIMA(history, order=arima_order) model_fit = model.fit(disp=0) yhat = model_fit.forecast()[0] predictions.append(yhat) history.append(test[t]) # calculate out of sample error error = mean_squared_error(test, predictions) return error继续定义一个evaluate_models()的函数，该函数为ARIMA指定(p,d,q)参数，并以网格循环到方式进行迭代。 # evaluate combinations of p, d and q values for an ARIMA model def evaluate_models(dataset, p_values, d_values, q_values): # 确保输入数据是浮点值（而不是整数或字符串） dataset = dataset.astype('float32') best_score, best_cfg = float("inf"), None for p in p_values: for d in d_values: for q in q_values: order = (p,d,q) try: mse = evaluate_arima_model(dataset, order) if mse ; best_score: best_score, best_cfg = mse, order print('ARIMA%s MSE=%.3f' % (order,mse)) except: continue print('Best ARIMA%s MSE=%.3f' % (best_cfg, best_score)) 参考资料[1] ARIMA 模型: https://people.duke.edu/~rnau/411arim.htm [2]autoRegressive I integrated Moving a average: https://wikipedia.org/wiki/Autoregressive_integrated_moving_average [3]AR自回归: https://wikipedia.org/wiki/Autoregressive_integrated_moving_average [4]综合: https://wikipedia.org/wiki/Order_of_integration [5]移动平均: https://wikipedia.org/wiki/Moving-average_model [6]ACF 和 PACF : https://people.duke.edu/~rnau/411arim3.htm [7]'pyramid'库: https://alkaline-ml.com/pmdarima/0.9.0/modules/generated/pyramid.arima.auto_arima.html [8]wikipedia: https://wikipedia.org/wiki/Mean_absolute_percentage_error 
往期精彩回顾
适合初学者入门人工智能的路线及资料下载机器学习及深度学习笔记等资料打印机器学习在线手册深度学习笔记专辑《统计学习方法》的代码复现专辑
AI基础下载黄海广老师《机器学习课程》视频课黄海广老师《机器学习课程》711页完整版课件
本站qq群554839127，加入微信群请扫码：
|
【本文地址】