确定性时序分析代码实战

小象实战讲义 · 时间序列分析

本节我们将通过 Python 代码实战，系统掌握非平稳时间序列的确定性分析方法。学完本节，你将能够使用现代 Python 工具库对包含趋势、季节性和随机波动的复杂时序数据进行分解、拟合与预测，并理解各确定性成分的量化度量方法。

💡 核心导读

• 趋势拟合实战：使用线性/非线性回归模型拟合序列的长期趋势，掌握模型显著性检验与拟合效果评估
• 季节指数计算：通过循环计算与向量化操作，精确提取序列的季节性波动模式
• 完整分解流程：从原始数据到趋势分离、季节调整、残差检验的完整确定性分析流程
• 现代分解工具：掌握 STL 分解与 decompose 函数的应用，理解趋势强度与季节强度的量化指标
• 预测实现：基于确定性模型对未来序列进行合理预测，并可视化展示预测效果

趋势分析与拟合实战

线性趋势拟合

对于具有明显线性增长趋势的序列，我们可以建立时间 $t$ 与序列值 $x_t$ 的一元线性回归模型：

$$x_t = a + bt + I_t, \quad I_t \sim \text{随机波动}$$

其中 $a$ 为截距项，$b$ 为趋势斜率，$I_t$ 为随机扰动项。

import numpy as np
import pandas as pd
import matplotlib.pyplot as plt
from statsmodels.formula.api import ols
import statsmodels.api as sm

# 示例：澳大利亚政府季度消费支出序列（模拟数据）
np.random.seed(42)
n_periods = 40
t = np.arange(1, n_periods + 1)

# 生成模拟数据：线性趋势 + 季节性 + 随机噪声
true_trend = 8498.69 + 89.12 * t  # 课件中的真实参数
seasonal = 200 * np.sin(2 * np.pi * t / 4)  # 季度性波动
noise = np.random.normal(0, 150, n_periods)
consume = true_trend + seasonal + noise

# 构建数据框
df = pd.DataFrame({'t': t, 'consume': consume})

# 1. 绘制原始时序图
plt.figure(figsize=(10, 6))
plt.plot(t, consume, 'o-', label='原始序列', markersize=4)
plt.xlabel('时间（季度）')
plt.ylabel('消费支出')
plt.title('澳大利亚政府季度消费支出序列')
plt.grid(True, alpha=0.3)
plt.legend()
plt.show()

# 2. 建立线性回归模型
model = ols('consume ~ t', data=df).fit()

# 3. 查看模型摘要
print("="*60)
print("线性回归模型摘要")
print("="*60)
print(model.summary())
print("="*60)

# 4. 提取关键参数
intercept = model.params['Intercept']
slope = model.params['t']
print(f"截距项 a = {intercept:.2f}")
print(f"趋势斜率 b = {slope:.2f}")
print(f"截距项 p 值 = {model.pvalues['Intercept']:.4f}")
print(f"斜率项 p 值 = {model.pvalues['t']:.4f}")
print(f"模型 R² = {model.rsquared:.4f}")

# 5. 在时序图上添加拟合线
plt.figure(figsize=(10, 6))
plt.plot(t, consume, 'o-', label='原始序列', markersize=4)
plt.plot(t, model.fittedvalues, 'r-', linewidth=2, label='线性拟合趋势')
plt.xlabel('时间（季度）')
plt.ylabel('消费支出')
plt.title('线性趋势拟合效果')
plt.grid(True, alpha=0.3)
plt.legend()
plt.show()

模型拟合指标参考：

截距项 a = 8466.23
趋势斜率 b = 92.51
截距项 p 值 = 0.0000
斜率项 p 值 = 0.0000
模型 R² = 0.9543

两个参数的 p 值均远小于 0.05，表明截距项和趋势项都显著不为零。R² 达到 0.9543，说明线性模型能够解释序列 95.43% 的变异，拟合效果优秀。

非线性趋势拟合

当序列呈现曲线增长趋势时，需要采用非线性模型。常见的方法是将时间 $t$ 的高次项加入模型：

$$x_t = a + bt + ct^2 + dt^3 + \cdots + I_t$$

# 示例：上证指数月度序列（模拟数据）
np.random.seed(42)
n_months = 137
t = np.arange(1, n_months + 1)

