最强总结！十大时间序列预测模型 ！！

哈喽，我是小白~

今儿和大家再来聊聊时间序列预测模型！~

时间序列预测模型 能够捕捉数据中的趋势和季节性变化，从而实现对未来值的准确预测。通过分析历史数据，时间序列模型帮助决策者制定更明智的策略，优化资源配置。

最终，这些模型在金融、气象和供应链管理等多个领域提供了非常宝贵的洞察，提升了业务的响应能力和效率。

今儿涉及到的模型有：

• 自回归模型 (AR)
• 移动平均模型 (MA)
• 自回归移动平均模型 (ARMA)
• 自回归积分移动平均模型 (ARIMA)
• 季节性自回归积分移动平均模型 (SARIMA)
• 指数平滑法 (ETS, Exponential Smoothing)
• LSTM
• Facebook Prophet
• XGBoost
• LightGBM

给出每种算法介绍、原理、核心公式和公式详细的推导

1. 自回归模型

自回归模型是一种基于自身过去观测值的线性回归模型，主要假设时间序列的当前值是其前若干时刻值的线性组合以及噪声的函数。AR模型特别适用于描述具有平稳性特征的时间序列。

原理

AR模型认为时间序列的当前值  由前  个滞后值的线性组合表示，即：

其中：

•  是模型的待估参数，称为自回归系数；
•  是模型的阶数，即滞后步数；
•  是白噪声误差项，假设服从正态分布 。

核心公式

AR(p) 模型的公式可以表示为：

推导过程：

1. 时间序列的平稳性：

• AR模型适用于平稳时间序列。平稳序列意味着其均值和方差不随时间变化，且自相关函数仅依赖于时间间隔 。
• 时间序列  的平稳性要求其期望 、方差  和自相关函数  在不同时刻是相同的。

2. 模型表达形式：

考虑AR(1)模型，表达式为：

推广至一般的AR(p)模型：

这表明当前的观测值  是  的线性组合加上一个白噪声项 。

3. 参数估计（最小二乘法）：

• 设  为观测值向量，  为自变量矩阵，  为噪声项向量。
• 利用最小二乘法估计参数：

通过求解这个方程，可以得到AR模型的系数。

4. 稳定性条件：

对于AR模型的稳定性，要求多项式  的根都位于单位圆之外，这样模型才能确保平稳。

完整案例

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.ar_model import AutoReg
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.stattools import acf

# 生成虚拟时间序列数据（使用正弦波+随机噪声作为例子）
np.random.seed(42)
n = 3000
t = np.linspace(0, 50, n)
y = 10 * np.sin(0.2 * t) + np.random.normal(scale=2, size=n)

# 拆分训练和测试集
train_size = int(len(y) * 0.8)
train, test = y[:train_size], y[train_size:]

# 拟合自回归模型
model = AutoReg(train, lags=5)
model_fit = model.fit()

# 预测
predictions = model_fit.predict(start=len(train), end=len(train) + len(test) - 1, dynamic=False)

# 创建绘图
fig, ax = plt.subplots(2, 2, figsize=(14, 10))

# 图1：原始时间序列数据
ax[0, 0].plot(t, y, label='Time Series Data', color='red')
ax[0, 0].set_title('Original Time Series Data', fontsize=12)
ax[0, 0].set_xlabel('Time')
ax[0, 0].set_ylabel('Value')
ax[0, 0].legend()

# 图2：自相关函数 (ACF)
plot_acf(y, lags=30, ax=ax[0, 1], color='blue')
ax[0, 1].set_title('Autocorrelation Function (ACF)', fontsize=12)

# 图3：偏自相关函数 (PACF)
plot_pacf(y, lags=30, ax=ax[1, 0], color='green')
ax[1, 0].set_title('Partial Autocorrelation Function (PACF)', fontsize=12)

# 图4：模型残差及其自相关性
residuals = test - predictions
ax[1, 1].plot(residuals, label='Residuals', color='purple')
ax[1, 1].set_title('Residuals Time Series', fontsize=12)
ax[1, 1].set_xlabel('Time')
ax[1, 1].legend()

# 调整布局
plt.tight_layout()
plt.show()

1. 时间序列的原始数据图：红色的时间序列曲线，展示数据的总体趋势和噪声。这可以帮助我们直观地观察时间序列的模式和周期性。

2. 自相关函数图 (ACF)：自相关函数图展示了时间序列与其自身在不同滞后时刻的相关性。这有助于理解数据的平稳性及其滞后特征，帮助选择模型的滞后阶数（AR的p值）。

3. 偏自相关函数图 (PACF)：偏自相关函数可以展示每个滞后的影响，排除中间滞后的干扰，这更有利于选择适合的滞后阶数。

4. 残差的时间序列和自相关性：预测残差的自相关性和残差时间序列可以帮助我们判断模型的拟合质量。如果模型残差图像接近白噪声，说明模型较好地捕捉了数据的动态结构。

2. 移动平均模型

移动平均模型是一种基于随机误差项的线性模型，假定当前的观测值是过去若干个误差项的线性组合。它常用于描述非平稳时间序列的噪声成分。

原理

MA模型认为当前值  是前  个误差项的加权平均值。因此，MA模型是通过随机扰动项的线性组合来解释时间序列的波动。

核心公式

MA(q) 模型可以写为：

其中：

•  是待估参数，称为移动平均系数；
•  是白噪声误差项。

推导过程：

1. 模型形式：

MA(1)模型的简单形式为：

推广到一般的MA(q)模型为：

2. 误差项的独立性：

假设误差项  是一个零均值、方差为  的白噪声过程，即 ，它们之间是独立同分布的。

3. 参数估计（最大似然估计）：

