Transformer时间序列：PatchTST引领时间序列预测进

引言
- 为什么transformer框架可以应用到时间序列呢
- 统计学模型
- 深度学习模型
PatchTST
PatchTST模型架构原理。
- 通道独立性
- Patching
- patching的优点
- Transformer编码器
利用表示学习改进PatchTST
使用PatchTST模型进行预测
- 初始设置
- 使用transformer模型建模
- 预测
- 评估

在这里插入图片描述

引言

为什么transformer框架可以应用到时间序列呢

我们不妨看下传统的时间序列数据构造

在这里插入图片描述

这个过程与vision Transformer模型导入图像做patch处理很相似吧
在这里插入图片描述

大佬说不定就是看到这个点，才有了创新的🤔

为了更好的了解时间序列先回顾一下做时间序列预测的二种方法，一种是基于统计学模型，一种是基于深度学习

统计学模型

ARIMA(p,d,q)

差分自回归滑动平均（ARIMA）模型是自回归模型AR§和移动平均模型MA(q)的结合，但是针对的是经过差分的序列。

ARMA(p,d,q)，其中p是自回归部分的阶数，d是积分的阶数，q是移动平均部分的阶数。
可以用于非平稳数据。方程如下：
$y′_t=C+\phi_1y′_{t-1}+...\phi_py′_{t-p}+\theta_1\epsilon_{t-1}+...\theta_q\epsilon_{t-q}+\epsilon_t$

from statsmodels.tsa.statespace.sarimax import SARIMAX

model = SARIMAX(data,order=(p,d,q))
res = model.fit()
predictions = res.get_prediction(start_index,end_index)

注：积分的阶数d是指为使序列平稳而进行差分的次数。

SARIMA(p,d,q)(P,D,Q)m
季节性自回归移动平均（SARIMA）模型在ARIMA模型的基础上增加了季节性组成部分。

SARIMA(p,d,q)(P,D,Q)m。在这里，p,d和q与ARIMA模型中的含义相同。
P是自回归部分的季节性顺序
D是积分的季节性顺序
Q是移动平均部分的季节性顺序
m是数据的频率（即一个季节内的数据点数量）

from statsmodels.tsa.statespace.sarimax import SARIMAX

model = SARIMAX(data,order=(p,d,q),seasonal_order(P,D,Q,m))

res = model.fit()
predictions = res.get_prediction(start_index,end_index)

VARMAX

矢量自回归滑动平均模型与外生变量（VARMAX）模型用于多变量预测，即同时预测两个时间序列。

假设格兰杰因果性。必须使用格兰杰因果性检验。如果检验失败，则无法使用VARMAX模型。

from statsmodels.tsa.statespace.varmax import VARMAX

model = VARMAX(target,exog,order=(p,q))
res = model.fit(disp=False)
predictions = res.get_prediction(start_index,end_index)

BATS and TBATS

BATS和TBATS是在时间序列具有多个季节周期时使用的方法。这种情况通常发生在高频数据，比如每日数据的情况下。

当存在多个季节周期时，无法使用SARIMA模型。此时可以选择使用BATS或TBATS方法。
BATS方法包括以下组成部分：Box-Cox转换、ARMA误差、趋势和季节性成分。
TBATS方法包括以下组成部分：三角函数季节性、Box-Cox转换、ARMA误差、趋势和季节性成分。

from sktime.forecasting.bats import BATS
forecaster = BATS(sp=[period_1,period_2])
forecaster.fit(data)

predictions = forecaster.predict(horizon)

# Same syntax for TBATS
from sktime.forecasting.tbats import TBATS

forecaster = TBATS(sp=[period_1,period_2])
forecaster.fit(data)

predictions = forecaster.predict(horizon)

Exponential smoothing

指数平滑法是一种时间序列预测方法，常用于预测具有一定规律性的数据趋势。该方法假设未来的数值可以通过过去的数值来估计，并且随着时间的推移，过去的数值对于预测的影响逐渐减弱。

简单指数平滑法：简单指数平滑法是指数平滑法的一种基本形式，它只考虑过去的数值，并生成平稳的预测值。这意味着预测结果没有趋势性，只是一个平均水平的估计，用来生成平稳的预测值
双重指数平滑法：双重指数平滑法在简单指数平滑法的基础上增加了趋势分量的考虑。它假设未来的数据不仅受到过去的数值的影响，还受到趋势的影响。因此，预测结果是一条直线，可以是递增或递减的趋势。
三重指数平滑法：三重指数平滑法进一步考虑了季节性的影响。它假设未来的数值除了受到过去的数值和趋势的影响外，还受到季节性的周期性影响。这种方法可以用于处理具有季节性变化的数据，例如每年的销售数据。
趋势可以是“加法”或“指数形式”
季节性可以是“加法”或“乘法形式”