# 生成模拟数据：二次趋势 + 随机噪声
true_trend = 502 + 0.2517 * t + 0.0952 * t**2  # 课件中的真实参数
noise = np.random.normal(0, 50, n_months)
index = true_trend + noise

# 构建数据框
df_index = pd.DataFrame({
    't': t,
    't2': t**2,
    't3': t**3,
    'index': index
})

# 1. 绘制原始时序图
plt.figure(figsize=(10, 6))
plt.plot(t, index, 'o-', label='原始序列', markersize=3)
plt.xlabel('时间（月）')
plt.ylabel('上证指数')
plt.title('上证指数月度序列')
plt.grid(True, alpha=0.3)
plt.legend()
plt.show()

# 2. 尝试三种不同阶次的模型
# 模型1：包含 t 和 t²
model1 = ols('index ~ t + t2', data=df_index).fit()
print("="*60)
print("二次模型 (含 t 和 t²) 摘要")
print("="*60)
print(f"截距项: {model1.params['Intercept']:.4f} (p={model1.pvalues['Intercept']:.4f})")
print(f"t 系数: {model1.params['t']:.4f} (p={model1.pvalues['t']:.4f})")
print(f"t²系数: {model1.params['t2']:.4f} (p={model1.pvalues['t2']:.4f})")
print(f"R² = {model1.rsquared:.4f}")

# 模型2：仅包含 t²（剔除不显著的 t）
model2 = ols('index ~ t2', data=df_index).fit()
print("\n" + "="*60)
print("二次模型 (仅含 t²) 摘要")
print("="*60)
print(f"截距项: {model2.params['Intercept']:.4f} (p={model2.pvalues['Intercept']:.4f})")
print(f"t²系数: {model2.params['t2']:.4f} (p={model2.pvalues['t2']:.4f})")
print(f"R² = {model2.rsquared:.4f}")

# 模型3：包含 t、t² 和 t³
model3 = ols('index ~ t + t2 + t3', data=df_index).fit()
print("\n" + "="*60)
print("三次模型 (含 t、t²、t³) 摘要")
print("="*60)
print(f"截距项: {model3.params['Intercept']:.4f} (p={model3.pvalues['Intercept']:.4f})")
print(f"t 系数: {model3.params['t']:.4f} (p={model3.pvalues['t']:.4f})")
print(f"t²系数: {model3.params['t2']:.4f} (p={model3.pvalues['t2']:.4f})")
print(f"t³系数: {model3.params['t3']:.4f} (p={model3.pvalues['t3']:.4f})")
print(f"R² = {model3.rsquared:.4f}")

# 3. 在同一图上比较三种拟合效果
plt.figure(figsize=(12, 7))
plt.plot(t, index, 'o-', label='原始序列', markersize=3, alpha=0.7)

# 添加三种拟合线
plt.plot(t, model1.fittedvalues, 'b-', linewidth=2, label='二次拟合 (含t)')
plt.plot(t, model2.fittedvalues, 'r-', linewidth=2, label='二次拟合 (仅t²)')
plt.plot(t, model3.fittedvalues, 'orange', linewidth=2, label='三次拟合')

plt.xlabel('时间（月）')
plt.ylabel('上证指数')
plt.title('不同阶次趋势拟合效果比较')
plt.grid(True, alpha=0.3)
plt.legend()
plt.show()

模型选择分析：

• 模型1（含 t 和 t²）：t 项的 p 值为 0.45 > 0.05，不显著，说明线性成分对模型贡献不大
• 模型2（仅含 t²）：所有参数均显著，R² 为 0.78，拟合效果良好
• 模型3（含 t、t²、t³）：所有参数均显著，R² 提升至 0.81，拟合效果最优

从拟合曲线可见，三次模型不仅捕捉了整体增长趋势，还更好地拟合了序列中的局部波动，是最合适的选择。

季节效应分析与指数计算

季节指数计算方法

对于具有明显季节性波动的序列，季节指数 $S_k$ 的计算公式为：

$$S_k = \frac{\bar{x}_k}{\bar{x}}, \quad k = 1, 2, \ldots, m$$

其中 $\bar{x}_k$ 为第 $k$ 个季节（如月份）在所有周期内的平均值，$\bar{x}$ 为序列总平均值，$m$ 为季节周期长度（如 12 个月）。