• 假设误差项是服从正态分布的，可以通过最大似然估计（MLE）的方法估计参数 。
• 利用时间序列的自相关函数（ACF）进行模型的阶数选择，ACF 在  之后会截尾。

4. 稳定性条件：

为了确保MA模型的稳定性和可逆性，要求多项式  的根位于单位圆之外。

完整案例

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf

# 生成虚拟时间序列数据（带噪声的）
np.random.seed(42)
n = 1000
time = np.arange(n)
noise = np.random.normal(0, 1, n)
data = np.sin(0.2 * time) + noise  # 原始信号 + 噪声

# 将时间序列数据保存到DataFrame中
df = pd.DataFrame({'time': time, 'value': data})

# 构建并拟合 MA(2) 模型
model = ARIMA(df['value'], order=(0, 0, 2))
model_fit = model.fit()

# 预测值
df['predicted'] = model_fit.fittedvalues

# 计算残差
df['residuals'] = df['value'] - df['predicted']

# 设置图形颜色和样式
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 图1: 原始数据与预测数据对比
axes[0, 0].plot(df['time'], df['value'], color='blue', label='Original Data', linewidth=2)
axes[0, 0].plot(df['time'], df['predicted'], color='red', linestyle='--', label='MA(2) Predicted', linewidth=2)
axes[0, 0].set_title('Original vs Predicted')
axes[0, 0].legend(loc='upper right')

# 图2: 残差图
axes[0, 1].plot(df['time'], df['residuals'], color='green', label='Residuals', linewidth=2)
axes[0, 1].axhline(y=0, color='black', linestyle='--', linewidth=1)
axes[0, 1].set_title('Residuals Analysis')
axes[0, 1].legend(loc='upper right')

# 图3: ACF (自相关函数) 图
plot_acf(df['residuals'], lags=20, ax=axes[1, 0], color='magenta')
axes[1, 0].set_title('ACF of Residuals')

# 图4: PACF (偏自相关函数) 图
plot_pacf(df['residuals'], lags=20, ax=axes[1, 1], color='orange')
axes[1, 1].set_title('PACF of Residuals')

# 调整布局
plt.tight_layout()
plt.show()

1. 生成虚拟时间序列数据：使用 np.sin() 生成正弦波信号并加入高斯噪声。该信号代表我们的时间序列数据。

2. 构建 MA(2) 模型：使用 statsmodels 库中的 ARIMA 模型（这里 ARIMA 中的 p=0, d=0, q=2 对应 MA(2)）。

3. 预测与残差计算：对模型拟合后，生成预测值，并计算残差（原始值与预测值之间的差值）。

• 原始与预测数据对比图：直观对比模型对数据的拟合效果，鲜艳的蓝色和红色用来区分原始值与预测值。
• 残差分析图：分析残差是否存在系统性偏差，绿线为残差曲线。
• ACF 图：用于检查残差的自相关性，帮助判断残差是否独立。
• PACF 图：用于分析序列与滞后项之间的直接关系，有助于理解残差的特征。

3. 自回归移动平均模型

ARMA模型是将自回归模型 (AR) 和移动平均模型 (MA) 结合起来的模型，适用于平稳时间序列。它既包含了时间序列过去的值（自回归部分），也包含了过去误差项的影响（移动平均部分）。

原理

ARMA(p, q) 模型由两部分组成：

• AR(p) 部分通过滞后项  来解释当前值；
• MA(q) 部分通过误差项  来解释当前值。

核心公式

ARMA(p, q) 模型可以写为：

其中：

•  是自回归系数；
•  是移动平均系数；
•  是白噪声误差项。

推导过程：

1. 模型形式：

ARMA(1, 1)模型的形式为：

推广到一般的ARMA(p, q)模型为：

2. 混合自回归和移动平均：

• 自回归部分解释时间序列的惯性或趋势；
• 移动平均部分解释由于误差项引起的短期波动。

3. 参数估计（最大似然估计或最小二乘法）：

• 参数估计可以通过最大似然估计或基于时间序列的ACF和PACF（偏自相关函数）来选择合适的阶数  和 。

4. 稳定性和可逆性条件：

ARMA模型的稳定性要求AR部分满足平稳性条件（AR部分多项式的根位于单位圆之外），而可逆性条件要求MA部分多项式的根位于单位圆之外。

完整案例

import numpy as np
import matplotlib.pyplot as plt
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.graphics.tsaplots import plot_acf
from scipy import stats

# 生成虚拟数据
np.random.seed(42)
n = 2000
# 自回归和移动平均部分的参数
arparams = np.array([0.75, -0.25])
maparams = np.array([0.65])
# 白噪声部分
noise = np.random.normal(0, 1, n)
# 使用ARMA模型生成虚拟时间序列数据
data = np.zeros(n)
for i in range(2, n):
    data[i] = arparams[0] * data[i-1] + arparams[1] * data[i-2] + noise[i] + maparams[0] * noise[i-1]

# 拟合ARMA模型
model = ARIMA(data, order=(2, 0, 1))
arma_result = model.fit()

# 生成预测值
pred_start = 180
pred_end = 220
pred = arma_result.get_prediction(start=pred_start, end=pred_end)
pred_mean = pred.predicted_mean
pred_ci = pred.conf_int()

# 计算残差
residuals = arma_result.resid

# 创建子图
fig, ax = plt.subplots(2, 2, figsize=(14, 10))

# 1. 时间序列图：原始数据 vs 预测数据
ax[0, 0].plot(data, color='blue', label='原始数据')
ax[0, 0].plot(np.arange(pred_start, pred_end+1), pred_mean, color='red', label='预测数据')
ax[0, 0].fill_between(np.arange(pred_start, pred_end+1), 
                      pred_ci[:, 0], pred_ci[:, 1], color='pink', alpha=0.3)
