多元时间序列数据的预测和建模

#### 简介

时间是决定企业是上升还是下降的最关键因素。 这就是为什么我们看到商店和电子商务平台的销售与节日保持一致。 这些企业分析了多年的支出数据，以了解打开大门并看到消费者支出增加的最佳时机。

但作为数据科学家，您如何进行此分析？ 别担心，你不需要建立时间机器！ 时间序列建模是一种强大的技术，可作为理解和预测趋势和模式的门户。

<center>![](https://s3-ap-south-1.amazonaws.com/av-blog-media/wp-content/uploads/2018/09/mts-850x478.jpg)</center>

但即使是时间序列模型也有不同的方面。 我们在网上看到的大多数例子都涉及单变量时间序列。 不幸的是，现实世界的用例并不像那样。 有多个变量在起作用，同时处理所有这些变量是数据科学家将获得其价值的地方。

在本文中，我们将了解多元变量时间序列是什么，以及如何处理它。 我们还将进行案例研究并在Python中实现它，以便您对该主题有实际的了解。

[TOC]

#### 一、单变量时间序列与多变量时间序列

本文假设您熟悉单变量时间序列，其属性以及用于预测的各种技术。 由于本文将重点关注多变量时间序列，我建议您阅读以下文章，作为单变量时间序列的一个很好的介绍：

- [创建时间序列预测的综合指南](https://www.analyticsvidhya.com/blog/2016/02/time-series-forecasting-codes-python/ "创建时间序列预测的综合指南")
- [使用Auto Arima构建高性能时间序列模型](https://www.analyticsvidhya.com/blog/2018/08/auto-arima-time-series-modeling-python-r/)

但在进入多变量时间序列的详细信息之前，我会给你快速回顾一下单变量时间序列。 让我们逐一看一下它们来理解它们之间的区别。

##### 1.1、单变量时间序列

顾名思义，单变量时间序列是一个具有单个时间相关变量的系列。

例如，查看下面的样本数据集，其中包含过去2年的温度值（每小时）。 这里，温度是因变量（取决于时间）。

<center>![](https://s3-ap-south-1.amazonaws.com/av-blog-media/wp-content/uploads/2018/09/var_3.png)</center>

如果要求我们预测接下来几天的温度，我们将查看过去的值并尝试测量和提取一种模式。 我们会注意到早上和晚上的温度较低，而下午达到峰值。 此外，如果您有过去几年的数据，您会发现它在11月到1月期间更冷，而在4月到6月则相对较热。

这些观察将有助于我们预测未来的价值观。 您是否注意到我们只使用了一个变量（过去2年的温度）？ 因此，这称为单变量时间序列分析/预测。

##### 1.2、多变量时间序列（MTS）

多变量时间序列具有多个时间相关变量。 每个变量不仅取决于其过去的值，而且还依赖于其他变量。 此依赖关系用于预测未来值。 听起来很复杂？ 让我解释。

考虑上面的例子。 现在假设我们的数据集包括出汗百分比、露点(dew point，在一定压力下出 现第一个液滴时的温度)、风速、云量百分比等，以及过去两年的温度值。 在这种情况下，需要考虑多个变量来最佳地预测温度。 像这样的系列将属于多变量时间序列的范畴。 以下是对此的说明：

<center>![](https://s3-ap-south-1.amazonaws.com/av-blog-media/wp-content/uploads/2018/09/var_4.png)</center>

现在我们了解了多变量时间序列的样子，让我们了解如何使用它来构建预测。

#### 二、处理多变量时间序列 - VAR
在本节中，我将向您介绍多变量时间序列预测中最常用的方法之一 ——即向量自回归模型（Vector Auto Regression，VAR）。

在VAR模型中，每个变量都是过去自身值和所有其他变量的过去值的线性函数。 为了更好地解释这一点，我将使用一个简单的可视化示例：

我们有两个变量，y1和y2。 我们需要在时间t预测这两个变量的值，从过去n值的给定数据开始。 为简单起见，我认为滞后值为1。

<center>![](https://s3-ap-south-1.amazonaws.com/av-blog-media/wp-content/uploads/2018/09/var_12.png)</center>

为了计算y1（t），我们将使用y1和y2的过去值。 类似地，为了计算y2（t），将使用y1和y2的过去值。 下面是一种表示这种关系的简单数学方法：

```math
y_1(t) = a_1 + w_{11} \times y_1(t-1) + w_{12} \times y_2(t-1) + e_1(t-1)
```

```math
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）过程：

```math
y(t) = a+w\times y(t-1) +e
```

在这个例子中，我们只有一个变量$y$，一个常数项$a$，一个误差项$e$和一个系数$w$。为了适应VAR的每个等式中的多个变量项，我们将使用向量。 我们可以用以下形式写出等式（1）和（2）：

<center>![](https://s3-ap-south-1.amazonaws.com/av-blog-media/wp-content/uploads/2018/09/vector_eqn1.png)</center>

上述形式扩展到多维之后：

<center>![](https://s3-ap-south-1.amazonaws.com/av-blog-media/wp-content/uploads/2018/09/vector_eqn2.png)</center>

```math
y(t) = a+ w_1\times y(t-1) + \cdots + w_p\times y(t-p) + \epsilon t
```

上述等式中的$\epsilon t$表示多变量白噪声。 对于多变量时间序列，$\epsilon t$应该是满足以下条件的连续随机向量：

- E(εt) = 0，即误差期望为0
- E(εt1,εt2‘) = σ12，εt和εt‘的期望是时间序列的标准差

#### 三、为什么需要VAR？
回想一下我们之前看到的温度预测的例子。 可以对其进行论证，将其视为多个单变量系列。 我们可以使用像AR这样的简单单变量预测方法来解决它。 由于目标是预测温度，我们可以简单地删除其他变量（温度除外）并在剩余的单变量系列上拟合模型。

另一个简单的想法是使用我们已知的技术单独预测每个系列的值。 这将使工作变得非常简单！ 那你为什么要学习另一种预测技术呢？ 这个话题不够复杂吗？

从上面的等式（1）和（2）可以清楚地看出，每个变量都使用每个变量的过去值来进行预测。 与AR不同，VAR能够理解和使用几个变量之间的关系。 这对于描述数据的动态行为非常有用，并且还可以提供更好的预测结果。 此外，实现VAR就像使用任何其他单变量技术一样简单（您将在上一节中看到）。

#### 四、多元时间序列的平稳性
我们从研究单变量概念得知，固定时间序列通常会给我们一组更好的预测。 如果您不熟悉平稳性的概念，请首先阅读本文：[处理非平稳时间序列的温度介绍](https://www.analyticsvidhya.com/blog/2018/09/non-stationary-time-series-python/)

从数学上讲，如果时间序列的数据的均值和方差没有系统变化，且严格消除了周期变化，那么这个时间序列就是平稳的，只有平稳的时间序列数据的预测才是有可能的。非平稳时间序列数据目前在理论上有缺陷，预测很难。

#### 五、Python的实现
预测大气质量，数据来源https://archive.ics.uci.edu/ml/datasets/Air+Quality#：

```python
#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
```
输出结果：
```python
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。 此外，为了准备数据，我们需要索引具有日期时间。 请遵循以下命令：

```python
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。

您可以选择使用之前几个值的平均值替换该值，或者选择前一天同时的值（您可以在下面的评论部分中分享您对缺失值进行估算的想法）。

```python
#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])
```

我们现在可以继续创建验证集以适应模型，并测试模型的性能：

```python
#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))
```

预测采用数组的形式，其中每个列表代表行的预测。 我们将把它转换成更易于呈现的格式。

```python
#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
```

预测代码：

```python
#make final predictions
model = VAR(endog=data)
model_fit = model.fit()
yhat = model_fit.forecast(model_fit.y, steps=1)
print(yhat)
```

本文翻译自:https://www.analyticsvidhya.com/blog/2018/09/multivariate-time-series-guide-forecasting-modeling-python-codes/

相关博客

最热博客