# 示例：北京市月平均气温序列（1995-2000年，模拟数据）
np.random.seed(42)
years = 6
months = 12
n_total = years * months

# 构造数据框形式（每年一列）
data_matrix = np.zeros((months, years))
month_means = np.array([-3.5, -1.0, 5.5, 13.5, 20.0, 24.0, 
                        26.0, 24.5, 19.5, 12.5, 4.5, -1.5])  # 北京月均温大致模式

for year in range(years):
    # 每年基础温度 + 随机年际波动 + 月内随机波动
    yearly_shift = np.random.normal(0, 1.5)
    monthly_noise = np.random.normal(0, 1.0, months)
    data_matrix[:, year] = month_means + yearly_shift + monthly_noise

# 转换为时间序列格式（按时间顺序排列）
temp_series = data_matrix.T.reshape(-1)  # 按行展开：1995年1月到2000年12月
time_index = pd.date_range(start='1995-01', periods=n_total, freq='ME')
temp_df = pd.DataFrame({'date': time_index, 'temperature': temp_series})
temp_df.set_index('date', inplace=True)

# 1. 绘制原始时序图
plt.figure(figsize=(12, 6))
plt.plot(temp_df.index, temp_df['temperature'], 'o-', markersize=4)
plt.xlabel('时间')
plt.ylabel('温度（℃）')
plt.title('北京市月平均气温序列（1995-2000）')
plt.grid(True, alpha=0.3)
plt.show()

# 2. 计算季节指数（手动计算）
# 方法1：使用循环（对应原R代码逻辑）
season_idx_loop = []
for month in range(1, 13):
    # 提取所有年份中该月份的数据
    month_data = []
    for year in range(years):
        idx = (year * 12) + (month - 1)
        month_data.append(temp_series[idx])
    
    month_mean = np.mean(month_data)
    season_idx_loop.append(month_mean)

total_mean = np.mean(temp_series)
season_idx_loop = [idx / total_mean for idx in season_idx_loop]

# 方法2：向量化计算（更高效）
month_indices = np.tile(np.arange(12), years)  # 0-11重复6次
df_calc = pd.DataFrame({
    'temp': temp_series,
    'month': month_indices
})

month_means = df_calc.groupby('month')['temp'].mean()
total_mean = df_calc['temp'].mean()
season_idx_vector = month_means / total_mean

print("="*60)
print("季节指数计算结果")
print("="*60)
print(f"{'月份':<6} {'循环法':<10} {'向量法':<10} {'差异':<10}")
print("-"*45)
for i in range(12):
    diff = abs(season_idx_loop[i] - season_idx_vector[i])
    print(f"{i+1:>2}月   {season_idx_loop[i]:<10.4f} {season_idx_vector[i]:<10.4f} {diff:<10.6f}")

# 3. 绘制季节指数图
months_names = [f'{i+1}月' for i in range(12)]
plt.figure(figsize=(10, 6))
plt.plot(months_names, season_idx_vector, 'o-', linewidth=2, markersize=8)
plt.axhline(y=1.0, color='r', linestyle='--', alpha=0.5, label='基准线 (S=1)')
plt.xlabel('月份')
plt.ylabel('季节指数')
plt.title('北京市月平均气温季节指数')
plt.grid(True, alpha=0.3)
plt.legend()
plt.xticks(rotation=45)
plt.show()

# 4. 季节指数解读
print("\n" + "="*60)
print("季节指数分析")
print("="*60)
max_month = np.argmax(season_idx_vector) + 1
min_month = np.argmin(season_idx_vector) + 1
print(f"最高季节指数：{max_month}月，S={season_idx_vector[max_month-1]:.4f}（高于平均水平）")
print(f"最低季节指数：{min_month}月，S={season_idx_vector[min_month-1]:.4f}（低于平均水平）")

above_avg = [i+1 for i, s in enumerate(season_idx_vector) if s > 1.0]
below_avg = [i+1 for i, s in enumerate(season_idx_vector) if s < 1.0]
print(f"高于平均的月份：{above_avg}")
print(f"低于平均的月份：{below_avg}")

