目录
💥1 概述
📚2 运行结果
🎉3 参考文献
👨💻4 Matlab代码
💥1 概述
本文用于通过最大可能性将马尔可夫切换状态空间模型 (SSM) 拟合到多变量时间序列数据。本文考虑了三种开关SSM:开关动力学,开关观察和切换矢量自回归(VAR)。最大似然估计器通过近似电磁算法计算。(精确的计算是不可处理的,因为可能的政权历史呈指数级,M^T与M是马尔可夫链的状态/政权数,T是时间序列的长度。为了使计算易于处理,本文在EM算法的E步中使用了Kim(1994)的过滤/平滑算法。
📚2 运行结果
主函数部分代码:
%=========================================================================%
% TOY EXAMPLE FOR MATLAB TOOLBOX 'switch-ssm' %
%=========================================================================%
clc; clearvars; close all;
% This example is not calibrated to produce good statistical results,
% it is simply meant to illustrate the toolbox functionalities.
% The example demonstrates the toolbox functions for switching dynamics
% models. The functions for the switching VAR and switching observations
% models work exactly in the same way.
%-------------------------------------------------------------------------%
% 1) Simulate SSM with switching dynamics %
%-------------------------------------------------------------------------%
fprintf(['\n\n---------------\n',...
'DATA SIMULATION',...
'\n---------------\n\n']);
% Model:
% y(t) = C(S(t)) x(t) + w(t) (observation equation)
% x(t) = A(1,S(t)) x(t-1) + ... + A(p,S(t)) x(t-p) + v(t) (state equation)
% P(S(t)=j | S(t-1)=i) = Z(i,j) (regime equation)
%@@@@@ Model dimensions
N = 5; % number of variables
T = 200; % number of time points
M = 2; % number of regimes
p = 2; % VAR order
r = 2; % factor dimension
fprintf('Model: switching dynamics\nM = %d\np = %d\nr = %d\n',...
M,p,r);
%@@@@@@ Model parameters
% VAR transition matrices A(l,j), l=1:p, j=1:M
A = zeros(r,r,p,M);
A(:,:,1,1) = 0.6 * eye(r);
A(:,:,1,2) = 0.9 * eye(r);
A(:,:,2,1) = 0.1 * eye(r);
A(:,:,2,2) = -0.2 * eye(r);
% Observation matrix
[C,~,~] = svd(randn(N,r),'econ');
% Innovation variance/covariance V(v(t)|S(t)=j), j=1:M
sigQ = .01;
Q = repmat(sigQ * eye(r),[1,1,M]);
% Noise variance/covariance matrix V(w(t))
sigR = .005;
R = sigR * eye(N);
% Initial mean of state vector E(x(1)|S(1)=j), j=1:M
mu = zeros(r,M);
% Initial variance/covariance V(x(1)|S(1)=j), j=1:M
Sigma = repmat(.02 * eye(r),[1,1,M]);
% Initial probabilities P(S(1)=j), j=1:M
Pi = repelem(1/M,M,1);
% Transition probabilities P(S(t)=j | S(t-1)=i), i,j=1:M
Z = [.98,.02;.02,.98];
% Collect all model parameters in single structure
theta = struct('A',A, 'C',C, 'Q',Q, 'R',R, 'mu',mu, 'Sigma',Sigma, ...
'Pi',Pi, 'Z',Z);
%@@@@@ Simulate data
fprintf('Generating model... ')
[y,S,x] = simulate_dyn(theta,T);
% Visualize the data
tiledlayout(3,1);
nexttile
plot(x');
title("Hidden state vector")
nexttile
plot(y');
title("Observed time series")
nexttile
plot(S,'*');
title("Regimes")
xlabel("Time")
fprintf('Done')
%%
%-------------------------------------------------------------------------%
% 2) Fit switching SSM to simulated data %
%-------------------------------------------------------------------------%
fprintf(['\n\n----------------\n',...
'MODEL ESTIMATION',...
'\n----------------\n\n']);
% Minimal example (no algorithm tuning, no estimation constraints)
% [Mf,Ms,Sf,Ss,xf,xs,theta,LL] = switch_dyn(y,M,p,r);
%@@@@@ More advanced example
% Set the relative tolerance for convergence to 1e-6 (The EM will stop if
% for 5 sucessive iterations, the relative change in log-likelihood is less
% than 1e-6) and the maximum number of EM iterations to 500. Disable
% progress display.
control = struct('eps',1e-6,'ItrNo',500,'verbose',false);
% Specify EM initialization method: use fixed intervals of length 50 to
% estimate VAR parameters A and Q (one pair of estimate for each interval),
% then cluster the parameter estimates by k-means and re-estimate A(j) and
% Q(j) (j=1:M) over associated time segmentation
opts = struct('segmentation','fixed','len',25);
% Other initialization method: dichotomic segmentation (binary search for
% change points in time series). Set minimal distance between change points
% to 20
% opts = struct('segmentation','binary','delta',20);
% Set up equality constraints: force mu(j) and Sigma(j) to be equal
% across regimes (j=1:M)
equal = struct('mu',true,'Sigma',true);
% Set up fixed coefficient constraints: make covariance matrices Q
% and R diagonal (set off-diagonal terms to zero and leave other terms
% unspecified as NaN)
fixed = struct('Q',repmat(diag(NaN(r,1)),[1,1,M]),'R',diag(NaN(N,1)));
% Force the columns of the observation matrix C to have Euclidean norm 1
% and force the eigenvalues associated with VAR processes x(t) to be less
% than .98 (for stationarity)
scale = struct('A',.98,'C',1);
% Initialize EM
fprintf('Calculating MLE... ')
[thetahat0,S0] = init_dyn(y,M,p,r,opts,[],equal,fixed,scale);
% Fix mu and Sigma to their pilot estimates (these parameters are usually
% not relevant and fixing them can help the numerical stability of the EM
fixed.mu = thetahat0.mu;
fixed.Sigma = thetahat0.Sigma;
% Run EM (first pass)
[Mf,Ms,Sf,Ss,xf,xs,thetahat1,LL] = ...
switch_dyn(y,M,p,r,thetahat0,control,equal,fixed,scale); %#ok<*ASGLU>
🎉3 参考文献
[1]李强. 时滞马尔科夫切换复值神经网络的动力学分析[D].东南大学,2022.
[1]叶志勇,王泽权,唐朝君等.基于离散时间观测带有马尔可夫切换的随机网络系统的反馈控制[J].重庆理工大学学报(自然科学),2019,33(08):222-226.
点击文章左下角【阅读全文】
部分理论引用网络文献,若有侵权联系博主删除。