用Python实现时间序列模型实战——Day 12: 状态空间模型

news2024/12/25 13:15:44

一、学习内容

1. 状态空间模型的基本概念

状态空间模型是一种用于时间序列分析的强大工具，能够描述具有潜在状态动态变化的系统。该模型通过显式地建模时间序列中的潜在状态（即隐藏变量），能够捕捉复杂的动态结构，适用于诸如动态系统、金融市场预测、信号处理等场景。

状态空间模型包含两个主要方程：

状态方程 (State Equation)：描述隐藏状态如何随时间演化。

$x_t = F x_{t-1} + \epsilon_t$

其中：

$x_t$ 是时间 $t$ 的隐藏状态向量。
$F$ 是状态转移矩阵，描述状态随时间的演化方式。
$\epsilon_t$ 是高斯噪声（白噪声）向量。

观测方程 (Observation Equation)：描述如何根据隐藏状态生成观测数据。

$y_t = H x_t + \nu_t$

其中：

$y_t$ 是时间 $t$ 的观测值。
$H$ 是观测矩阵，描述隐藏状态如何生成观测数据。
$\nu_t$ 是观测噪声（白噪声）向量。

2. 状态空间模型的构建与卡尔曼滤波

卡尔曼滤波 (Kalman Filter) 是一种递归算法，用于估计状态空间模型中的隐藏状态。它通过更新当前的状态估计，结合新的观测数据，生成最优估计。

卡尔曼滤波主要分为两个步骤：

预测步骤：根据当前的状态估计和状态方程预测下一时刻的状态。

$\hat{x}_t = F \hat{x}_{t-1}$

$P_t = F P_{t-1} F^T + Q$

其中：

$\hat{x}_t$ 是时间 $t$ 的预测状态估计值。
$P_t$ 是状态估计的协方差矩阵。
$Q$ 是过程噪声协方差矩阵。

更新步骤：根据观测方程结合新观测值来更新状态估计。

$K_t = P_t H^T (H P_t H^T + R)^{-1}$

$\hat{x}_t = \hat{x}_t + K_t (y_t - H \hat{x}_t)$

$P_t = (I - K_t H) P_t$

其中：

$K_t$ 是卡尔曼增益矩阵。
$R$ 是观测噪声协方差矩阵。
$I$ 是单位矩阵。

状态空间模型及卡尔曼滤波特别适用于处理噪声较大的系统，并能有效地进行平滑和预测。

3. 基于状态空间的时间序列分析

状态空间模型可以用于多种时间序列分析任务，如：

趋势和季节性捕捉：通过隐藏状态变量建模趋势和季节性。
平滑与预测：卡尔曼滤波提供了对噪声数据的平滑效果，并可以对未来时刻进行预测。
异常检测：通过分析残差，状态空间模型能够有效识别异常行为。

二、实战案例

我们将使用 statsmodels 库中的状态空间模型来构建并分析一个简单的时间序列数据。

1. 数据生成与可视化

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.tsa.statespace.sarimax import SARIMAX

# 生成模拟的时间序列数据
np.random.seed(42)
n = 100
time = pd.date_range(start='2020-01-01', periods=n, freq='M')
trend = 0.5 * np.arange(n)  # 线性趋势
noise = np.random.normal(0, 1, size=n)  # 噪声
observations = trend + noise

# 创建数据框
data = pd.DataFrame({'Date': time, 'Observation': observations})
data.set_index('Date', inplace=True)

# 绘制原始数据
plt.figure(figsize=(10, 6))
plt.plot(data['Observation'], label='Observed')
plt.title('Simulated Time Series with Trend')
plt.xlabel('Date')
plt.ylabel('Value')
plt.legend()
plt.grid(True)
plt.show()

程序解释：

我们生成了一个带有线性趋势和噪声的模拟时间序列数据，并绘制了原始的观测数据。这是我们要进行状态空间建模的数据。

结果输出：

2. 状态空间模型构建

# 构建状态空间模型
model = SARIMAX(data['Observation'], order=(0, 1, 1), trend='c')
results = model.fit()

# 输出模型结果摘要
print(results.summary())

程序解释：

使用 SARIMAX 函数构建状态空间模型，设置模型为 (0,1,1)(0, 1, 1)(0,1,1) 的 ARIMA 结构，这表示一阶差分和一阶移动平均模型，并包含一个常数项表示趋势。

结果输出：

SARIMAX Results                                
==============================================================================
Dep. Variable:            Observation   No. Observations:                  100
Model:               SARIMAX(0, 1, 1)   Log Likelihood                -133.143
Date:                Thu, 05 Sep 2024   AIC                            272.286
Time:                        21:27:20   BIC                            280.072
Sample:                    01-31-2020   HQIC                           275.436
                         - 04-30-2028                                         
Covariance Type:                  opg                                         
==============================================================================
                 coef    std err          z      P>|z|      [0.025      0.975]
------------------------------------------------------------------------------
intercept      0.5014      0.004    132.527      0.000       0.494       0.509
ma.L1         -0.9998      5.939     -0.168      0.866     -12.641      10.641
sigma2         0.8236      4.902      0.168      0.867      -8.784      10.432
===================================================================================
Ljung-Box (L1) (Q):                   0.07   Jarque-Bera (JB):                 0.47
Prob(Q):                              0.80   Prob(JB):                         0.79
Heteroskedasticity (H):               0.95   Skew:                            -0.15
Prob(H) (two-sided):                  0.89   Kurtosis:                         2.85
===================================================================================

3. 模型拟合与平滑

# 获取平滑后的状态
smoothed_state = results.smoothed_state

# 检查平滑状态的维度
print("Smoothed state shape:", smoothed_state.shape)

# 提取平滑的趋势状态
# 如果 smoothed_state 的维度匹配，我们提取第一个平滑状态作为趋势部分
smoothed_trend = smoothed_state[:, 0] if smoothed_state.shape[0] == len(data) else np.zeros(len(data))

程序解释：

通过卡尔曼滤波，模型对隐藏的状态进行平滑。模型输出了平滑的状态参数（即去噪后的趋势）和未来的预测值。

4. 预测与结果可视化

# 获取预测值
forecast = results.get_forecast(steps=12)
forecast_ci = forecast.conf_int()

# 绘制平滑后的状态与预测结果
plt.figure(figsize=(10, 6))
plt.plot(data.index, data['Observation'], label='Observed')
plt.plot(data.index, smoothed_trend, color='red', label='Smoothed State (Trend)')
plt.plot(forecast.predicted_mean.index, forecast.predicted_mean, color='green', label='Forecast')
plt.fill_between(forecast_ci.index, forecast_ci.iloc[:, 0], forecast_ci.iloc[:, 1], color='green', alpha=0.3)
plt.title('State Space Model: Smoothed State and Forecast')
plt.xlabel('Date')
plt.ylabel('Value')
plt.legend()
plt.grid(True)
plt.show()

程序解释：