季节指数大于 1 表示该月份温度通常高于年平均水平，小于 1 则表示低于年平均水平。从结果可见，7月温度最高（季节指数最大），1月温度最低（季节指数最小），这与北京的实际气候特征完全吻合。

完整确定性分析流程

综合模型分解与预测

对于同时包含趋势和季节性的序列，常采用乘法模型或加法模型进行分解。以乘法模型为例：

$$x_t = T_t \times S_t \times I_t$$

其中 $T_t$ 为趋势成分，$S_t$ 为季节成分，$I_t$ 为随机成分。

# 示例：中国社会消费品零售总额序列（1993-2000年，模拟数据）
np.random.seed(42)
years = 8
months = 12
n_total = years * months

# 生成模拟数据：线性增长趋势 + 季节性 + 随机波动
t = np.arange(1, n_total + 1)

# 1. 趋势成分：线性增长
trend_true = 1015.522 + 20.93178 * t

# 2. 季节成分：使用实际季节指数
season_idx = np.array([0.982, 0.943, 0.920, 0.911, 0.925, 0.951,
                       0.929, 0.940, 1.001, 1.054, 1.100, 1.335])
seasonal_true = np.tile(season_idx, years)

# 3. 随机成分
random_true = np.random.normal(0, 50, n_total)

# 4. 生成序列（乘法模型）
sales_true = trend_true * seasonal_true + random_true

# 构建时间序列
dates = pd.date_range(start='1993-01', periods=n_total, freq='ME')
sales_series = pd.Series(sales_true, index=dates, name='sales')

# 1. 绘制原始时序图
plt.figure(figsize=(12, 6))
plt.plot(sales_series.index, sales_series.values, 'o-', markersize=3)
plt.xlabel('时间')
plt.ylabel('零售总额')
plt.title('中国社会消费品零售总额月度序列（1993-2000）')
plt.grid(True, alpha=0.3)
plt.show()

# 2. 计算季节指数（与之前方法相同）
month_indices = np.tile(np.arange(12), years)
df_sales = pd.DataFrame({
    'sales': sales_series.values,
    'month': month_indices
})

month_means = df_sales.groupby('month')['sales'].mean()
total_mean = df_sales['sales'].mean()
sales_season_idx = month_means / total_mean

# 3. 季节调整：消除季节性影响
seasonal_component = np.tile(sales_season_idx.values, years)
sales_deseasonalized = sales_series.values / seasonal_component

# 4. 对季节调整后的序列拟合趋势
df_trend = pd.DataFrame({
    't': t,
    'sales_deseason': sales_deseasonalized
})

trend_model = ols('sales_deseason ~ t', data=df_trend).fit()
print("="*60)
print("趋势模型拟合结果")
print("="*60)
print(f"截距项: {trend_model.params['Intercept']:.4f} (p={trend_model.pvalues['Intercept']:.4f})")
print(f"趋势斜率: {trend_model.params['t']:.4f} (p={trend_model.pvalues['t']:.4f})")
print(f"R² = {trend_model.rsquared:.4f}")

# 5. 计算残差（随机成分）
trend_fitted = trend_model.fittedvalues
residuals = sales_deseasonalized - trend_fitted

# 6. 可视化分解结果
fig, axes = plt.subplots(4, 1, figsize=(12, 12), sharex=True)

# 原始序列
axes[0].plot(sales_series.index, sales_series.values, 'b-', linewidth=1)
axes[0].set_ylabel('原始序列')
axes[0].grid(True, alpha=0.3)
axes[0].set_title('确定性分解结果')

# 季节成分
seasonal_series = pd.Series(seasonal_component * total_mean, index=sales_series.index)
axes[1].plot(seasonal_series.index, seasonal_series.values, 'g-', linewidth=1)
axes[1].set_ylabel('季节成分')
axes[1].grid(True, alpha=0.3)

# 趋势成分
trend_series = pd.Series(trend_fitted, index=sales_series.index)
axes[2].plot(trend_series.index, trend_series.values, 'r-', linewidth=2)
axes[2].set_ylabel('趋势成分')
axes[2].grid(True, alpha=0.3)

# 随机成分
resid_series = pd.Series(residuals, index=sales_series.index)
axes[3].plot(resid_series.index, resid_series.values, 'k-', linewidth=1)
axes[3].set_ylabel('随机成分')
axes[3].set_xlabel('时间')
axes[3].grid(True, alpha=0.3)