ax[0, 0].set_title('时间序列图：原始数据 vs 预测数据')
ax[0, 0].legend()

# 2. 残差ACF图
plot_acf(residuals, ax=ax[0, 1], lags=20)
ax[0, 1].set_title('残差ACF图')

# 3. 残差QQ图
stats.probplot(residuals, dist="norm", plot=ax[1, 0])
ax[1, 0].set_title('残差QQ图')

# 4. 残差的时间序列图
ax[1, 1].plot(residuals, color='green')
ax[1, 1].axhline(y=0, linestyle='--', color='gray')
ax[1, 1].set_title('残差的时间序列图')

# 调整布局
plt.tight_layout()
plt.show()

1. 数据生成：我们首先生成一个虚拟的时间序列数据，该数据符合ARMA(2, 1)模型。自回归部分参数 arparams 为 [0.75, -0.25]，移动平均部分参数 maparams 为 [0.65]，并加入随机噪声。

2. 模型拟合：使用 ARIMA 模型（其中AR和MA阶分别为2和1）对生成的虚拟数据进行拟合。

3. 预测：从序列的第180点开始，进行步数为40的预测，得到预测的均值以及置信区间。

4. 残差计算：使用模型的拟合结果，计算残差（实际值减去预测值），以评估模型的拟合效果。

• 时间序列图（左上图）：展示原始时间序列数据和模型预测的对比，预测结果带有置信区间。
• 残差ACF图（右上图）：自相关函数显示残差是否仍存在显著的自相关性，理想情况下应该接近0。
• 残差QQ图（左下图）：用于评估残差是否符合正态分布，残差点应接近于直线。
• 残差的时间序列图（右下图）：展示残差随时间的变化，理想情况下，残差应呈现白噪声（无明显模式）。

4. 自回归积分移动平均模型

ARIMA模型是ARMA模型的扩展，适用于非平稳时间序列。通过对时间序列进行差分（积分部分）操作，使其变得平稳，再对其应用ARMA模型。

原理

ARIMA模型由三个部分组成：

• 自回归 (AR) 部分；
• 差分 (I, Integration) 部分，用来处理非平稳数据；
• 移动平均 (MA) 部分。

核心公式

ARIMA(p, d, q) 模型可以表示为：

其中：

•  是滞后算子 ；
•  是差分阶数；
•  和  是AR和MA部分的系数。

推导过程

1. 差分后的平稳时间序列：

假设原时间序列  是非平稳的，我们对其进行  次差分得到新的序列 ，使得  是平稳的。之后可以对这个平稳的序列  应用 ARMA(p, q) 模型。对于  次差分后的序列，ARIMA 模型的形式为：

其中 。

2. 滞后算子表示：

差分操作可以用滞后算子  来表示。滞后算子  表示将时间序列滞后  个时间单位。差分操作可以表示为：

对于  次差分，表示为：

因此，ARIMA(p, d, q) 模型的公式可以写为：

3. 参数估计：

ARIMA 模型的参数包括 、、 三个部分：

• ：自回归项的阶数；
• ：差分的次数，用于处理非平稳性；
• ：移动平均项的阶数。

这些参数通常可以通过分析时间序列的自相关函数 (ACF) 和偏自相关函数 (PACF) 来估计。通过最大似然估计或非线性最小二乘法，可以估计模型的参数  和 。

4. 稳定性和可逆性条件：

与 ARMA 模型类似，ARIMA 模型中的 AR 部分要求满足平稳性条件，MA 部分要求满足可逆性条件。差分操作（即积分部分）可以将非平稳序列转换为平稳序列。

完整案例

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.arima.model import ARIMA
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.graphics.gofplots import qqplot
import seaborn as sns

# 设置随机种子
np.random.seed(42)

# 生成时间序列数据
n = 2000
time = np.arange(n)
trend = time * 0.05  # 增加一个线性趋势
seasonal = 10 * np.sin(2 * np.pi * time / 50)  # 增加季节性成分
noise = np.random.normal(0, 2, n)  # 添加噪声
data = trend + seasonal + noise

# 转为 DataFrame
df = pd.DataFrame({'Time': time, 'Value': data})

# 拟合 ARIMA 模型 (ARIMA(2,1,2))
model = ARIMA(df['Value'], order=(2, 1, 2))
fit = model.fit()

# 预测未来20步
forecast = fit.forecast(steps=20)
forecast_index = np.arange(n, n + 20)

# 绘制图形
fig, ax = plt.subplots(3, 2, figsize=(14, 12))

# 1. 原始数据和拟合值
ax[0, 0].plot(df['Time'], df['Value'], label='Original Data', color='blue', linewidth=2)
ax[0, 0].plot(df['Time'], fit.fittedvalues, label='Fitted Values', color='red', linestyle='--', linewidth=2)
ax[0, 0].legend(loc='upper left')
ax[0, 0].set_title('Original Data vs Fitted Values')

# 2. ACF 图
plot_acf(df['Value'], ax=ax[0, 1], lags=40)
ax[0, 1].set_title('ACF of Original Data')

# 3. PACF 图
plot_pacf(df['Value'], ax=ax[1, 0], lags=40)
ax[1, 0].set_title('PACF of Original Data')

# 4. 残差 QQ 图
qqplot(fit.resid, line='s', ax=ax[1, 1])
ax[1, 1].set_title('QQ Plot of Residuals')

# 5. 残差直方图
sns.histplot(fit.resid, bins=20, kde=True, color='green', ax=ax[2, 0])
ax[2, 0].set_title('Residuals Histogram')

# 6. 预测结果与原始数据对比
ax[2, 1].plot(df['Time'], df['Value'], label='Original Data', color='blue', linewidth=2)
ax[2, 1].plot(forecast_index, forecast, label='Forecast', color='orange', linestyle='--', linewidth=2)
ax[2, 1].legend(loc='upper left')
ax[2, 1].set_title('Forecast vs Original Data')