在指数平滑法中，趋势可以是“加法”或“指数形式”。加法趋势意味着趋势变化是线性的，而指数趋势则表示趋势变化是指数级别的。

同样地，季节性可以是“加法”或“乘法形式”。加法季节性表示季节性的影响是固定的，与整体水平无关。而乘法季节性则表示季节性的影响与整体水平成比例变化。

from statsmodels.tsa.holtwinters import ExponentialSmoothing

model = ExponentailSmoothing(data,trend='add',seasonal='add',seasoanl_periods=m,initialzation_method='estimated').fit()

predictions = model.forecast(horizon)

深度学习模型

DNN
深度神经网络（Deep Neural Network）由多个全连接层（Fully Connected Layers）堆叠而成，如果激活函数是非线性的，它可以对时间序列中的非线性关系进行建模。

首先从一个简单模型开始，只有几个隐藏层。在添加更多层之前，尝试增加训练时的迭代次数（epochs）进行实验。

from tensorflow.keras import Sequential
from tensorflow.keras.layers import Dense

dnn = Sequential([Dense(units=64,activation='relu'),
				  Dense(units=64,activation='relu'),
				  Dense(units=1)])

dnn.compile(loss='mean_squared_error',
			optimizer='adam',
			metrics=['mean_absolute_error'])

LSTM
LSTM（长短期记忆网络）在处理文本和时间序列等数据序列方面表现出色。它的架构允许过去的信息仍然被用于后续的预测。

您可以在模型中堆叠多个LSTM层。
您可以尝试将LSTM与CNN结合使用。
由于数据按序列处理，因此训练LSTM需要更长的时间

from tensorflow.keras import Sequential
from tensorflow.keras.layers import Dense,LSTM

lstm_model = Sequential([LSTM(32,return_sequences=True),
						 Dense(units1=1)])

lstm_model.compile(loss='mean_squared_error',
				   optimizer='adam',
				   metrics=['mean_absolute_error'])

CNN
卷积神经网络（CNN）可以作为我们时间序列的过滤器，这是由于卷积运算的特性，它可以减小特征空间。

CNN的训练速度比LSTM更快。
可以与LSTM结合使用。将CNN层放置在LSTM之前

from tensorflow.keras import Sequential
from tensorflow.keras.layers import Dense,Conv1D

cnn_model = Sequential([
			Conv1D(filter=32,
				   kernel_size=(KERNEL_WIDTH,),
				   activation='relu'),
			Dense(units=32,activation='relu'),
			Dense(units=1)])

cnn_model.compile(loss='mean_squared_error',
				  optimizer='adam',
				  metrics=['mean_absolute_error'])

🎙🕺🎵🏀现在回到说transformer的时间序列模型

基于Transformer的模型已经成功应用于许多领域，比如自然语言处理（如BERT或GPT模型）和计算机视觉等。

然而，在时间序列方面，最先进的结果主要是由MLP模型（多层感知器）如N-BEATS和N-HiTS实现的。最近的一篇论文甚至表明，简单的线性模型在许多基准数据集上优于复杂的基于Transformer的预测模型（参见Zheng等人，2022）。

然而，已经提出了一种新的基于Transformer的模型，该模型在长期预测任务中实现了最先进的结果：PatchTST。

PatchTST代表patch时间序列Transformer，它首次提出于2023年3月，由Nie、Nguyen等人在他们的论文《A TIME SERIES IS WORTH 64 WORDS:LONG-TERM FORECASTING WITH TRANSFORMERS》中。与其他基于Transformer的模型相比，他们提出的方法取得了最先进的结果。

在本文中，我们首先通过直观的方式探讨PatchTST的内部工作原理，不使用公式。然后，我们将该模型应用于一个预测项目，并将其性能与MLP模型（如N-BEATS和N-HiTS）进行比较和评估。

论文标题：A TIME SERIES IS WORTH 64 WORDS:
LONG-TERM FORECASTING WITH TRANSFORMERS
论文地址：https://arxiv.org/pdf/2211.14730.pdf

PatchTST

PatchTST：顾名思义，它利用了 Patch 和 Transformer 架构。它还包括通道独立性，用于处理多变量时间序列。下面是该模型的总体架构。
在这里插入图片描述

在这里，关键要点 PatchTST 使用通道独立性来预测多变量时间序列。然后，在其 Transformer 主干中，模型使用 Patch，这些 Patch 由小的垂直矩形表示。此外，该模型有两个版本：有监督和自监督。

PatchTST模型架构原理。

通道独立性

在这里，多变量时间序列被视为多通道信号。每个时间序列基本上是包含信号的通道。也就是说有多少通道其实就是有多少条时间序列。

