平滑时变扩散状态空间模型
这个例子展示了如何从已知模型生成数据,为数据拟合扩散状态空间模型,然后平滑状态。
假设一个潜在过程包括一个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);amdl = 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);
潜伏过程和观测方程共同构成了一个状态空间模型。如果系数为未知参数,则状态空间模型为
在前25个时段,
对于第26期,和
在过去的24个时期。
写一个函数来指定参数的输入方式参数个数
映射到状态空间模型矩阵、初始状态值和状态类型。
版权所有The MathWorks, Inc.函数[A,B,C,D,Mean0,Cov0,StateType] = diffuseAR2MAParamMap(params,T)时变扩散状态空间模型参数%映射函数%这个函数将向量参数映射到状态空间矩阵(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;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 =估计(Mdl,y,params0);
方法:最大似然(fminunc)有效样本量: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.49222 0.52515 0.59948 c(4) | 1.64260 0.66736 2.46134 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:4)显示前四个变量的前十行
ans = 10x4 table LogLikelihood SmoothedStates SmoothedStatesCov SmoothedStateDisturb _____________ ______________ _________________ ____________________ {0x0 double} {0x0 double} {0x0 double} {0x0 double} {0x0 double} {0x0 double} {0x0 double} {0x0 double} {0x0 double} {0x0 double} {0x0 double} {0x0 double} {0x0 double} {0x0 double} {0x0 double} {0x0 double} {0x0 double} {0x0 double} {[-2.3218]} {4x1 double} {4x4 double} {2x1 double} {[-2.4464]} {4x1 double} {4x4 double} {2x1 double} {[-3.8758]} {4x1 double} {4x4 double} {2x1 double} {[- 2.8758]} {4x1 double} {4x4 double} {2x1 double} {2x1 double} {[- 2.512]} {4x1 double} {4x4 double} {2x1 double} {2x1 double} {[- 2.512]} {4x1 double} {4x4 double} {2x1 double} {2x1 double} {2x1 double} {[- 2.512]} {4x1 double} {4x4 double} {2x1 double} {2x1 double} {[-1.9016]} {4x1 double} {4x4 double}Double} {2x1 Double} {[-1.9284]} {4x1 Double} {4x4 Double} {2x1 Double} {[-2.4110]} {4x1 Double} {4x4 Double} {2x1 Double} {[-2.6502]} {4x1 Double} {4x4 Double} {2x1 Double} {2x1 Double}
表的前两行包含空单元格或零,它们对应于初始化扩散卡尔曼滤波器所需的观察值。也就是说,SwitchTime
是2。
SwitchTime = 2;
从中提取平滑状态输出
,并计算他们的95%个体wald型置信区间。回想一下,位置1和3是两个不同的状态。位置2和4中的状态有助于指定感兴趣的过程。
statidx = [1 3];表示感兴趣的指数smooththedstates = nan(T,numel(stateIdx));CI = nan(T,2,numel(stateIdx));为t = (SwitchTime + 1): t maxInd = size(Output(t).SmoothedStates,1);mask = statidx <= maxInd;smooththedstates (t,mask) = Output(t).SmoothedStates(statidx (mask),1);CovX = Output(t).SmoothedStatesCov(statidx (mask), statidx (mask));CI (t): 1) = SmoothedStates (t, 1) + 1.96 * sqrt (CovX (1,1)) * (1);如果(max (stateIdx(面具))> 1)CI (t): 2) = SmoothedStates (t, 2) + 1.96 * sqrt (CovX (2, 2)) * (1);结束结束smooththedstates (1:SwitchTime,:) = 0;CI(1:SwitchTime,:,:) = 0;
绘制每个模型的真实状态值、平滑状态和95%平滑状态置信区间。
图的阴谋(1:T, x (: 1),“b”1: T, SmoothedStates (: 1),“k”1: T, CI (:,: 1),“——r”);标题(“AR(2)状态值”)包含(“时间”) ylabel (“国家价值”)({传奇“真实状态值”,平滑的状态值,95%置信区间的});图的阴谋(1:T, x (:, 2),“b”1: T, SmoothedStates (:, 2),“k”1: T, CI (:,:, 2),“——r”);标题(MA(1)状态值)包含(“时间”) ylabel (“国家价值”)({传奇“真实状态值”,平滑的状态值,95%置信区间的});