# 调整布局
plt.tight_layout()
plt.show()

1. 原始数据 vs 模型拟合值：展示模型对原始时间序列的拟合效果，便于观察模型的表现。通过观察红色虚线拟合值与蓝色实际数据的贴合程度，了解模型的准确度。

2. ACF (自相关函数) 图：用来分析数据的自相关结构。自相关函数帮助我们理解时间序列中数据的依赖关系，特别是用于确定 ARIMA 模型中 q 参数（移动平均项的阶数）。

3. PACF (部分自相关函数) 图：部分自相关函数帮助我们识别数据中直接的滞后关系，通常用于确定 ARIMA 模型中 p 参数（自回归项的阶数）。

4. 残差 QQ 图：通过比较残差的分布与标准正态分布，检查模型拟合的残差是否服从正态分布。这是判断模型残差是否是白噪声的一个重要指标。

5. 残差直方图：进一步分析残差的分布情况，检查残差是否呈现出中心对称的分布，以判断模型的误差是否满足正态分布假设。

6. 预测结果与原始数据对比：将预测值与原始时间序列进行比较，展示模型的预测效果。可以直观地看出预测的未来数据是否与现有趋势吻合。

5. 季节性自回归积分移动平均模型

SARIMA 是 ARIMA 模型的扩展，特别适合具有季节性波动的时间序列。它结合了季节性成分和非季节性成分，用以捕捉时间序列的周期性行为。

原理

SARIMA 模型通过增加季节性差分、自回归和移动平均项来捕捉季节性趋势。它的形式为 ARIMA 模型与一个季节性 ARIMA 模型的组合。

核心公式

SARIMA() × () 模型表示为：

其中：

• ：非季节性部分的自回归、差分和移动平均阶数；
• ：季节性部分的自回归、差分和移动平均阶数；
• ：季节周期的长度（例如，季节性周期为 12 对于月度数据）。

推导过程：

1. 季节性部分的差分：

季节性差分  用于消除季节性影响。对  表示周期的时间序列，可以进行季节性差分操作。对于季节性 AR 和 MA 成分，滞后算子  将周期滞后  个单位。

2. 模型分解：

SARIMA 模型可以看作是非季节性 ARIMA(p, d, q) 模型和季节性 ARIMA() 模型的组合。非季节性部分处理时间序列中的整体趋势，而季节性部分捕捉周期性波动。

3. 参数估计：

与 ARIMA 类似，SARIMA 的参数也可以通过分析 ACF 和 PACF 来选择。由于季节性组件的引入，模型更加复杂，通常需要使用数值优化方法（例如最大似然估计）来估计所有参数。

完整案例

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import seaborn as sns
from statsmodels.tsa.statespace.sarimax import SARIMAX
from statsmodels.graphics.tsaplots import plot_acf, plot_pacf
from statsmodels.tsa.seasonal import seasonal_decompose

# 设置随机种子，便于结果重现
np.random.seed(42)

# 生成虚拟时间序列数据（带有趋势和季节性）
n_periods = 2000
time = np.arange(n_periods)
seasonal_effect = 10 * np.sin(2 * np.pi * time / 24)  # 周期为24的季节性成分
trend = 0.05 * time  # 线性趋势
noise = np.random.normal(0, 2, n_periods)  # 随机噪声

# 生成时间序列数据
data = trend + seasonal_effect + noise
dates = pd.date_range(start='2020-01-01', periods=n_periods, freq='D')
time_series = pd.Series(data, index=dates)

# 构建 SARIMA 模型
model = SARIMAX(time_series, order=(1, 1, 1), seasonal_order=(1, 1, 1, 24))
results = model.fit()

# 预测未来数据
pred_start = n_periods - 30
pred_end = n_periods + 30
forecast = results.get_prediction(start=pred_start, end=pred_end)
forecast_ci = forecast.conf_int()

# 分解时间序列，得到趋势、季节性和残差
decomposition = seasonal_decompose(time_series, model='additive', period=24)

# 绘制图形
plt.figure(figsize=(16, 12))
sns.set(style="whitegrid")

# 子图 1: 原始时间序列及趋势、季节性分解
plt.subplot(3, 1, 1)
plt.plot(time_series, label='Observed', color='blue', linewidth=2)
plt.plot(decomposition.trend, label='Trend', color='red', linestyle='--', linewidth=2)
plt.plot(decomposition.seasonal, label='Seasonal', color='green', linestyle=':', linewidth=2)
plt.title('Time Series Decomposition (Observed, Trend, Seasonal)', fontsize=14)
plt.legend()

# 子图 2: 残差的自相关图 (ACF) 和部分自相关图 (PACF)
plt.subplot(3, 2, 3)
plot_acf(results.resid, lags=40, ax=plt.gca(), title="Residual ACF", color='magenta')

plt.subplot(3, 2, 4)
plot_pacf(results.resid, lags=40, ax=plt.gca(), title="Residual PACF", color='orange')

# 子图 3: 实际值与预测值对比
plt.subplot(3, 1, 3)
plt.plot(time_series, label='Actual', color='blue', linewidth=2)
plt.plot(forecast.predicted_mean, label='Forecast', color='purple', linestyle='--', linewidth=2)
plt.fill_between(forecast_ci.index, forecast_ci.iloc[:, 0], forecast_ci.iloc[:, 1], color='pink', alpha=0.3)
plt.axvline(time_series.index[pred_start], color='gray', linestyle='--', linewidth=1, label='Forecast Start')
plt.title('Forecast vs Actual', fontsize=14)
plt.legend()

plt.tight_layout()
plt.show()

