时间序列信号分解涉及将原始时间序列分解为其组成部分。通过分解时间序列，我们可以深入了解其潜在模式、趋势-周期和季节性影响，从而提高理解和预测准确性。

本笔记将展示如何使用 NHITS/NBEATSx 提取这些时间序列的组成部分。我们将
- 安装 NeuralForecast。
- 模拟谐波信号。
- NHITS 的预测分解。
- NBEATSx 的预测分解。

您可以使用 Google Colab 在 GPU 上运行这些实验。

​
1. 安装 NeuralForecast

!pip install neuralforecast

​
2. 模拟谐波信号

在此示例中，我们将考虑一个包含两个频率的谐波信号：一个低频和一个高频。

import numpy as np
import pandas as pd

N = 10_000
T = 1.0 / 800.0 # sample spacing
x = np.linspace(0.0, N*T, N, endpoint=False)

y1 = np.sin(10.0 * 2.0*np.pi*x) 
y2 = 0.5 * np.sin(100 * 2.0*np.pi*x)
y = y1 + y2

import matplotlib.pyplot as plt
plt.rcParams["axes.grid"]=True

fig, ax = plt.subplots(figsize=(6, 2.5))
plt.plot(y[-80:], label='True')
plt.plot(y1[-80:], label='Low Frequency', alpha=0.4)
plt.plot(y2[-80:], label='High Frequency', alpha=0.4)
plt.ylabel('Harmonic Signal')
plt.xlabel('Time')
plt.legend()
plt.show()
plt.close()


# Split dataset into train/test
# Last horizon observations for test
horizon = 96
Y_df = pd.DataFrame(dict(unique_id=1, ds=np.arange(len(x)), y=y))
Y_train_df = Y_df.groupby('unique_id').head(len(Y_df)-horizon)
Y_test_df = Y_df.groupby('unique_id').tail(horizon)
Y_test_df
unique_iddsy
990419904-0.951057
990519905-0.570326
990619906-0.391007
990719907-0.499087
990819908-0.809017
…………
999519995-0.029130
999619996-0.309017
999719997-0.586999
999819998-0.656434
999919999-0.432012

​
3. NHITS 分解

我们将利用 NHITS 的堆栈特化来恢复潜在的谐波函数。

NHITS 是一种受小波启发的算法，可以将时间序列分解为各种尺度或分辨率，有助于识别局部模式或特征。每一层的表达比例可以控制模型的堆栈特化。

from neuralforecast.models import NHITS, NBEATSx
from neuralforecast import NeuralForecast
from neuralforecast.losses.pytorch import HuberLoss, MQLoss

models = [NHITS(h=horizon,                           # Forecast horizon
                input_size=2 * horizon,              # Length of input sequence
                loss=HuberLoss(),                    # Robust Huber Loss
                max_steps=1000,                      # Number of steps to train
                dropout_prob_theta=0.5,
                interpolation_mode='linear',
                stack_types=['identity']*2,
                n_blocks=[1, 1],
                mlp_units=[[64, 64],[64, 64]],
                n_freq_downsample=[10, 1],           # Inverse expressivity ratios for NHITS' stacks specialization
                val_check_steps=10,                  # Frequency of validation signal (affects early stopping)
              )
          ]
nf = NeuralForecast(models=models, freq=1)
nf.fit(df=Y_train_df)

from neuralforecast.tsdataset import TimeSeriesDataset

# NHITS decomposition plot
model = nf.models[0]
dataset, *_ = TimeSeriesDataset.from_df(df = Y_train_df)
y_hat = model.decompose(dataset=dataset)

GPU available: True (cuda), used: True
TPU available: False, using: 0 TPU cores
HPU available: False, using: 0 HPUs
LOCAL_RANK: 0 - CUDA_VISIBLE_DEVICES: [0]

Predicting: |          | 0/? [00:00<?, ?it/s]

fig, ax = plt.subplots(3, 1, figsize=(6, 7))

ax[0].plot(Y_test_df['y'].values, label='True', linewidth=4)
ax[0].plot(y_hat.sum(axis=1).flatten(), label='Forecast', color="#7B3841")
ax[0].legend()
ax[0].set_ylabel('Harmonic Signal')

ax[1].plot(y_hat[0,1]+y_hat[0,0], label='stack1', color="green")
ax[1].set_ylabel('NHITS Stack 1')

ax[2].plot(y_hat[0,2], label='stack2', color="orange")
ax[2].set_ylabel('NHITS Stack 2')
ax[2].set_xlabel(r'Prediction $\tau \in \{t+1,..., t+H\}$')
plt.show()

​
4. NBEATSx 分解

此处我们将利用 NBEATSx 可解释的基础投影来恢复潜在的谐波函数。

NBEATSx 的可解释变体顺序地将信号投影到多项式和谐波基上，以学习趋势 TT 和季节性组成部分：y^[t+1:t+H]=θ1T+θ2S\hat{y}_{[t+1:t+H]} = \theta_{1} T + \theta_{2} S

与 NHITS 的小波状投影不同，基础极大地决定了投影的行为。并且傅里叶投影无法立即分解为单个频率。

models = [NBEATSx(h=horizon,                           # Forecast horizon
                  input_size=2 * horizon,              # Length of input sequence
                  loss=HuberLoss(),                    # Robust Huber Loss
                  max_steps=1000,                      # Number of steps to train
                  dropout_prob_theta=0.5,
                  stack_types=['trend', 'seasonality'], # Harmonic/Trend projection basis
                  n_polynomials=0,                      # Lower frequencies can be captured by polynomials
                  n_blocks=[1, 1],
                  mlp_units=[[64, 64],[64, 64]],
                  val_check_steps=10,                  # Frequency of validation signal (affects early stopping)
              )
          ]
nf = NeuralForecast(models=models, freq=1)
nf.fit(df=Y_train_df)

# NBEATSx decomposition plot
model = nf.models[0]
dataset, *_ = TimeSeriesDataset.from_df(df = Y_train_df)
y_hat = model.decompose(dataset=dataset)

GPU available: True (cuda), used: True
TPU available: False, using: 0 TPU cores
HPU available: False, using: 0 HPUs
LOCAL_RANK: 0 - CUDA_VISIBLE_DEVICES: [0]

Predicting: |          | 0/? [00:00<?, ?it/s]

fig, ax = plt.subplots(3, 1, figsize=(6, 7))

ax[0].plot(Y_test_df['y'].values, label='True', linewidth=4)
ax[0].plot(y_hat.sum(axis=1).flatten(), label='Forecast', color="#7B3841")
ax[0].legend()
ax[0].set_ylabel('Harmonic Signal')

ax[1].plot(y_hat[0,1]+y_hat[0,0], label='stack1', color="green")
ax[1].set_ylabel('NBEATSx Trend Stack')

ax[2].plot(y_hat[0,2], label='stack2', color="orange")
ax[2].set_ylabel('NBEATSx Seasonality Stack')
ax[2].set_xlabel(r'Prediction $\tau \in \{t+1,..., t+H\}$')
plt.show()

​
参考文献