前言

系列专栏:【深度学习：算法项目实战】✨︎
涉及医疗健康、财经金融、商业零售、食品饮料、运动健身、交通运输、环境科学、社交媒体以及文本和图像处理等诸多领域，讨论了各种复杂的深度神经网络思想，如卷积神经网络、循环神经网络、生成对抗网络、门控循环单元、长短期记忆、自然语言处理、深度强化学习、大型语言模型和迁移学习。

在数据科学领域中，时间序列预测一直是一个重要且挑战性的任务。随着深度学习技术的兴起，特别是循环神经网络（RNN）及其变种如长短期记忆网络（LSTM）和门控循环单元（GRU）的引入，时间序列预测的性能得到了显著提升。在这些模型中，双向长短期记忆网络（BiLSTM）因其能够同时捕获序列中的前向和后向信息，而在许多应用中表现出色。

在本文中，我们将探讨如何使用PyTorch框架实现基于双向长短期记忆网络（BiLSTM）的多变量时间序列预测模型。多变量时间序列预测涉及多个输入变量和一个或多个输出变量，这些变量随时间变化并可能相互影响。通过利用BiLSTM的双向特性，我们的模型将能够同时考虑历史数据和未来趋势，从而更准确地预测时间序列的未来值。

1. 股票时间序列数据

股票时间序列数据是指按照时间顺序排列的关于股票市场的历史数据。这些数据记录了股票市场在不同时间点上的各种信息，如股票的开盘价、最高价、最低价、收盘价、成交量等。时间序列数据是金融市场分析的重要基础，因为它们反映了市场参与者的行为、市场供求关系以及宏观经济和政策的影响等因素。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as snsfrom sklearn.preprocessing import MinMaxScaler
from sklearn.model_selection import TimeSeriesSplitimport torch
import torch.nn as nn
from torch.utils.data import DataLoader, TensorDataset
from torchinfo import summary
from tqdm import tqdm

np.random.seed(0)
torch.manual_seed(0)

1.1 数据预处理

pandas.to_datetime 函数将标量、数组、Series 或 DataFrame/dict-like 转换为 pandas datetime 对象。

AAPL = pd.read_csv('AAPL.csv')
print(type(AAPL['Close'].iloc[0]),type(AAPL['Date'].iloc[0]))
# Let's convert the data type of timestamp column to datatime format
AAPL['Date'] = pd.to_datetime(AAPL['Date'])
print(type(AAPL['Close'].iloc[0]),type(AAPL['Date'].iloc[0]))# Selecting subset
cond_1 = AAPL['Date'] >= '2021-04-23 00:00:00'
cond_2 = AAPL['Date'] <= '2024-04-23 00:00:00'
AAPL = AAPL[cond_1 & cond_2].set_index('Date')
print(AAPL.shape)

<class 'numpy.float64'> <class 'str'>
<class 'numpy.float64'> <class 'pandas._libs.tslibs.timestamps.Timestamp'>
(755, 6)

---

1.2 数据可视化

收盘价是股票在正常交易日的最后交易价格。股票的收盘价是投资者跟踪其长期表现的标准基准。

# plt.style.available
plt.style.use('_mpl-gallery')
plt.figure(figsize=(18,6))
plt.title('Close Price History')
plt.plot(AAPL['Close'],label='AAPL')
plt.ylabel('Close Price USD ($)', fontsize=18)
plt.legend()
plt.show()

2. 时间数据特征工程(APPL)

在时间序列分析中，时间窗口通常用于描述在训练模型时考虑的连续时间步 time steps 的数量。这个时间窗口的大小，即 window_size，对于模型预测的准确性至关重要。

具体来说，window_size 决定了模型在做出预测时所使用的历史数据的长度。例如，如果我们想要用前60天的股票数据来预测未来7天的收盘价，那么window_size 就是60。

# 设置时间窗口大小
window_size = 60

除此之外还需构建序列函数，该函数需要两个参数：dataset 和 lookback，前者是要转换成数据集的 NumPy 数组，后者是用作预测下一个时间段的输入变量的前一时间步数，默认设为 1。

# 构造序列数据函数，若target为 Close，即索引为 3
def create_dataset(dataset, lookback= 1):X, y = [], []for i in range(len(dataset)-lookback): feature = dataset[i:(i+lookback), :]target = dataset[i + lookback, 3]X.append(feature)y.append(target)return np.array(X), np.array(y)