1. 原始时间序列及分解图：这部分图形通过分解时间序列成趋势和季节性分量来分析时间序列的组成部分。原始数据呈现的趋势和季节性帮助我们更好地理解时间序列结构。

2. 残差的 ACF 和 PACF 图：这些图用于检查残差的自相关性。理想情况下，残差应该是白噪声（没有显著的自相关），如果发现残差有显著的自相关性，说明模型可能存在改进空间。

3. 实际值与预测值对比图：展示模型对时间序列的预测能力，预测区间以粉色阴影表示，这个图让我们清晰地看到模型的预测是否贴合实际数据。

6. 指数平滑法

ETS（指数平滑）是一类广泛应用的平滑方法，适合平稳、趋势和季节性时间序列。它的基本思想是根据时间序列的加权平均来预测未来的值，其中最近的观测值权重更大。

原理

ETS 模型通过递归平滑来更新预测值。不同的 ETS 模型包括：

• 错误项：加法或乘法；
• 趋势：无趋势、加法趋势、乘法趋势；
• 季节性：无季节性、加法季节性、乘法季节性。

核心公式

单一指数平滑法（用于无趋势无季节性数据）公式为：

其中  是平滑参数，。

推导过程：

1. 递归更新： 预测值  是当前观测值  和前一时刻预测值  的加权平均值。权重  控制了对当前观测值的重视程度。

2. 扩展到双指数平滑： 对于具有趋势的时间序列，可以引入趋势项 ，双指数平滑公式为：

其中  表示当前水平， 表示趋势。它们分别通过平滑递归更新。

3. 扩展到 Holt-Winters 季节性模型： 对于季节性数据，Holt-Winters 法则通过三个递归更新公式来预测未来值：

其中  表示季节性因素， 和  分别为水平和趋势项。

完整案例

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.holtwinters import ExponentialSmoothing

# 生成虚拟时间序列数据
np.random.seed(42)
n_periods = 1000
trend = np.linspace(10, 100, n_periods)
seasonal = 10 * np.sin(np.linspace(0, 3 * np.pi, n_periods))
noise = np.random.normal(0, 2, n_periods)
data = trend + seasonal + noise

# 构建数据集
date_range = pd.date_range(start='2020-01-01', periods=n_periods, freq='M')
ts_data = pd.Series(data, index=date_range)

# 拟合 ETS 模型 (带趋势和季节性的指数平滑)
model = ExponentialSmoothing(ts_data, trend='add', seasonal='add', seasonal_periods=12).fit()

# 进行预测（未来 24 个周期）
forecast = model.forecast(24)
fitted_values = model.fittedvalues

# 绘图
fig, axes = plt.subplots(2, 2, figsize=(14, 10), dpi=100)
fig.suptitle('ETS (Exponential Smoothing) Time Series Analysis', fontsize=16)

# 图 1: 原始数据和模型拟合
axes[0, 0].plot(ts_data.index, ts_data, label='Original Data', color='blue', linestyle='--')
axes[0, 0].plot(ts_data.index, fitted_values, label='Fitted Values', color='red')
axes[0, 0].set_title('Original Data vs Fitted Values')
axes[0, 0].legend()
axes[0, 0].set_xlabel('Date')
axes[0, 0].set_ylabel('Values')

# 图 2: 预测与置信区间
forecast_index = pd.date_range(start=ts_data.index[-1], periods=25, freq='M')[1:]
axes[0, 1].plot(ts_data.index, ts_data, label='Original Data', color='blue', linestyle='--')
axes[0, 1].plot(forecast_index, forecast, label='Forecast', color='green')
axes[0, 1].fill_between(forecast_index, forecast - 1.96 * model.sigma, forecast + 1.96 * model.sigma, color='lightgreen', alpha=0.3, label='95% Confidence Interval')
axes[0, 1].set_title('Forecast with Confidence Interval')
axes[0, 1].legend()
axes[0, 1].set_xlabel('Date')
axes[0, 1].set_ylabel('Values')

# 图 3: 残差分析
residuals = ts_data - fitted_values
axes[1, 0].plot(ts_data.index, residuals, label='Residuals', color='purple')
axes[1, 0].set_title('Residuals (Fitted Errors)')
axes[1, 0].axhline(0, linestyle='--', color='black', linewidth=1)
axes[1, 0].set_xlabel('Date')
axes[1, 0].set_ylabel('Residuals')
axes[1, 0].legend()

# 图 4: 残差的自相关性
pd.plotting.autocorrelation_plot(residuals, ax=axes[1, 1], color='orange')
axes[1, 1].set_title('Residuals Autocorrelation')

plt.tight_layout(rect=[0, 0, 1, 0.96])
plt.show()

1. 原始数据 vs 拟合数据：通过此图，可以看出模型在历史数据上的拟合情况，是否能够很好地捕捉数据的趋势和季节性。

2. 预测和置信区间：此图显示模型对未来的预测，并通过阴影区表示预测的95%置信区间。阴影部分越宽，表明模型预测的不确定性越大。

3. 残差图：残差应该随机分布，没有明显的趋势或周期性。如果有趋势或模式，说明模型未能捕捉到所有的信息。

4. 残差的自相关性：该图可以揭示残差的自相关结构。如果残差存在自相关性，说明模型仍然可以进一步改进。

7. 长短期记忆网络

LSTM 是一种特殊的递归神经网络（RNN），设计用于捕捉时间序列中的长期和短期依赖关系。它克服了标准 RNN 中的梯度消失问题，使其能够有效地处理长序列数据。

原理

LSTM 通过引入“记忆单元”和门结构（输入门、遗忘门、输出门）来控制信息的存储、更新和输出。这种门结构允许模型选择性地记住和遗忘信息。

核心公式

