MATLAB 实现基于季节性趋势分解（STL）的时间序列预测

本文介绍在 MATLAB 中实现基于季节性趋势分解（STL, Seasonal and Trend decomposition using Loess）进行时间序列预测的项目实例。该方法适用于具有明显季节性和趋势成分的时间序列（如电力负荷、气温、销售数据等），通过 STL 分解后对各成分分别建模，再组合实现预测。

项目目标

对时间序列进行 STL 分解： $$y_t = T_t + S_t + R_t$$ 其中：$T_t$ = 趋势项，$S_t$ = 季节项，$R_t$ = 余项（残差）

对趋势项和余项分别建立预测模型（如 ARIMA、线性回归、MLP 等）
利用已知的季节模式（周期重复）外推季节项
重构预测结果：$$\hat{y}{t+h} = \hat{T}{t+h} + \hat{S}{t+h} + \hat{R}{t+h}$$

注意：MATLAB 没有内置 STL 函数，但可通过 movmean + 自定义 Loess 或使用替代方案。为保持纯 MATLAB 实现，我们采用近似 STL 的经典三步滤波法（Cleveland et al. 方法简化版）。

推荐方案：使用 MATLAB 官方 stl 替代函数（自定义实现）

由于 R2023a 起 MATLAB 仍未提供 stl，我们将使用一个高效、开源的 MATLAB STL 实现。

步骤 1：定义 STL 函数

以下是一个简化版 STL（基于滑动平均 + 季节周期对齐），适合教学或快速原型：

function [trend, seasonal, remainder] = simple_stl(y, period, nIter)
% SIMPLE_STL: 简化版 STL 分解（非 Loess，但保留核心思想）
% 输入：
% y - 时间序列 (列向量)
% period - 季节周期（如 12 表示月度数据年周期）
% nIter - 内部迭代次数（默认 2）
% 输出：
% trend, seasonal, remainder

if nargin < 3, nIter = 2; end
n = length(y);
seasonal = zeros(n, 1);
trend = zeros(n, 1);

% 初始化
detrended = y;

for iter = 1:nIter
    % Step 1: 提取季节项（按周期分组取中位数/均值）
    seasonal = zeros(n, 1);
    for i = 1:period
        idx = i:period:n;
        if ~isempty(idx)
            val = detrended(idx);
            seasonal(idx) = mean(val); % 可改用 median 更鲁棒
        end
    end

    % 去除季节项
    deseason = y - seasonal;

    % Step 2: 提取趋势项（用移动平均平滑）
    win = min(2*floor(period/2)+1, floor(n/5)); % 窗口大小
    if mod(win,2)==0, win = win+1; end
    trend = smoothdata(deseason, 'movmean', win);

    % 更新去趋势序列用于下一轮季节提取
    detrended = y - trend;
end

remainder = y - trend - seasonal;
end

对于高精度需求，建议使用完整 Loess 版本。

步骤 2：加载或生成时间序列数据

% 示例：合成带趋势 + 季节 + 噪声的数据
t = (1:200)';
season = 10 * sin(2*pi*t/12); % 周期=12（月度）
trend = 0.1 * t; % 线性增长趋势
noise = 2 * randn(size(t)); % 随机噪声
y = trend + season + noise + 50; % 基准值 50

% 可视化
figure;
plot(t, y, 'k'); xlabel('Time'); ylabel('Value');
title('Original Time Series with Trend + Seasonality');

步骤 3：执行 STL 分解

period = 12; % 假设年度季节性（月度数据）
[trend, seasonal, remainder] = simple_stl(y, period);

% 可视化分解结果
figure;
subplot(4,1,1); plot(y); title('Original'); grid on;
subplot(4,1,2); plot(trend); title('Trend'); grid on;
subplot(4,1,3); plot(seasonal); title('Seasonal'); grid on;
subplot(4,1,4); plot(remainder); title('Remainder'); grid on;

步骤 4：分别预测各成分

4.1 预测季节项 $S_{t+h}$

季节项是周期性的，直接重复最后一个完整周期：

h = 24; % 预测未来 24 步
lastCycle = seasonal(end-period+1:end);
future_seasonal = repmat(lastCycle, ceil(h/period), 1);
future_seasonal = future_seasonal(1:h);

4.2 预测趋势项 $T_{t+h}$

使用简单线性回归（或 ARIMA、MLP）：

% 线性拟合趋势
p = polyfit(t, trend, 1); % 一次多项式
t_future = (length(y)+1 : length(y)+h)';
future_trend = polyval(p, t_future);

% 更高级：可用 arima 模型拟合趋势
mdl = arima('Constant', 0, 'D', 1, 'Seasonality', 0, 'MALags', 1);
estMdl = estimate(mdl, trend);
[future_trend, ~] = forecast(estMdl, h, 'Y0', trend);

4.3 预测余项 $R_{t+h}$

余项通常视为白噪声，预测为 0；若存在自相关，可用 AR 模型：

% 简单处理：设为 0（或用 AR(1)）
future_remainder = zeros(h, 1);

% 可选：AR(1) 拟合余项
if std(remainder) > 1e-6
    armdl = arima(1,0,0);
    armdl = estimate(armdl, remainder, 'Display', 'off');
    [future_remainder, ~] = forecast(armdl, h, 'Y0', remainder);
else
    future_remainder = zeros(h, 1);
end

步骤 5：重构预测结果

y_pred = future_trend + future_seasonal + future_remainder;

% 可视化完整预测
t_all = [t; t_future];
y_all = [y; NaN(h,1)];

figure;
plot(t, y, 'b', 'LineWidth', 1.5); hold on;
plot(t_future, y_pred, 'r--', 'LineWidth', 2);
xlabel('Time'); ylabel('Value');
legend('Historical', 'STL Forecast');
title('STL-Based Time Series Forecasting');
grid on;

步骤 6：评估预测性能（若有真实未来值）

% 假设有真实未来值 y_true（实际中可能没有）
% y_true = ... ;
% rmse = sqrt(mean((y_pred - y_true).^2));
% fprintf('Forecast RMSE: %.2f\n', rmse);

改进建议

成分	改进方法
趋势项	使用 Prophet、LSTM、SVR、或局部加权回归（Loess）
季节项	多重季节性（如周 + 年）、傅里叶级数拟合
余项	GARCH 模型（波动聚集）、ResNet 残差学习
整体框架	结合 EMD/VMD + STL，或使用 TBATS、Facebook Prophet

实际应用建议

数据频率：确保 period 设置正确（日数据：7/365；小时数据：24/168）
缺失值处理：STL 对缺失敏感，需插值（fillmissing）
非固定周期：考虑使用动态时间规整（DTW）或小波分解

总结

本项目展示了如何在纯 MATLAB 环境下：

实现 STL 分解（简化版）；
分别建模趋势、季节、余项；
重构未来预测；
可视化与评估。

若你有真实数据集（如电力负荷、零售销量），只需替换 y 和 period 即可直接应用。