2.1 特征缩放（归一化）

MinMaxScaler() 函数主要用于将特征数据按比例缩放到指定的范围。默认情况下，它将数据缩放到[0, 1]区间内，但也可以通过参数设置将数据缩放到其他范围。在机器学习中，MinMaxScaler()函数常用于不同尺度特征数据的标准化，以提高模型的泛化能力。

# 选取 AAPL[['Open', 'High', 'Low', 'Close']]作为特征, 归一化数据
scaler = MinMaxScaler(feature_range=(0, 1))
scaled_data = scaler.fit_transform(AAPL[['Open', 'High', 'Low', 'Close']].values)# 获取反归一化参数(即原始数据的最小值和最大值)
original_min = scaler.data_min_
original_max = scaler.data_max_# scale_params是一个包含所有特征反归一化参数的列表或数组，
# 其中第四个元素是Close价格的反归一化参数
scale_params = original_max - original_min

2.2 数据集划分（TimeSeriesSplit）

TimeSeriesSplit() 函数与传统的交叉验证方法不同，TimeSeriesSplit 特别适用于需要考虑时间顺序的数据集，因为它确保测试集中的所有数据点都在训练集数据点之后，并且可以分割多个训练集和测试集。

# 创建数据集,数据形状为 [samples, time steps, features]
X, y = create_dataset(scaled_data, lookback = window_size)
print(X.shape, y.shape)

(695, 60, 4) (695,)

# 使用TimeSeriesSplit划分数据集，根据需要调整n_splits
tscv = TimeSeriesSplit(n_splits=3, test_size=90)
# 遍历所有划分进行交叉验证
for i, (train_index, test_index) in enumerate(tscv.split(X)):X_train, X_test = X[train_index], X[test_index]y_train, y_test = y[train_index], y[test_index]# print(f"Fold {i}:")# print(f"  Train: index={train_index}")# print(f"  Test:  index={test_index}")# 查看最后一个 fold 数据帧的维度
print(X_train.shape, X_test.shape, y_train.shape, y_test.shape)

(605, 60, 4) (90, 60, 4) (605,) (90,)

2.3 数据集张量（TensorDataset）

张量是一个多维数组或矩阵的数学对象，可以看作是向量和矩阵的推广。在深度学习中，张量通常用于表示输入数据、模型参数以及输出数据

# 将 NumPy数组转换为 tensor张量
X_train_tensor = torch.from_numpy(X_train).type(torch.Tensor)
X_test_tensor = torch.from_numpy(X_test).type(torch.Tensor)
y_train_tensor = torch.from_numpy(y_train).type(torch.Tensor).view(-1, 1)
y_test_tensor = torch.from_numpy(y_test).type(torch.Tensor).view(-1, 1)print(X_train_tensor.shape, X_test_tensor.shape, y_train_tensor.shape, y_test_tensor.shape)

view() 函数用于重塑张量对象，它等同于 NumPy 中的 reshape() 函数，允许我们重组数据，以匹配 BiLSTM 模型所需的输入形状。以这种方式重塑数据可确保 BiLSTM 模型以预期格式接收数据。

torch.Size([605, 60, 4]) torch.Size([90, 60, 4]) torch.Size([605, 1]) torch.Size([90, 1])

接下来，我们将使用 TensorDataset 和 DataLoader创建数据集和数据加载器

train_dataset = TensorDataset(X_train_tensor, y_train_tensor)
train_loader = DataLoader(train_dataset, batch_size=32, shuffle=True)
test_dataset = TensorDataset(X_test_tensor, y_test_tensor)
test_loader = DataLoader(test_dataset, batch_size=32, shuffle=False)

shuffle=True 表示在每个epoch开始时，数据集将被随机打乱，这有助于防止模型在训练时过拟合。与训练数据加载器类似，shuffle=False 表示在测试时不需要打乱数据集。因为测试集通常用于评估模型的性能，而不是用于训练，所以不需要打乱。

3. 构建时间序列模型（BiLSTM）

双向LSTM（Bi-directional Long Short-Term Memory）网络可以通过使用 nn.LSTM 类并设置 bidirectional=True 参数来完成