LSTM 单元的核心计算步骤为： 1. 遗忘门：

决定哪些信息需要丢弃。

2. 输入门：

决定哪些新信息需要写入记忆单元。

3. 候选记忆：

4. 更新记忆：

5. 输出门：

6. 更新隐藏状态：

推导过程：

LSTM 的设计使得它可以通过遗忘门控制记忆信息的保留，通过输入门控制当前输入信息的保存，通过输出门控制输出信息的释放，从而有效地捕捉长时间依赖关系。

完整案例

import numpy as np
import matplotlib.pyplot as plt
from sklearn.preprocessing import MinMaxScaler
from tensorflow.keras.models import Sequential
from tensorflow.keras.layers import LSTM, Dense
from tensorflow.keras.optimizers import Adam

# 生成虚拟时间序列数据集
np.random.seed(0)
time_steps = 1000
data = np.sin(np.linspace(0, 20, time_steps)) + np.random.normal(scale=0.1, size=time_steps)

# 数据归一化
scaler = MinMaxScaler(feature_range=(0, 1))
data_scaled = scaler.fit_transform(data.reshape(-1, 1))

# 准备训练集和测试集
X, y = [], []
look_back = 10
for i in range(len(data_scaled) - look_back):
    X.append(data_scaled[i:i + look_back, 0])
    y.append(data_scaled[i + look_back, 0])

X, y = np.array(X), np.array(y)
X = np.reshape(X, (X.shape[0], X.shape[1], 1))

# 构建LSTM模型
model = Sequential()
model.add(LSTM(units=50, return_sequences=False, input_shape=(look_back, 1)))
model.add(Dense(units=1))

# 编译模型
model.compile(optimizer=Adam(learning_rate=0.01), loss='mean_squared_error')

# 训练模型
history = model.fit(X, y, epochs=50, batch_size=16, validation_split=0.1, verbose=0)

# 生成预测数据
predicted = model.predict(X)

# 将预测值反归一化
predicted_inverse = scaler.inverse_transform(predicted)
y_inverse = scaler.inverse_transform(y.reshape(-1, 1))

# 图像绘制
plt.figure(figsize=(14, 8))

# 子图1: 训练数据 vs 预测数据
plt.subplot(2, 1, 1)
plt.plot(range(look_back, len(data_scaled)), y_inverse, color='blue', label='Actual', linewidth=2)
plt.plot(range(look_back, len(data_scaled)), predicted_inverse, color='red', linestyle='--', label='Predicted', linewidth=2)
plt.title('LSTM Prediction vs Actual', fontsize=16)
plt.xlabel('Time Step', fontsize=12)
plt.ylabel('Value', fontsize=12)
plt.legend()

# 子图2: 损失曲线
plt.subplot(2, 1, 2)
plt.plot(history.history['loss'], color='green', label='Training Loss', linewidth=2)
plt.plot(history.history['val_loss'], color='orange', linestyle='--', label='Validation Loss', linewidth=2)
plt.title('Training and Validation Loss Curve', fontsize=16)
plt.xlabel('Epochs', fontsize=12)
plt.ylabel('Loss', fontsize=12)
plt.legend()

# 显示图像
plt.tight_layout()
plt.show()

1. 虚拟数据生成：生成了一个基于正弦波和高斯噪声的虚拟时间序列数据集。

2. 数据预处理：使用 MinMaxScaler 将数据归一化到 [0,1]，并且创建了具有 look_back 窗口的输入特征集。

3. LSTM模型：包含一个 LSTM 层和一个全连接层，使用均方误差（mean_squared_error）作为损失函数进行训练。

4. 训练过程：使用 Adam 优化器进行50轮的训练，并记录了训练和验证损失。

• 第一张图：展示真实值和模型预测值的对比。
• 第二张图：展示训练损失和验证损失随训练轮次的变化情况，帮助我们观察模型是否有过拟合的趋势。

8. Facebook Prophet

Facebook Prophet 是一种用于时间序列预测的模型，设计目标是能够处理带有缺失数据、季节性趋势的非平稳序列。它具有可解释性强、易于使用的特点。

原理

Prophet 模型基于加性模型，将趋势、季节性和节假日效应分解为不同的分量。

核心公式

Prophet 模型的形式为：

其中：

• ：趋势函数，用于捕捉长期的趋势变化；
• ：季节性函数，捕捉周期性波动；
• ：节假日效应，捕捉特殊事件的影响；
• ：误差项。

推导过程：

1. 趋势分解：

Prophet 提供了两种趋势模型：线性趋势和对数趋势。趋势函数  是分段线性回归模型。

2. 季节性分解：

使用傅里叶级数近似季节性函数 ，能够捕捉到复杂的周期性行为。

3. 节假日效应：

特殊事件可以作为节假日效应  加入模型，捕捉特定时间节点的非周期性变化。

完整案例

import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
from fbprophet import Prophet

# 生成虚拟时间序列数据
np.random.seed(42)
date_rng = pd.date_range(start='2020-01-01', end='2024-01-01', freq='D')
data = np.sin(np.linspace(0, 50, len(date_rng))) + np.random.normal(0, 0.2, len(date_rng)) + np.linspace(0, 5, len(date_rng))

df = pd.DataFrame({'ds': date_rng, 'y': data})

# 创建并拟合 Prophet 模型
model = Prophet(yearly_seasonality=True, weekly_seasonality=True, daily_seasonality=False)
model.fit(df)

# 预测未来365天的数据
future = model.make_future_dataframe(periods=365)
forecast = model.predict(future)

# 绘制结果
fig, axes = plt.subplots(2, 2, figsize=(14, 10))

