【MATLAB第117期】#源码分享 | 基于MATLAB的SSM状态空间模型多元时间序列预测方法（多输入单输出）

news2025/4/23 6:42:23

【MATLAB第117期】#源码分享 | 基于MATLAB的SSM状态空间模型多元时间序列预测方法（多输入单输出）

引言

本文使用状态空间模型实现失业率递归预测，状态空间模型（State Space Model, SSM）是一种用于描述动态系统行为的数学模型，通过状态变量、输入和输出的关系来刻画系统的时变特性。
在这里插入图片描述

本示例演示如何使用Econometrics Toolbox中的状态空间模型，实现美国年度失业率的滚动窗口预测

% 模型结构：
%   x(t) = A*x(t-1) + B*u(t) + w(t)  (状态方程)
%   y(t) = C*x(t)   + D*u(t) + v(t)  (观测方程)
% 其中：
%   y(t) - 观测失业率
%   x(t) - 潜在状态变量
%   u(t) - 名义GNP增长率（外生变量）
%   w(t) ~ N(0,Q) - 过程噪声
%   v(t) ~ N(0,R) - 观测噪声

一、关键思路

‌1、滚动窗口机制‌
使用31年的窗口进行参数估计，逐年滚动预测下一年失业率
窗口滑动方式：1940-1970年间共生成31个预测窗口
‌2、状态空间模型结构‌

x(t)=[α1;01]∗x(t−1)+[β;0]∗ΔlnGNP(t)+w(t)y(t)=[10]∗x(t)+γ∗ΔlnGNP(t)+v(t)
包含AR(2)动态和GNP的外生影响
参数估计目标：α(自回归系数), β(GNP对状态的系数), γ(GNP对观测的直接系数)
3、‌数据预处理细节‌

对原始GNP数据取对数差分处理：ΔlnGNP = diff(log(GNP))
失业率使用一阶差分：ΔUR = diff(UR)

二、核心代码

%% 1. 导入年度经济数据
load Data_NelsonPlosser; % 加载Nelson-Plosser宏观经济数据集
% 数据集包含美国1909-1970年间的经济指标，包括：
%   GNPN - 名义国民生产总值 
%   UR   - 失业率
%   其他指标：CPI、工资率等

%% 2. 数据预处理
isNaN = any(ismissing(DataTable),2); % 标记包含缺失值的行
Z = DataTable.GNPN(~isNaN);  % 提取有效期的名义GNP数据（61×1）
y = DataTable.UR(~isNaN);    % 提取有效期的失业率数据（61×1）

%% 3. 创建时间序列数组
WindowSize = 31; % 滚动窗口大小（31年）
ForecastPeriod = numel(y) - WindowSize + 1; % 预测期数（61-31+1=31）

% 初始化存储矩阵
ZZ = zeros(ForecastPeriod, WindowSize); % GNP窗口数据（31×31）
yy = zeros(ForecastPeriod, WindowSize); % 失业率窗口数据（31×31）

% 创建滚动窗口数据集
m = 1;
for nYear = 1:ForecastPeriod
    ZZ(nYear,:) = Z(m:m+WindowSize-1); % 当前窗口的GNP数据
    yy(nYear,:) = y(m:m+WindowSize-1); % 当前窗口的失业率数据
    m = m + 1;
end

% 提取时间戳（1940-1970年）
Time = str2double(DataTable.Properties.RowNames(~isNaN));
Time = Time(end-ForecastPeriod+1:end); 

% 构建带时间戳的数组
ObsUnemployRate = [Time, yy]; % 失业率时间序列（31×32）
nGNP = [Time, ZZ];           % GNP时间序列（31×32）

%% 4. MATLAB状态空间模型递归估计
eUR = zeros(numel(Time),1);    % 存储预测失业率
param0 = [0.5; 0.1; -20];      % 初始参数 [A; B; D]

% 滚动窗口参数估计循环
for t = 0:numel(Time)-1
    % 数据准备
    dlZ = diff(log(ZZ(t+1,:)))';  % GNP对数差分（增长率）
    dy = diff(yy(t+1,:))';        % 失业率差分
    
    % 状态空间模型定义（调用rwAR2ParamMap函数）
    Mdl = ssm(@(c)rwAR2ParamMap(c,dy,dlZ));
    
    % 参数估计（最大似然估计）
    [Mdl, param0] = estimate(Mdl, dy, param0, 'Display', 'off');
    
    % 1步超前预测
    dyhat = forecast(Mdl, 1, dy,...
        'Predictors0', dlZ,...
        'PredictorsF', dlZ(end),...
        'Beta', param0(end));
    
    % 整合预测结果
    eUR(t+1) = dyhat + yy(t+1,end); % 预测值 = 差分预测 + 当前水平值
end

%% 5. 可视化MATLAB预测结果
figure;
axH = axes;
plot(axH, Time, y(end-numel(Time)+1:end),... % 实际失业率
    'Color', [0.9290 0.6940 0.1250], 'LineWidth', 1.2);
hold on; grid on;
plot(axH, Time, eUR,... % 预测失业率
    'Color', [0 0.4470 0.7410], 'LineWidth', 1.2);
axH.XLim(1) = Time(1);
axH.Color = [0.5020 0.5020 0.5020];
title('失业率预测对比 (%)');
legend(["实际值", "预测值"]);