PyTorch所提供的数学公式及解释如下：🔗

---

Apply a multi-layer long short-term memory (LSTM) RNN to an input sequence. For each element in the input sequence, each layer computes the following function:
$$\begin{array}{l}\\ i_t=\sigma (W_{ii}{x}_t+b_{ii}+W_{hi}{h}_{t-1}+b_{hi})\\ f_t=\sigma (W_{if}{x}_t+b_{if}+W_{hf}{h}_{t-1}+b_{hf})\\ g_t=\tanh (W_{ig}{x}_t+b_{ig}+W_{hg}{h}_{t-1}+b_{hg})\\ o_t=\sigma (W_{io}{x}_t+b_{io}+W_{ho}{h}_{t-1}+b_{ho})\\ c_t=f_t\odot c_{t-1}+i_t\odot g_t\\ h_t=o_t\odot \tanh (c_t)\end{array}$$
where $h_t$ is the hidden state at time $t$, $c_t$ is the cell state at time $t$, $x_t$ is the input at time $t$, $h_{t-1}$ is the hidden state of the layer at time $t-1$ or the initial hidden state at time $0$, and $i_t$, $f_t$, $g_t$, $o_t$ are the input, forget, cell, and output gates, respectively. $\sigma$ is the sigmoid function, and $\odot$ is the Hadamard product.

In a multilayer LSTM, the input $x_t^{(l)}$ of the $l$ -th layer ($l\ge 2$) is the hidden state $h_t^{(l-1)}$ of the previous layer multiplied by dropout ${\delta }_t^{(l-1)}$ where each ${\delta }_t^{(l-1)}$ is a Bernoulli random variable which is $0$ with probability $dropout$.

If proj_size > 0 is specified, LSTM with projections will be used. This changes the LSTM cell in the following way. First, the dimension of $h_t$ will be changed from hidden_size to proj_size (dimensions of $W_{hi}$ will be changed accordingly). Second, the output hidden state of each layer will be multiplied by a learnable projection matrix: $h_t=W_{hr}{h}_t$. Note that as a consequence of this, the output of LSTM network will be of different shape as well. See Inputs/Outputs sections below for exact dimensions of all variables. You can find more details in https://arxiv.org/abs/1402.1128.

---

3.1 构建Bi-LSTM 模型

class BiLSTM(nn.Module):def __init__(self, input_size, hidden_size, num_layers, output_size):# input_size 是输入特征的维度，hidden_size 是LSTM隐藏层神经单元维度（或称为隐藏状态的大小），# num_layers 是LSTM网络层数，output_size 是输出维度super(BiLSTM, self).__init__()# 通过调用 super(BiLSTM, self).__init__() 初始化父类 nn.Moduleself.hidden_size = hidden_sizeself.num_layers = num_layersself.lstm = nn.LSTM(input_size, hidden_size, num_layers, batch_first=True, bidirectional=True)# 定义 LSTM 层，使用 batch_first=True 表示输入数据的形状是 [batch_size, seq_len(time_steps), input_size]# 使用 bidirectional=True 表示的是 bidirectional LSTM，即双向 LSTMself.fc = nn.Linear(hidden_size * 2, output_size)  # 定义全连接层，将 LSTM 的最后一个隐藏状态映射到输出维度 output_size，双向LSTM的输出维度是隐藏层维度的两倍def forward(self, x):# 初始化隐藏状态(hidden state)和单元状态(cell state)为全零张量h0 = torch.zeros(self.num_layers * 2, x.size(0), self.hidden_size).to(x.device)  # 双向LSTM，所以是2倍c0 = torch.zeros(self.num_layers * 2, x.size(0), self.hidden_size).to(x.device)  # LSTM输出, (output, (hn, cn))out, _ = self.lstm(x, (h0, c0))  # out: tensor of shape (batch_size, seq_length(time_steps), hidden_size*2)# 取最后一个时间步的输出（或者对输出进行平均/最大池化等操作）# 这里我们选择最后一个时间步的输出out = self.fc(out[:, -1, :])return out

在时间序列模型中 seq_len 代表的就是时间步 time_steps
batch_first=True 表示输入数据的形状是 [batch_size, seq_len, input_size]，输出数据形状是 [batch_size, seq_len, num_directions * hidden_size]
batch_first=False 表示输入数据的形状是 [seq_len, batch_size, input_size]，输出数据形状是 [seq_len, batch_size, num_directions * hidden_size]