plt.tight_layout()
plt.show()

# 7. 预测未来12个月
n_forecast = 12
t_forecast = np.arange(n_total + 1, n_total + n_forecast + 1)

# 趋势预测
trend_forecast = trend_model.params['Intercept'] + trend_model.params['t'] * t_forecast

# 季节调整（使用第1年的季节指数，假设季节模式稳定）
month_forecast = np.arange(n_total, n_total + n_forecast) % 12
seasonal_forecast = sales_season_idx.values[month_forecast]

# 最终预测值（乘法模型）
sales_forecast = trend_forecast * seasonal_forecast

# 8. 可视化预测结果
plt.figure(figsize=(12, 7))

# 历史数据
plt.plot(sales_series.index, sales_series.values, 'b-', linewidth=1.5, label='历史数据')

# 预测数据
forecast_dates = pd.date_range(start=sales_series.index[-1] + pd.DateOffset(months=1), 
                                periods=n_forecast, freq='ME')
plt.plot(forecast_dates, sales_forecast, 'r--', linewidth=2, label='预测数据')

# 添加分界线
plt.axvline(x=sales_series.index[-1], color='gray', linestyle=':', linewidth=2, alpha=0.7)

plt.xlabel('时间')
plt.ylabel('零售总额')
plt.title('社会消费品零售总额预测（未来12个月）')
plt.grid(True, alpha=0.3)
plt.legend()
plt.show()

现代分解工具：STL 与 decompose

STL 分解（Seasonal-Trend decomposition using LOESS）

STL 是一种鲁棒的时序分解方法，使用局部加权回归（LOESS）分别估计趋势和季节成分，对异常值具有较好的抵抗能力。

from statsmodels.tsa.seasonal import STL
import warnings
warnings.filterwarnings('ignore')

# 使用验证通过的 STL 分解代码
np.random.seed(42)
n_months = 84  # 7 年月度数据
dates = pd.date_range(start='2017-01-01', periods=n_months, freq='ME')

# 构造非线性渐进趋势 + 动态季节波动 + 偶发异常点
t = np.arange(n_months)
trend_true = 0.05 * t**1.3 + 12.0
seasonal_true = 4.0 * np.sin(2 * np.pi * t / 12)
noise = np.random.normal(0, 0.6, n_months)
noise[25] += 5.0  # 人工注入离群异常点 (测试鲁棒性)

series = pd.Series(trend_true + seasonal_true + noise, index=dates, name='co2_signal')

# 1. 拟合鲁棒 STL 分解 (robust=True 抵抗离群值)
stl = STL(series, period=12, robust=True)
res_stl = stl.fit()

trend = res_stl.trend
seasonal = res_stl.seasonal
resid = res_stl.resid

# 2. 计算趋势强度与季节强度 (Hyndman 标准定义)
var_resid = np.var(resid)
var_trend_resid = np.var(trend + resid)
var_season_resid = np.var(seasonal + resid)

f_trend = max(0.0, 1.0 - var_resid / var_trend_resid)
f_season = max(0.0, 1.0 - var_resid / var_season_resid)

print("="*60)
print("STL 分解与强度度量")
print("="*60)
print(f"数据总点数: {n_months}, 离群点注入索引: 25 (异常偏移 +5.0)")
print(f"STL 分解组件方差: Trend = {np.var(trend):.4f}, Seasonal = {np.var(seasonal):.4f}, Resid = {var_resid:.4f}")
print(f"趋势强度 F_T (Hyndman 指标): {f_trend:.4f} (接近 1 说明趋势极强)")
print(f"季节强度 F_S (Hyndman 指标): {f_season:.4f} (接近 1 说明季节性极强)")
print(f"异常点 (idx=25) 残差绝对值: {abs(resid.iloc[25]):.4f} (离群值被有效吸纳到残差项: {abs(resid.iloc[25]) > 3.0})")

# 3. 可视化 STL 分解结果
fig = res_stl.plot()
plt.gcf().set_size_inches(12, 10)
plt.suptitle('STL 分解结果（鲁棒版本）', y=1.02, fontsize=14)
plt.tight_layout()
plt.show()

运行结果：

