多元时间序列数据的预测和建模
23,808 阅读
加载中...
加载中...
时间是决定企业是上升还是下降的最关键因素。 这就是为什么我们看到商店和电子商务平台的销售与节日保持一致。 这些企业分析了多年的支出数据，以了解打开大门并看到消费者支出增加的最佳时机。
但作为数据科学家，您如何进行此分析？ 别担心，你不需要建立时间机器！ 时间序列建模是一种强大的技术，可作为理解和预测趋势和模式的门户。
但即使是时间序列模型也有不同的方面。 我们在网上看到的大多数例子都涉及单变量时间序列。 不幸的是，现实世界的用例并不像那样。 有多个变量在起作用，同时处理所有这些变量是数据科学家将获得其价值的地方。
在本文中，我们将了解多元变量时间序列是什么，以及如何处理它。 我们还将进行案例研究并在Python中实现它，以便您对该主题有实际的了解。
本文假设您熟悉单变量时间序列，其属性以及用于预测的各种技术。 由于本文将重点关注多变量时间序列，我建议您阅读以下文章，作为单变量时间序列的一个很好的介绍：
但在进入多变量时间序列的详细信息之前，我会给你快速回顾一下单变量时间序列。 让我们逐一看一下它们来理解它们之间的区别。
顾名思义，单变量时间序列是一个具有单个时间相关变量的系列。
例如，查看下面的样本数据集，其中包含过去2年的温度值（每小时）。 这里，温度是因变量（取决于时间）。
如果要求我们预测接下来几天的温度，我们将查看过去的值并尝试测量和提取一种模式。 我们会注意到早上和晚上的温度较低，而下午达到峰值。 此外，如果您有过去几年的数据，您会发现它在11月到1月期间更冷，而在4月到6月则相对较热。
这些观察将有助于我们预测未来的价值观。 您是否注意到我们只使用了一个变量（过去2年的温度）？ 因此，这称为单变量时间序列分析/预测。
多变量时间序列具有多个时间相关变量。 每个变量不仅取决于其过去的值，而且还依赖于其他变量。 此依赖关系用于预测未来值。 听起来很复杂？ 让我解释。
考虑上面的例子。 现在假设我们的数据集包括出汗百分比、露点(dew point，在一定压力下出 现第一个液滴时的温度)、风速、云量百分比等，以及过去两年的温度值。 在这种情况下，需要考虑多个变量来最佳地预测温度。 像这样的系列将属于多变量时间序列的范畴。 以下是对此的说明：
现在我们了解了多变量时间序列的样子，让我们了解如何使用它来构建预测。
在本节中，我将向您介绍多变量时间序列预测中最常用的方法之一 ——即向量自回归模型（Vector Auto Regression，VAR）。
在VAR模型中，每个变量都是过去自身值和所有其他变量的过去值的线性函数。 为了更好地解释这一点，我将使用一个简单的可视化示例：
我们有两个变量，y1和y2。 我们需要在时间t预测这两个变量的值，从过去n值的给定数据开始。 为简单起见，我认为滞后值为1。
为了计算y1（t），我们将使用y1和y2的过去值。 类似地，为了计算y2（t），将使用y1和y2的过去值。 下面是一种表示这种关系的简单数学方法：
y_1(t) = a_1 + w_{11} \times y_1(t-1) + w_{12} \times y_2(t-1) + e_1(t-1)
y_2(t) = a_2 + w_{21} \times y_1(t-1) + w_{22} \times y_2(t-1) + e_1(t-1)
其中，$a_1$和$a_2$是常数项，$w_{*}$是系数矩阵中的系数，$e_{*}$是误差项。
这些方程类似于AR过程的方程。 由于AR过程用于单变量时间序列数据，因此未来值仅是它们自己的过去值的线性组合。 考虑AR（1）过程：
y(t) = a+w\times y(t-1) +e
在这个例子中，我们只有一个变量$y$，一个常数项$a$，一个误差项$e$和一个系数$w$。为了适应VAR的每个等式中的多个变量项，我们将使用向量。 我们可以用以下形式写出等式（1）和（2）：
上述形式扩展到多维之后：
y(t) = a+ w_1\times y(t-1) + \cdots + w_p\times y(t-p) + \epsilon t
上述等式中的$\epsilon t$表示多变量白噪声。 对于多变量时间序列，$\epsilon t$应该是满足以下条件的连续随机向量：
回想一下我们之前看到的温度预测的例子。 可以对其进行论证，将其视为多个单变量系列。 我们可以使用像AR这样的简单单变量预测方法来解决它。 由于目标是预测温度，我们可以简单地删除其他变量（温度除外）并在剩余的单变量系列上拟合模型。
另一个简单的想法是使用我们已知的技术单独预测每个系列的值。 这将使工作变得非常简单！ 那你为什么要学习另一种预测技术呢？ 这个话题不够复杂吗？
从上面的等式（1）和（2）可以清楚地看出，每个变量都使用每个变量的过去值来进行预测。 与AR不同，VAR能够理解和使用几个变量之间的关系。 这对于描述数据的动态行为非常有用，并且还可以提供更好的预测结果。 此外，实现VAR就像使用任何其他单变量技术一样简单（您将在上一节中看到）。
我们从研究单变量概念得知，固定时间序列通常会给我们一组更好的预测。 如果您不熟悉平稳性的概念，请首先阅读本文：处理非平稳时间序列的温度介绍
从数学上讲，如果时间序列的数据的均值和方差没有系统变化，且严格消除了周期变化，那么这个时间序列就是平稳的，只有平稳的时间序列数据的预测才是有可能的。非平稳时间序列数据目前在理论上有缺陷，预测很难。
预测大气质量，数据来源https://archive.ics.uci.edu/ml/datasets/Air+Quality#：
#import required packages
import pandas as pd
import matplotlib.pyplot as plt
%matplotlib inline
#read the data
df = pd.read_csv("AirQualityUCI.csv", parse_dates=[['Date', 'Time']])
#check the dtypes
df.dtypes
输出结果：
ate_Time object
CO(GT) int64
PT08.S1(CO) int64
NMHC(GT) int64
C6H6(GT) int64
PT08.S2(NMHC) int64
NOx(GT) int64
PT08.S3(NOx) int64
NO2(GT) int64
PT08.S4(NO2) int64
PT08.S5(O3) int64
T int64
RH int64
AH int64
dtype: object
Date_Time列的数据类型是object，我们需要将其更改为datetime。 此外，为了准备数据，我们需要索引具有日期时间。 请遵循以下命令：
df['Date_Time'] = pd.to_datetime(df.Date_Time , format = '%d/%m/%Y %H.%M.%S')
data = df.drop(['Date_Time'], axis=1)
data.index = df.Date_Time
下一步是处理缺失的值。 由于数据中的缺失值被替换为值-200，因此我们必须使用更好的数字来估算缺失值。 考虑一下 - 如果缺少当前的露点值，我们可以安全地假设它将接近前一个小时的值。 有道理，对吗？ 在这里，我将用前一个值来估算-200。
您可以选择使用之前几个值的平均值替换该值，或者选择前一天同时的值（您可以在下面的评论部分中分享您对缺失值进行估算的想法）。
#missing value treatment
cols = data.columns
for j in cols:
for i in range(0,len(data)):
if data[j][i] == -200:
data[j][i] = data[j][i-1]
#checking stationarity
from statsmodels.tsa.vector_ar.vecm import coint_johansen
#since the test works for only 12 variables, I have randomly dropped
#in the next iteration, I would drop another and check the eigenvalues
johan_test_temp = data.drop([ 'CO(GT)'], axis=1)
coint_johansen(johan_test_temp,-1,1).eig
结果：
array([ 0.17806667, 0.1552133 , 0.1274826 , 0.12277888, 0.09554265, 0.08383711, 0.07246919, 0.06337852, 0.04051374, 0.02652395, 0.01467492, 0.00051835])
我们现在可以继续创建验证集以适应模型，并测试模型的性能：
#creating the train and validation set
train = data[:int(0.8*(len(data)))]
valid = data[int(0.8*(len(data))):]
#fit the model
from statsmodels.tsa.vector_ar.var_model import VAR
model = VAR(endog=train)
model_fit = model.fit()
# make prediction on validation
prediction = model_fit.forecast(model_fit.y, steps=len(valid))
预测采用数组的形式，其中每个列表代表行的预测。 我们将把它转换成更易于呈现的格式。
#converting predictions to dataframe
pred = pd.DataFrame(index=range(0,len(prediction)),columns=[cols])
for j in range(0,13):
for i in range(0, len(prediction)):
pred.iloc[i][j] = prediction[i][j]
#check rmse
for i in cols:
print('rmse value for', i, 'is : ', sqrt(mean_squared_error(pred[i], valid[i])))
输出:
rmse value for CO(GT) is : 1.4200393103392812
rmse value for PT08.S1(CO) is : 303.3909208229375
rmse value for NMHC(GT) is : 204.0662895081472
rmse value for C6H6(GT) is : 28.153391799471244
rmse value for PT08.S2(NMHC) is : 6.538063846286176
rmse value for NOx(GT) is : 265.04913993413805
rmse value for PT08.S3(NOx) is : 250.7673347152554
rmse value for NO2(GT) is : 238.92642219826683
rmse value for PT08.S4(NO2) is : 247.50612831072633
rmse value for PT08.S5(O3) is : 392.3129907890131
rmse value for T is : 383.1344361254454
rmse value for RH is : 506.5847387424092
rmse value for AH is : 8.139735443605728
预测代码：
#make final predictions
model = VAR(endog=data)
model_fit = model.fit()
yhat = model_fit.forecast(model_fit.y, steps=1)
print(yhat)