3.2 定义模型、损失函数与优化器

在神经网络中，模型通常指的是神经网络的架构，它定义了输入数据如何通过网络进行转换和计算，从而得到输出。在本案例中我们可以通过调用 BiLSTM 类来定义 Bi-directional Long Short-Term Memory。

损失函数（又称代价函数）是衡量模型预测值与真实值之间不一致性的指标。在训练神经网络时，我们的目标是找到一组网络参数，它使得损失函数的值最小化。损失函数的选择取决于问题的类型，例如，对于回归问题，常用的损失函数包括均方误差（MSE）；对于分类问题，则常使用交叉熵损失等。在PyTorch中，可以使用torch.nn模块中的损失函数类，如nn.MSELoss用于回归问题，nn.CrossEntropyLoss用于分类问题

优化器的任务是通过更新网络的权重和偏置来最小化损失函数。这通常是通过反向传播算法完成的，该算法计算损失相对于每个参数的梯度，并使用梯度下降或其变体来更新参数。在PyTorch中，torch.optim模块提供了多种优化器实现，如随机梯度下降（SGD）、带动量的SGD（Momentum）、RMSProp、Adam等。这些优化器都是基于梯度下降算法进行改进的，旨在更有效地更新模型参数。

model = BiLSTM(input_size = X_train_tensor.size(2), # 输入数据的特征数量 X_train.shape[2]hidden_size = 64,num_layers = 3,output_size = 1)
criterion = torch.nn.MSELoss() # 定义均方误差损失函数
optimizer = torch.optim.Adam(model.parameters(), lr=0.01) # 定义优化器

summary(model, (32, 60, 4)) # batch_size, seq_len(time_steps), input_size

==========================================================================================
Layer (type:depth-idx)                   Output Shape              Param #
==========================================================================================
BiLSTM                                   [32, 1]                   --
├─LSTM: 1-1                              [32, 60, 128]             234,496
├─Linear: 1-2                            [32, 1]                   129
==========================================================================================
Total params: 234,625
Trainable params: 234,625
Non-trainable params: 0
Total mult-adds (Units.MEGABYTES): 450.24
==========================================================================================
Input size (MB): 0.03
Forward/backward pass size (MB): 1.97
Params size (MB): 0.94
Estimated Total Size (MB): 2.94
==========================================================================================

4. 模型训练与可视化

4.1 模型训练与权重更新（保存最佳模型）

train_losses = [] # 初始化列表来存储损失值
test_losses = [] # 初始化列表来存储损失值
best_loss = float('100') # 初始化为一个很大的值
best_model = None # 初始化最佳模型为None# 训练循环
num_epochs = 20  # 假设我们要训练 20个 epoch
for epoch in range(num_epochs):model.train()  # 初始化训练进程train_loss = 0.0  # 初始化每个epoch的损失总和pbar = tqdm(train_loader, desc=f"Epoch {epoch+1}/{num_epochs}")for batch_idx, (data, target) in enumerate(pbar):optimizer.zero_grad()  # 将优化器中所有参数的梯度归零output = model(data)  # 每个批次的预测值loss = criterion(output, target)  # 计算模型损失值(y_train_pred 与 y_train)# 累加每个批次的损失train_loss += loss.item() * data.size(0)# 反向传播和优化loss.backward()optimizer.step()# 更新进度条(显示当前批次的损失)pbar.update()# 计算当前epoch的平均损失average_loss = train_loss / len(train_loader.dataset)train_losses.append(average_loss)model.eval()  # 将模型设置为评估模式with torch.no_grad():  # 不计算梯度以节省内存和计算资源test_loss = 0.0for data, target in test_loader:output = model(data) # 每个批次的预测值loss = criterion(output, target)test_loss += loss.item() * data.size(0)# 计算测试集的平均损失test_loss = test_loss / len(test_loader.dataset)test_losses.append(test_loss)# 如果测试损失比之前的最小值小，则更新最佳模型权重if test_loss < best_loss:best_loss = test_lossbest_model = model.state_dict() # 训练循环结束后保存最佳模型
if best_model is not None:torch.save(best_model, 'model_best.pth')print(f'Best model saved with test loss {best_loss:.4f}')
else:print('No model saved as no improvement was made.')

loss.item()：loss是一个PyTorch的Tensor对象，而.item()方法用于从只包含一个元素的Tensor中提取Python数值（通常是浮点数）。由于损失通常是一个标量（scalar），因此我们可以安全地使用.item()来获取其数值。
data.size(0)：这返回data这个Tensor在第一维度（通常是批次维度）的大小。在深度学习中，我们经常使用mini-batch梯度下降，这意味着我们将整个数据集分成较小的批次，并在每个批次上计算梯度并更新模型。data.size(0)就告诉我们当前批次中有多少个数据样本。
loss.item() * data.size(0)：这实际上是在对损失进行“加权”。我们乘以批次大小是为了稍后能够计算平均损失（average loss）。这是因为在累积损失时，我们不仅要考虑损失的大小，还要考虑该损失是基于多少数据样本计算出来的。
train_loss += ...：最后，我们将加权后的损失累加到train_loss变量中。这样做是为了在整个训练循环结束时计算平均损失。

Epoch 1/20: 100%|███████████████████████████████| 19/19 [00:01<00:00, 13.24it/s]
Epoch 2/20: 100%|███████████████████████████████| 19/19 [00:01<00:00, 14.14it/s]
Epoch 3/20: 100%|███████████████████████████████| 19/19 [00:01<00:00, 13.10it/s]
Epoch 4/20: 100%|███████████████████████████████| 19/19 [00:01<00:00, 14.17it/s]
Epoch 5/20: 100%|███████████████████████████████| 19/19 [00:01<00:00, 14.00it/s]
Epoch 6/20: 100%|███████████████████████████████| 19/19 [00:01<00:00, 14.26it/s]
Epoch 7/20: 100%|███████████████████████████████| 19/19 [00:01<00:00, 14.48it/s]
Epoch 8/20: 100%|███████████████████████████████| 19/19 [00:01<00:00, 13.82it/s]
Epoch 9/20: 100%|███████████████████████████████| 19/19 [00:01<00:00, 13.80it/s]
Epoch 10/20: 100%|██████████████████████████████| 19/19 [00:01<00:00, 14.18it/s]
Epoch 11/20: 100%|██████████████████████████████| 19/19 [00:01<00:00, 14.16it/s]
Epoch 12/20: 100%|██████████████████████████████| 19/19 [00:01<00:00, 14.30it/s]
Epoch 13/20: 100%|██████████████████████████████| 19/19 [00:01<00:00, 13.36it/s]
Epoch 14/20: 100%|██████████████████████████████| 19/19 [00:01<00:00, 12.76it/s]
Epoch 15/20: 100%|██████████████████████████████| 19/19 [00:01<00:00, 15.39it/s]
Epoch 16/20: 100%|██████████████████████████████| 19/19 [00:01<00:00, 14.80it/s]
Epoch 17/20: 100%|██████████████████████████████| 19/19 [00:01<00:00, 13.82it/s]
Epoch 18/20: 100%|██████████████████████████████| 19/19 [00:01<00:00, 12.27it/s]
Epoch 19/20: 100%|██████████████████████████████| 19/19 [00:01<00:00, 14.48it/s]
Epoch 20/20: 100%|██████████████████████████████| 19/19 [00:01<00:00, 13.97it/s]
Best model saved with test loss 0.0010

best_epoch_idx = test_losses.index(min(test_losses))
print(f"Best epoch: {best_epoch_idx+1}, with test loss: {min(test_losses):.4f}")

Best epoch: 20, with test loss: 0.0010

4.2 可视化训练过程（Loss损失）

# 绘制损失图
plt.figure(figsize=(18, 6))
plt.plot(train_losses, label='Train Loss')
plt.plot(test_losses, label='Test Loss')
plt.title('Training and Test Loss')
plt.xlabel('Epoch')
plt.ylabel('Loss')
plt.legend()
plt.show()

5. 模型评估与可视化

5.1 评估指标（MAE、RMSE、MAPE）

RMSE：该指标计算预测值与实际值之间平均平方差的平方根。它对较大误差给予较高的惩罚。
MAPE：该指标计算实际值与预测值之间绝对百分比差异的平均值。它用百分比来表示平均绝对误差，这对了解相对预测误差非常有用。

# 加载最佳模型
model.load_state_dict(torch.load('model_best.pth'))

def mean_absolute_percentage_error(y_true, y_pred):y_true, y_pred = np.array(y_true), np.array(y_pred)mask = y_true != 0return np.mean(np.abs((y_true[mask] - y_pred[mask]) / y_true[mask])) * 100 if mask.any() else 0# 将模型设置为评估模式
model.eval()
y_pred_all = []
y_true_all = []with torch.no_grad():pbar = tqdm(test_loader, desc='Evaluating')for data, target in pbar:y_pred = model(data).detach().cpu().numpy()  # 确保张量在CPU上y_true = target.detach().cpu().numpy()  # 确保张量在CPU上y_pred_all.append(y_pred)y_true_all.append(y_true)pbar.update()# 合并所有批次的预测和真实值
y_pred_all = np.concatenate(y_pred_all)
y_true_all = np.concatenate(y_true_all)mae = np.mean(np.abs(y_pred_all - y_true_all))
print(f"MAE: {mae:.4f}")rmse = np.sqrt(np.mean((y_pred_all - y_true_all) ** 2))
print(f"RMSE: {rmse:.4f}")mape = mean_absolute_percentage_error(y_true_all, y_pred_all)
print(f"MAPE: {mape:.4f}%")

Evaluating: 100%|████████████████████████████| 3/3 [00:00<00:00, 63.98it/s]
MAE: 0.0241
RMSE: 0.0321
MAPE: 3.2193%

5.2 反归一化

.inverse_transform 函数将经过转换或缩放的数据转换回其原始形式或接近原始形式，若在这里使用会发生错误，是因为用一个形状为(605,4)的数组去广播（broadcast）到一个形状为(605, 1)的数组上，这两个形状不兼容，它们的第二个维度不匹配，这里我们将使用下面的函数进行反归一化

# 提取Close价格的反归一化参数
close_denorm_param = scale_params[3]  # 假设Close是第四个特征# 如果反归一化包括一个偏移量shift，需要再加上 shift_param
# denormalized_predictions = (predictions[:, 0] * scale_param) + shift_param

# 将训练预测 Close的价格反归一化
y_train_pred = model(X_train_tensor).detach().numpy()
y_train_denormalized_predictions = (y_train_pred * scale_params[3]) + original_min[3]# 将测试预测 Close的价格反归一化
y_test_pred = model(X_test_tensor).detach().numpy()
y_test_denormalized_predictions = (y_test_pred * scale_params[3]) + original_min[3]

5.3 结果可视化

计算训练预测与测试预测的绘图数据

# shift train predictions for plotting
trainPredict = AAPL[window_size:X_train.shape[0]+X_train.shape[1]]
trainPredictPlot = trainPredict.assign(TrainPrediction= y_train_denormalized_predictions)testPredict = AAPL[X_train.shape[0]+X_train.shape[1]:]
testPredictPlot = testPredict.assign(TestPrediction= y_test_denormalized_predictions)

绘制模型收盘价格的原始数据与预测数据

# Visualize the data
plt.figure(figsize=(18,6))
plt.title('BiLSTM Close Price Validation')
plt.plot(AAPL['Close'], color='blue', label='original')
plt.plot(trainPredictPlot['TrainPrediction'], color='orange',label='Train Prediction')
plt.plot(testPredictPlot['TestPrediction'], color='red', label='Test Prediction')
plt.legend()
plt.show()

6. 模型预测

6.1 转换最新时间步收盘价的数组为张量

# 假设latest_closes是一个包含最新window_size个收盘价的列表或数组
latest_closes = AAPL[['Open', 'High', 'Low', 'Close']][-window_size:].values
scaled_latest_closes = scaler.transform(latest_closes)
tensor_latest_closes = torch.from_numpy(scaled_latest_closes).type(torch.Tensor).view(1, window_size, 4) 
print(tensor_latest_closes.shape)

torch.Size([1, 60, 4])

6.2 预测下一个时间点的收盘价格

# 使用模型预测下一个时间点的收盘价
next_close_pred = model(tensor_latest_closes).detach().numpy()
next_close_denormalized_pred = (next_close_pred * scale_params[3]) + original_min[3]
next_close_denormalized_pred

array([[165.9167]], dtype=float32)