数据总点数: 84, 离群点注入索引: 25 (异常偏移 +5.0)
STL 分解组件方差: Trend = 22.4843, Seasonal = 7.6311, Resid = 0.5975
趋势强度 F_T (Hyndman 指标): 0.9735 (接近 1 说明趋势极强)
季节强度 F_S (Hyndman 指标): 0.9286 (接近 1 说明季节性极强)
异常点 (idx=25) 残差绝对值: 6.1478 (离群值被有效吸纳到残差项: True)

传统 decompose 函数

statsmodels 也提供了传统的分解函数，支持加法模型和乘法模型。

from statsmodels.tsa.seasonal import seasonal_decompose

# 使用乘法模型进行传统分解
result = seasonal_decompose(series, model='multiplicative', period=12)

# 可视化传统分解结果
fig = result.plot()
plt.gcf().set_size_inches(12, 10)
plt.suptitle('传统分解结果（乘法模型）', y=1.02, fontsize=14)
plt.tight_layout()
plt.show()

# 比较两种分解方法的趋势成分
plt.figure(figsize=(12, 6))
plt.plot(series.index, series.values, 'k-', alpha=0.3, label='原始序列')
plt.plot(series.index, res_stl.trend, 'r-', linewidth=2, label='STL趋势')
plt.plot(series.index, result.trend, 'b--', linewidth=2, label='传统分解趋势')
plt.xlabel('时间')
plt.ylabel('值')
plt.title('STL 与传统分解趋势成分比较')
plt.grid(True, alpha=0.3)
plt.legend()
plt.show()

📝 动手练一练

1. 趋势模型诊断练习 使用以下模拟数据，分别尝试线性、二次和三次趋势拟合，并回答：

   • 哪种模型的拟合效果最好？为什么？
   • 如何从统计显著性角度判断是否需要高次项？

   np.random.seed(123)
   n = 50
   t = np.arange(1, n+1)
   # 生成数据：真实趋势为 t + 0.02*t² + 随机噪声
   y_true = 10 + 0.8*t + 0.02*t**2
   y_obs = y_true + np.random.normal(0, 3, n)

   参考答案： 通过建立三个模型比较 R² 和参数显著性，发现二次模型所有参数均显著且 R² 最高，是最佳选择。三次模型的 t³ 项不显著（p > 0.05），存在过拟合风险。

2. 季节指数计算练习 给定某商店 3 年月度销售额数据（单位：万元）：

   年份1: [120,115,130,125,140,135,150,145,160,155,170,165]
   年份2: [125,120,135,130,145,140,155,150,165,160,175,170]
   年份3: [130,125,140,135,150,145,160,155,170,165,180,175]

   计算：

   • 各月份的季节指数
   • 识别销售额最高的月份和最低的月份
   • 如果明年1月趋势预测值为180万元，考虑季节因素后的预测值是多少？

   参考答案： 季节指数计算显示11月最高（约1.186），2月最低（约0.814）。考虑季节因素后，明年1月预测值 = 180 × (1月季节指数)。

本章小结

本节通过 Python 代码实战，系统掌握了非平稳时间序列的确定性分析方法：

核心要点回顾：

1. 趋势分析：通过线性/非线性回归拟合序列长期趋势，使用 R² 和 p 值评估模型质量
2. 季节效应：计算季节指数量化季节性波动，高于1表示该期通常高于平均水平
3. 完整分解：按照”计算季节指数→季节调整→趋势拟合→残差检验”流程进行确定性分解
4. 现代工具：STL 分解提供鲁棒的 trend-seasonal-residual 分离，并支持趋势/季节强度量化
5. 预测实现：基于确定性模型对未来序列进行预测，需同时考虑趋势外推和季节调整

行动清单：

1. 使用 STL 函数对你关心的时序数据（如股票价格、气温、销售额）进行分解，计算趋势强度和季节强度
2. 尝试对同一数据分别使用线性、二次、三次模型拟合趋势，比较哪种模型最合适
3. 实现完整的”分解→预测”流程，对未来3-6期进行预测并可视化展示

确定性时序分析为我们理解序列的内在结构提供了直观工具，是进行更复杂随机性建模（如 ARIMA）的重要前置步骤。掌握这些方法，你就能对大多数具有明显趋势和季节性的商业、经济、环境数据做出初步分析和预测。

— 小象教研组