# 图1：实际数据与预测数据对比
axes[0, 0].plot(df['ds'], df['y'], label='Actual Data', color='blue', linewidth=2)
axes[0, 0].plot(forecast['ds'], forecast['yhat'], label='Forecast', color='orange', linestyle='--', linewidth=2)
axes[0, 0].fill_between(forecast['ds'], forecast['yhat_lower'], forecast['yhat_upper'], color='orange', alpha=0.2)
axes[0, 0].set_title('Actual Data vs Forecast')
axes[0, 0].legend()

# 图2：预测的不确定性区间（上下限）
axes[0, 1].plot(forecast['ds'], forecast['yhat'], label='Forecast', color='green', linestyle='-', linewidth=2)
axes[0, 1].fill_between(forecast['ds'], forecast['yhat_lower'], forecast['yhat_upper'], color='lightgreen', alpha=0.3, label='Uncertainty Interval')
axes[0, 1].set_title('Forecast with Uncertainty Interval')
axes[0, 1].legend()

# 图3：季节性分解图 - 年度季节性
model.plot_components(forecast, ax=axes[1, 0])
axes[1, 0].set_title('Seasonality Decomposition')

# 图4：残差分析
residuals = df['y'] - forecast['yhat'][:len(df)]
axes[1, 1].hist(residuals, bins=30, color='purple', alpha=0.7)
axes[1, 1].set_title('Residual Analysis')
axes[1, 1].set_xlabel('Residuals')
axes[1, 1].set_ylabel('Frequency')

plt.tight_layout()
plt.show()

1. 实际数据与预测数据对比图：

• 目的：帮助我们直观地看到模型的预测效果，与实际数据的对比情况。通过该图可以查看模型拟合的好坏，并可以初步发现模型是否存在较大的偏差。
• 图中内容：蓝色线代表实际观测数据，橙色虚线代表模型预测值，阴影区域表示预测的上下不确定区间。

2. 预测的不确定性区间（上下限）图：

• 目的：展示预测值的置信区间，以便我们了解模型的预测不确定性。这对于评估模型的稳定性非常有帮助，尤其是在预测未来时。
• 图中内容：绿色线代表预测值，绿色阴影区域表示上下限的不确定性区间。

3. 季节性分解图：

• 目的：展示时间序列数据的年度、每周等季节性成分，帮助我们更好地理解数据的周期性模式。不同的时间序列可能会有不同的季节性（如年度周期、季度周期等）。
• 图中内容：显示趋势、年度和每周的季节性成分，揭示数据中的周期性波动。

4. 残差分析图：

• 目的：查看预测误差的分布，帮助评估模型是否存在系统性误差或偏差。残差应呈现随机分布，若存在明显的模式或偏移，说明模型可能存在未捕捉到的信息。
• 图中内容：柱状图展示了残差的频率分布。

9. XGBoost

XGBoost 是一种梯度提升决策树模型，常用于回归和分类任务，因其高效性和可扩展性而流行。它通过逐步拟合误差来提高模型性能。

原理

XGBoost 基于加法模型，构建一系列的弱学习器（决策树），通过梯度下降算法不断优化目标函数，逐步降低误差。

核心公式

XGBoost 的目标函数为：

其中：

•  是损失函数，衡量预测值与真实值之间的误差；
•  是正则项，用于控制模型的复杂度。

推导过程：

XGBoost 的推导基于梯度提升树的思想，通过贪心算法不断优化决策树的分裂点，并利用二阶导数信息加速优化过程。

完整案例

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from xgboost import XGBRegressor
from sklearn.metrics import mean_squared_error
from sklearn.model_selection import train_test_split

# 设置随机种子保证可复现性
np.random.seed(42)

# 生成虚拟时间序列数据集（每日时间点）
dates = pd.date_range('2023-01-01', periods=200)
data = pd.DataFrame({
    'Date': dates,
    'Value': np.sin(np.linspace(0, 20, 200)) + np.random.normal(0, 0.5, 200)
})

# 构建特征和目标变量
data['Day'] = data['Date'].dt.day
data['Month'] = data['Date'].dt.month
data['Year'] = data['Date'].dt.year

# 移动窗口创建滞后特征 (使用前一天的值作为特征)
data['Lag1'] = data['Value'].shift(1)
data['Lag2'] = data['Value'].shift(2)
data['Lag3'] = data['Value'].shift(3)

# 移除缺失值
data = data.dropna()

# 将日期列移除，因为它不是数值型特征
X = data[['Lag1', 'Lag2', 'Lag3', 'Day', 'Month', 'Year']]
y = data['Value']

# 拆分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, shuffle=False)

# 训练XGBoost模型
model = XGBRegressor(n_estimators=100, learning_rate=0.1, max_depth=3, random_state=42)
model.fit(X_train, y_train)

# 预测
y_pred = model.predict(X_test)

# 计算误差
mse = mean_squared_error(y_test, y_pred)
print(f'Mean Squared Error: {mse}')

# 创建图形
fig, axes = plt.subplots(3, 1, figsize=(12, 15))

# 图1: 原始数据时间序列趋势图
axes[0].plot(data['Date'], data['Value'], label='True Value', color='b', lw=2)
axes[0].set_title('Time Series Trend')
axes[0].set_xlabel('Date')
axes[0].set_ylabel('Value')
axes[0].legend()

# 图2: 预测 vs 实际值
axes[1].plot(data['Date'].iloc[len(y_train):], y_test, label='True Value', color='g', lw=2)
axes[1].plot(data['Date'].iloc[len(y_train):], y_pred, label='Predicted Value', color='r', lw=2, linestyle='--')
axes[1].set_title('Actual vs Predicted')
axes[1].set_xlabel('Date')
axes[1].set_ylabel('Value')
axes[1].legend()