在上图中，我们可以看到多变量时间序列通过（Instance Norm + Patching）分隔为各个序列，并将每个序列作为输入token提供给 Transformer 主干。然后，对每个序列进行预测，并将结果连接以得到最终的预测结果。

Patching

大多数基于 Transformer 的预测模型的研究工作都集中在构建新的机制，以简化原始的注意力机制。然而，它们仍然依赖于逐个点计算注意力，这在处理时间序列时并不理想。

在时间序列预测中，我们希望提取过去时间步和未来时间步之间的关系，以进行预测。使用逐点注意力时，我们尝试从单个时间步中检索信息，而不考虑该点周围的内容。

换句话说，我们计算的注意力机制是隔一段距离取一个点来计算注意力机制，不查看之前或之后的点。

为什么时间序列计算注意力机制不是逐个点计算呢？

时间序列数据具有时序性，各个时间步之间存在着相关性和依赖关系。在时间序列预测任务中，我们需要从过去的时间步中提取信息，以预测未来的时间步。而传统的逐点计算注意力机制只关注当前时间步，无法充分考虑到时间序列中前后时间步的联系。

举个例子

假设我们要预测某股票的未来价格。在时间序列中，前一天的价格可能对接下来几天的价格走势有一定的影响。如果仅仅使用逐点计算的注意力机制，模型只能关注当前时间步的价格，而无法获取到前一天的价格信息。这样会导致模型无法捕捉到时间序列中的长期依赖性和趋势，从而影响预测的准确性。

这就像在不看句子中单词周围的词汇的情况下理解一个单词的意义。

因此，PatchTST 利用 Patch 来提取时间序列中的局部语义信息。

patching 的工作原理

每个输入序列被划分为 Patch，它们只是来自原始序列的较短序列。
在这里插入图片描述

在这里，Patch 之间可以是重叠的或非重叠的。

在这里插入图片描述

Patch 的数量取决于 Patch 的长度 P 和步长 S。这里的步长类似于卷积，它只是连续 Patch 的起始时间步之间的时间步数差。

在这里插入图片描述
在上图中，我们可以直观地看到patch的结果。这里，我们有一个15个时间步长的序列长度（L），patch长度（P）为5，步长（S）为5。结果是该系列被分离为3个patch。

patching的优点

通过patching，模型可以通过观察一组时间步骤而不是单个时间步骤来提取局部语义含义。

这还有一个额外的好处，就是大大减少了馈送到Transformer编码器的标记数量。在这里，每个patch都变成了输入到Transformer的一个token。这样，我们可以将token的数量从L减少到大约L/S。

大大降低了模型的空间和时间复杂度。这反过来意味着我们可以向模型输入更长的输入序列以提取有意义的时间关系。

Transformer编码器

一旦对序列进行了path，就将其馈送到Transformer编码器中。这是经典的Transformer架构，没有进行修改。

然后，输出被馈送到线性层，并进行预测。

利用表示学习改进PatchTST

论文的作者还提出了通过使用表示学习来改进模型的方法。

在这里插入图片描述

从上面的图中，我们可以看到PatchTST可以使用自监督的表示学习来捕捉数据的抽象表示。这可能会导致预测性能的潜在改进。

改进过程相当简单，随机的将patch进行mask屏蔽（bert的操作），意味着它们将被设为0。如上图中的空白垂直矩形所示。然后，模型被训练以重新创建原始的patch，这是图上方灰色垂直矩形所表示的。

使用PatchTST模型进行预测

在论文中，PatchTST与其他基于Transformer的模型进行了比较。包括之前发布了一些基于多层感知器（MLP）的模型，如N-BEATS和N-HiTS

结果效果显示PatchTST在长期预测任务上也表现出了最先进的性能。

它的源代码在github上有的：https://github.com/marcopeix/datasciencewithmarco/blob/master/PatchTST.ipynb

现在我们来看下这个项目代码吧

在这个实例中，他们使用了Exchange数据集，这是一个在研究中用于长期预测的常见基准数据集。该数据集包含了从1990年到2016年的八个国家相对于美元的每日汇率。该数据集通过MIT许可证提供。

初始设置

导入所需的库。在这里，将使用neuralforecast，因为它们有一个即插即用的PatchTST实现。对于数据集，使用了datasetsforecast库，其中包括用于评估预测算法的所有流行数据集。

import torch
import numpy as np
import pandas as pd
import matplotlib.pyplot as plt

from neuralforecast.core import NeuralForecast
from neuralforecast.models import NHITS, NBEATS, PatchTST

from neuralforecast.losses.pytorch import MAE
from neuralforecast.losses.numpy import mae, mse

from datasetsforecast.long_horizon import LongHorizon

