滤波时变扩散状态空间模型
此示例展示了如何从已知模型生成数据，将扩散状态空间模型拟合到数据中，然后过滤状态。
假设一个潜在过程包括一个AR(2)和一个MA(1)模型。有50个周期，MA(1)流程在最后25个周期退出模型。因此，前25个周期的状态方程为
在过去的25个周期里，确实如此
在哪里
和
均为高斯分布，均值为0，标准差为1。
假设序列分别从1.5和1开始，生成由50个观测值组成的随机序列
和
。
T = 50;ARMdl = arima(基于“增大化现实”技术的{0.7, -0.2},“不变”0,“方差”1);maml = arima;“马”, 0.6,“不变”0,“方差”1);X0 = [1.5 1;1.5 - 1];rng (1);x =[模拟(ARMdl,T，“Y0”x0 (: 1)),…模拟(MAMdl T / 2,“Y0”x0(: 2));南(T / 2, 1)];
模拟MA(1)数据的后25个值为南值。
潜在过程的测量使用
前25个学期，和
在过去的25个时期里
是高斯分布，均值为0，标准差为1。
使用随机潜在状态过程(x)和观测方程生成观测值。
Y = 2*sum(x'，“omitnan”)' + randn(T,1);
潜过程和观测方程共同构成一个状态空间模型。系数为未知参数，状态空间模型为
![数组$ $ \开始{}{1}& # xA;左\[{\开始数组{}{* {20}{c}} & # xA;{{间{1,t}}} \ \ & # xA;{{间{2,t}}} \ \ & # xA;{{间{3 t}}} \ \ & # xA;{{间{4 t}}} & # xA;结束\{数组}}\右]左= \[{\开始数组{}{* {20}{c}} & # xA;{{\φ_1}}和{{\φ_2}}& # 38;0 & # 38;0 \ \ & # xA; 1 & # 38; 0 & # 38; 0 & # 38; 0 \ \ & # xA; 0 & # 38; 0 & # 38; 0 &{{\θ_1}}\ \ & # xA; 0 & # 38; 0 & # 38; 0 & # 38; 0 & # xA;结束\{数组}}\右]\离开[{\开始{数组}{* {20}{c}} & # xA;{{间{1,t - 1}}} \ \ & # xA;{{间{2,t - 1}}} \ \ & # xA;{{间{3 t - 1}}} \ \ & # xA;{{间{4,t - 1}}} & # xA;结束\{数组}}\右]+左\ [数组{\开始{}{* {20}{c}} & # xA; 1 & # 38; 0 \ \ & # xA; 0 & # 38; 0 \ \ & # xA; 0 & # 38; 1 \ \ & # xA; 0 & # 38; 1 & # xA;结束\{数组}}\右]\离开[{\开始{数组}{* {20}{c}} & # xA; {{u_ {1, t}}} \ \ & # xA; {{u_ {2, t}}} & # xA;结束\{数组}}\右]\ \ & # xA; {y_t} = a({间{1,t}} +{间{3 t}}) + {\ varepsilon _t} & # xA; \{数组}$ $](http://www.tatmou.com/help/examples/econ/win64/FilterTimeVaryingDiffuseStateSpaceModelExample_eq08939596656225341196.png)
在前25个学期，
![数组$ $ \开始{}{1}& # xA;左\[{\开始数组{}{* {20}{c}} & # xA;{{间{1,t}}} \ \ & # xA;{{间{2,t}}} & # xA;结束\{数组}}\右]左= \[{\开始数组{}{* {20}{c}} & # xA;{{\φ_1}}和{{\φ_2}}& # 38;0 & # 38;0 \ \ & # xA; 1 & # 38; 0 & # 38; 0 & # 38; 0 & # xA;结束\{数组}}\右]\离开[{\开始{数组}{* {20}{c}} & # xA;{{间{1,t - 1}}} \ \ & # xA;{{间{2,t - 1}}} \ \ & # xA;{{间{3 t - 1}}} \ \ & # xA;{{间{4,t - 1}}} & # xA;结束\{数组}}\右]+ \离开[{\开始{数组}{* {20}{c}} & # xA; 1 \ \ & # xA; 0 & # xA;结束\{数组}}\右]{u_ {1, t}} \ \ & # xA; b {y_t} ={间{1,t}} +{\ varepsilon _t} & # xA; \{数组}$ $](http://www.tatmou.com/help/examples/econ/win64/FilterTimeVaryingDiffuseStateSpaceModelExample_eq14434019549381870163.png)
对于第26期，和
![数组$ $ \开始{}{1}& # xA;左\[{\开始数组{}{* {20}{c}} & # xA;{{间{1,t}}} \ \ & # xA;{{间{2,t}}} & # xA;结束\{数组}}\右]左= \[{\开始数组{}{* {20}{c}} & # xA;{{\φ_1}}和{{\φ_2}}\ \ & # xA; 1 & # 38; 0 & # xA;结束\{数组}}\右]\离开[{\开始{数组}{* {20}{c}} & # xA;{{间{1,t - 1}}} \ \ & # xA;{{间{2,t - 1}}} & # xA;结束\{数组}}\右]+ \离开[{\开始{数组}{* {20}{c}} & # xA; 1 \ \ & # xA; 0 & # xA;结束\{数组}}\右]{u_ {1, t}} \ \ & # xA; b {y_t} ={间{1,t}} + {\ varepsilon _t} & # xA; \{数组}$ $](http://www.tatmou.com/help/examples/econ/win64/FilterTimeVaryingDiffuseStateSpaceModelExample_eq14967989555778478098.png)
在过去的24个周期里。
编写一个函数，指定参数如何在参数个数映射到状态空间模型矩阵、初始状态值和状态类型。
版权所有2015 MathWorks, Inc。函数[A,B,C,D,Mean0,Cov0,StateType] = diffuseAR2MAParamMap(params,T)%diffuseAR2MAParamMap时变扩散状态空间模型参数%映射函数％这个函数将向量参数映射到状态空间矩阵(A, B，% C，和D)和状态类型(StateType)。从周期1到T/2%状态模型为AR(2)和MA(1)模型，观测模型为两种状态的和。从周期T/2 + 1到T，状态模型为%只是AR(2)模型。AR(2)模型是漫反射的。A1 = {[params(1) params(2) 0 0;1 0 0 0;0 0 0参数(3);0 0 0 0]};B1 = {[10];0 0;0 1;0 1]};C1 = {params(4)*[1 0 1 0]};Mean0 = []; Cov0 = []; StateType = [2 2 0 0]; A2 = {[params(1) params(2) 0 0; 1 0 0 0]}; B2 = {[1; 0]}; A3 = {[params(1) params(2); 1 0]}; B3 = {[1; 0]}; C3 = {params(5)*[1 0]}; A = [repmat(A1,T/2,1);A2;repmat(A3,(T-2)/2,1)]; B = [repmat(B1,T/2,1);B2;repmat(B3,(T-2)/2,1)]; C = [repmat(C1,T/2,1);repmat(C3,T/2,1)]; D = 1;结束
将此代码保存为名为diffuseAR2MAParamMap在您的MATLAB®路径。
通过传递函数创建扩散状态空间模型diffuseAR2MAParamMap的函数句柄dssm。
Mdl = dssm(@(params)diffuseAR2MAParamMap(params,T));
dssm隐式地创建扩散状态空间模型。通常，您无法验证隐式创建的扩散状态空间模型。
为了估计参数，将观察到的响应(y)估计。为未知参数指定任意一组正初始值。
Params0 = 0.1*ones(5,1);EstMdl = estimate(Mdl,y,params0);
方法:最大似然(fminc)有效样本量:48对数似然:-110.313赤池信息准则:230.626贝叶斯信息准则:240.186 | Coeff Std Err t Stat Prob --------------------------------------------------- c(1) | 0.44041 0.27687 1.59069 0.11168 c(2) | 0.03949 0.29585 0.13349 0.89380 c(3) | 0.78364 1.49223 0.52515 0.59948 c(4) | 1.64260 0.66737 2.46133 0.01384 c(5) | 1.90409 0.49374 3.85648 0.00012 | |最终状态Std Dev t Stat Prob x(1) | -0.81932 0.46706 -1.75420 0.07940 x(2) | -0.29909 0.45939 -0.65107 0.51500
EstMdl是一个dssm包含估计系数的模型。状态空间模型的似然面可能包含局部最大值。因此，尝试几个初始参数值，或者考虑使用完善。
过滤状态并通过传递获得状态预测EstMdl观察到的反应过滤器。
[~, ~,输出]=过滤器(EstMdl y);
输出是一个T-by-1结构数组，包含过滤状态和状态预测。
转换输出到一张桌子旁。
OutputTbl = struct2table(输出);OutputTbl (1:10, 1:5)%显示前5个变量的前10行
ans = 10 x5表LogLikelihood FilteredStates FilteredStatesCov ForecastedStates ForecastedStatesCov _____________ ______________ _________________ ________________ ___________________ { 0 x0双}{0 x0双}{0 x0双}{0 x0双}{0 x0双}{0 x0双}{0 x0双}{0 x0双}{0 x0双}{0 x0双}{[-2.3218]}{4 x1双}{4 x4双}{4 x1双}{4 x4双}{[-2.4464]}{4 x1双}{4 x4双}{4 x1双}{4 x4双}{[-3.8758]}{4 x1双}{4 x4双}{4 x1双}{4x4双}{[-2.5212]}{4x1双}{4x4双}{4x1双}{4x4双}{4x1双}{4x4双}{[-1.9016]}{4x1双}{4x4双}{4x1双}{4x4双}{[-1.9284]}{4x1双}{4x4双}{4x1双}{4x4双}{[-2.4110]}{4x1双}{4x4双}{4x1双}{4x4双}{[-2.6502]}{4x1双}{4x4双}{4x4双}{4x1双}{4x1双}{4x1双}{4x4双}{4x1双}{4x1双}{4x4双}{4x1双}{4x1双}{4x4双}{4x1双}{4x4双}{4x1双}{4x4双}{4x1双}{4x4双}{4x1双}{4x4双}{4x1双}{4x4双}{4x1双}{4x4双}{4x1双}{4x4双}{4x1双}{4x4双}}
表的前两行包含空单元格或零，它们对应于初始化漫反射卡尔曼滤波器所需的观测值。也就是说,SwitchTime是2。
SwitchTime = 2;
从表中提取过滤和预测的状态。回想一下，这两个不同的状态分别位于位置1和3。位置2和4中的状态有助于指定感兴趣的过程。
statidx = [1 3];%表示感兴趣的指数FilteredStates = NaN(T,numel(stateIdx));predicstedstates = NaN(T,numel(stateIdx));为t = (SwitchTime + 1): t maxInd = size(Output(t).FilteredStates,1);mask = statidx <= maxInd;FilteredStates(t,mask) = Output(t).FilteredStates(statidx (mask)，1);predicstedstates (t,mask) = Output(t).ForecastedStates(statidx (mask)，1);结束FilteredStates(1:SwitchTime，:) = 0;predicstedstates (1:SwitchTime，:) = 0;
为每个模型绘制真实状态值、过滤状态和状态预测。
图的阴谋(1:T, x (: 1),“- k”1: T, FilteredStates (: 1),“:r”，…1: T, ForecastedStates (: 1),“——g”，“线宽”2);标题(“AR(2)状态值”)包含(“时间”) ylabel (“国家价值”)({传奇“真实状态值”，“过滤状态值”，“国家预测”});图的阴谋(1:T, x (:, 2),“- k”1: T, FilteredStates (:, 2),“:r”，…1: T, ForecastedStates (:, 2),“——g”，“线宽”2);标题(“MA(1)状态值”)包含(“时间”) ylabel (“国家价值”)({传奇“真实状态值”，“过滤状态值”，“国家预测”});