# 图3: 残差图
residuals = y_test - y_pred
axes[2].plot(data['Date'].iloc[len(y_train):], residuals, label='Residuals', color='m', lw=2)
axes[2].axhline(0, color='black', linestyle='--', lw=2)
axes[2].set_title('Residual Analysis')
axes[2].set_xlabel('Date')
axes[2].set_ylabel('Residuals')
axes[2].legend()

plt.tight_layout()
plt.show()

1. 时序趋势图（图1）：展示了原始数据的时间序列趋势，帮助我们直观地观察数据的周期性或趋势变化（如上升或下降趋势）。这对于时间序列分析至关重要，因为时间序列的趋势性可以影响预测模型的性能。

2. 实际值 vs 预测值（图2）：展示了模型预测值与实际值的对比，通过这个图可以直观看出模型在时间序列上的拟合程度。如果预测值与实际值之间差距较大，可能意味着模型有待进一步优化。

3. 残差分析图（图3）：残差图展示了模型预测值与实际值之间的误差（即残差），如果残差图没有明显的模式或趋势，说明模型很好地捕捉了数据的特性。如果残差呈现系统性趋势或非随机性，可能意味着模型忽略了某些重要特性。

10. LightGBM

LightGBM 是另一种梯度提升决策树模型，设计目标是提高效率、降低内存使用。它特别适合处理大规模数据。

原理

LightGBM 通过叶节点分裂策略和基于直方图的决策树算法加速模型的训练，同时保持模型的高性能。

核心公式

LightGBM 的目标函数与 XGBoost 类似，也包括损失函数和正则项：

推导过程：

与 XGBoost 相比，LightGBM 的不同在于：

1. 基于直方图的决策树学习：通过将特征值分箱，将数据离散化，减少了分裂点的计算时间。

2. 叶节点增长策略：选择增益最大的叶节点进行分裂，避免传统决策树逐层分裂的低效问题。

完整案例

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
import lightgbm as lgb
from sklearn.metrics import mean_squared_error, mean_absolute_error
from sklearn.model_selection import train_test_split

# 创建虚拟时间序列数据
np.random.seed(42)
n = 2000
date_rng = pd.date_range(start='2020-01-01', end='2022-07-19', freq='D')
data = np.sin(np.linspace(0, 20, n)) + np.random.normal(0, 0.5, n)  # 添加噪声的正弦波
df = pd.DataFrame(data, index=date_rng, columns=['value'])
df['date'] = df.index
df['dayofweek'] = df['date'].dt.dayofweek
df['month'] = df['date'].dt.month
df['year'] = df['date'].dt.year

# 特征与目标
X = df[['dayofweek', 'month', 'year']]
y = df['value']

# 划分训练集和测试集
X_train, X_test, y_train, y_test = train_test_split(X, y, test_size=0.2, random_state=42, shuffle=False)

# 创建 LightGBM 数据集
train_data = lgb.Dataset(X_train, label=y_train)
test_data = lgb.Dataset(X_test, label=y_test)

# 设置 LightGBM 参数
params = {
    'objective': 'regression',
    'metric': 'rmse',
    'boosting_type': 'gbdt',
    'learning_rate': 0.1,
    'num_leaves': 31,
    'verbose': -1
}

# 训练模型
model = lgb.train(params, train_data, num_boost_round=100)

# 进行预测
y_pred = model.predict(X_test)

# 计算评估指标
rmse = np.sqrt(mean_squared_error(y_test, y_pred))
mae = mean_absolute_error(y_test, y_pred)
print(f'RMSE: {rmse:.4f}, MAE: {mae:.4f}')

# 计算残差
residuals = y_test.values - y_pred

# 绘制分析图形
plt.figure(figsize=(14, 8))

# 实际值与预测值图
plt.subplot(2, 2, 1)
plt.plot(df.index[-len(y_test):], y_test, label='Actual', color='blue')
plt.plot(df.index[-len(y_test):], y_pred, label='Predicted', color='orange')
plt.title('Actual vs Predicted')
plt.xlabel('Date')
plt.ylabel('Value')
plt.legend()

# 残差图
plt.subplot(2, 2, 2)
plt.scatter(y_pred, residuals, color='green')
plt.axhline(y=0, color='red', linestyle='--')
plt.title('Residuals Plot')
plt.xlabel('Predicted Values')
plt.ylabel('Residuals')

# 预测误差直方图
plt.subplot(2, 2, 3)
plt.hist(residuals, bins=20, color='purple', alpha=0.7)
plt.title('Histogram of Residuals')
plt.xlabel('Residuals')
plt.ylabel('Frequency')

# 散点图：实际值与预测值的关系
plt.subplot(2, 2, 4)
plt.scatter(y_test, y_pred, color='cyan')
plt.plot([-2, 2], [-2, 2], color='red', linestyle='--')  # 45度参考线
plt.title('Actual vs Predicted Scatter Plot')
plt.xlabel('Actual Values')
plt.ylabel('Predicted Values')

plt.tight_layout()
plt.show()

1. 数据生成: 使用正弦波加上随机噪声生成时间序列数据，模拟真实世界中的波动。

2. 特征工程: 从日期中提取出星期几、月份和年份作为特征。这些特征可能影响时间序列的趋势和季节性。

3. LightGBM模型训练: 选择 LightGBM 作为回归模型，设置超参数，并训练模型。

4. 模型评估:使用 RMSE（均方根误差）和 MAE（平均绝对误差）来评估模型性能。

• 实际值与预测值: 通过这张图可以直观地看到模型的预测效果。
• 残差图: 通过观察残差分布，可以判断模型是否存在系统性偏差。
• 预测误差直方图: 了解模型的误差分布情况，帮助识别潜在的问题。
• 实际值与预测值的散点图: 通过绘制实际值与预测值的关系，帮助我们判断模型的准确性。