下载Exchange数据集。

Y_df, X_df, S_df = LongHorizon.load(directory="./data", group="Exchange")

在这里，我们获得了三个数据框。第一个数据框包含每个国家的每日汇率数据。第二个数据框包含外部时间序列数据。第三个数据框包含静态外部变量，例如日期、月份、年份、小时或我们已知的其他未来信息。

在这个实例中，我们只使用Y_df。

Y_df['ds'] = pd.to_datetime(Y_df['ds'])

Y_df.head()

在这里插入图片描述

从上面的图中，我们可以看到我们有三列。第一列是一个唯一的标识符，当使用neuralforecast时，拥有一个id列是必要的。然后，ds列包含日期，y列包含汇率。

Y_df['unique_id'].value_counts()

在这里插入图片描述

从上面的图片中，我们可以看到每个唯一的id对应一个国家，并且每个国家有7588个观测值。

现在，我们定义验证集和测试集的大小。在这里，我选择了760个时间步长作为验证集，1517个时间步长作为测试集，这是根据datasets库的规定。

val_size = 760
test_size = 1517

print(n_time, val_size, test_size)

然后，让我们绘制其中一个序列，以了解我们正在处理的数据。在这里，绘制第一个国家（unique_id = 0）的序列.

u_id = '0'

x_plot = pd.to_datetime(Y_df[Y_df.unique_id==u_id].ds)
y_plot = Y_df[Y_df.unique_id==u_id].y.values

x_plot

x_val = x_plot[n_time - val_size - test_size]
x_test = x_plot[n_time - test_size]

fig, ax = plt.subplots(figsize=(12,8))
ax.plot(x_plot, y_plot)
ax.set_xlabel('Date')
ax.set_ylabel('Exhange rate')
ax.axvline(x_val, color='black', linestyle='--')
ax.axvline(x_test, color='black', linestyle='--')

plt.text(x_val, -2, 'Validation', fontsize=12)
plt.text(x_test,-2, 'Test', fontsize=12)

plt.tight_layout()

在这里插入图片描述
从上图可以看出，数据非常嘈杂，没有明显的季节性。

使用transformer模型建模

开始使用neuralforecast进行建模。

首先，我们需要设置时间跨度。在这种情况下，我使用了96个时间步长，因为这个时间跨度也在PatchTST论文中使用过。

然后，为了对每个模型进行公正评估，我决定将输入大小设置为时间跨度的两倍（即192个时间步长），并将最大迭代次数设置为50。其他超参数保持默认值。

horizon = 96

models = [NHITS(h=horizon,
               input_size=2*horizon,
               max_steps=50),
         NBEATS(h=horizon,
               input_size=2*horizon,
               max_steps=50),
         PatchTST(h=horizon,
                 input_size=2*horizon,
                 max_steps=50)

接下来，我们通过指定要使用的模型和预测频率（在本例中是每日）来初始化NeuralForecast对象。

nf = NeuralForecast(models=models, freq='D')

预测

为了生成预测结果，我们使用交叉验证方法利用验证集和测试集。它将返回一个包含所有模型预测和相关真实值的数据框。

preds_df = nf.cross_validation(df=Y_df, val_size=val_size, test_size=test_size, n_windows=None)

在这里插入图片描述
为了评估模型，我们将预测的值与真实值都放在时间维度上进行对比

fig, ax = plt.subplots(figsize=(12,8))

ax.plot(y_true[0, 0, :], label='True')
ax.plot(y_pred_nhits[0, 0, :], label='N-HiTS', ls='--')
ax.plot(y_pred_nbeats[0, 0, :], label='N-BEATS', ls=':')
ax.plot(y_pred_patchtst[0, 0, :], label='PatchTST', ls='-.')
ax.set_ylabel('Exchange rate')
ax.set_xlabel('Forecast horizon')
ax.legend(loc='best')

plt.tight_layout()

在这里插入图片描述
从上图看N-BEATS和N-HiTS的预测结果与实际值相差很大。然而，PatchTST虽然也有误差，但似乎是最接近实际值的。

评估

因此，让我们评估每个模型的性能。为了复制论文中的方法，我们使用平均绝对误差（MAE）和均方误差（MSE）作为性能指标。

data = {'N-HiTS': [mae(y_pred_nhits, y_true), mse(y_pred_nhits, y_true)],
       'N-BEATS': [mae(y_pred_nbeats, y_true), mse(y_pred_nbeats, y_true)],
       'PatchTST': [mae(y_pred_patchtst, y_true), mse(y_pred_patchtst, y_true)]}

metrics_df = pd.DataFrame(data=data)
metrics_df.index = ['mae', 'mse']

metrics_df.style.highlight_min(color='lightgreen', axis=1)