用蒙特卡罗方法预测状态空间模型

打开实时脚本

这个例子展示了如何使用蒙特卡罗方法预测一个状态空间模型，并将蒙特卡罗预测与理论预测进行比较。

假设失业率的变化( $x_{1 ， t}$ )和名义国民生产总值增长率( $x_{3. ， t}$ )可以用下面的状态空间模型形式表示。

$［ \begin{array}{cccccccccccccccccccc} x_{1 ， t} \\ x_{2 ， t} \\ x_{3. ， t} \\ x_{4 ， t} \end{array} ］＝［ \begin{array}{cccccccccccccccccccc} ϕ_{1} & θ_{1} & γ_{1} & 0 \\ 0 & 0 & 0 & 0 \\ γ_{2} & 0 & ϕ_{2} & θ_{2} \\ 0 & 0 & 0 & 0 \end{array} ］［ \begin{array}{cccccccccccccccccccc} x_{1 ， t - 1} \\ x_{2 ， t - 1} \\ x_{3. ， t - 1} \\ x_{4 ， t - 1} \end{array} ］ + ［ \begin{array}{cccccccccccccccccccc} 1 & 0 \\ 1 & 0 \\ 0 & 1 \\ 0 & 1 \end{array} ］［ \begin{array}{cccccccccccccccccccc} u_{1 ， t} \\ u_{2 ， t} \end{array} ］$

$［ \begin{array}{cccccccccccccccccccc} y_{1 ， t} \\ y_{2 ， t} \end{array} ］＝［ \begin{array}{cccccccccccccccccccc} 1 & 0 & 0 & 0 \\ 0 & 0 & 1 & 0 \end{array} ］［ \begin{array}{cccccccccccccccccccc} x_{1 ， t} \\ x_{2 ， t} \\ x_{3. ， t} \\ x_{4 ， t} \end{array} ］ + ［ \begin{array}{cccccccccccccccccccc} σ_{1} & 0 \\ 0 & σ_{2} \end{array} ］［ \begin{array}{cccccccccccccccccccc} ε_{1 ， t} \\ ε_{2 ， t} \end{array} ］，$

地点:

$x_{1 ， t}$ 是当时失业率的变化吗t．
$x_{2 ， t}$ 是MA(1)效应的虚拟状态 $x_{1 ， t}$ ．
$x_{3. ， t}$ 当时的国民生产总值增长率是多少t．
$x_{4 ， t}$ 是MA(1)效应的虚拟状态 $x_{3. ， t}$ ．
$y_{1 ， t}$ 是观察到的失业率变化。
$y_{2 ， t}$ 为观察到的nGNP增长率。
$u_{1 ， t}$ 而且 $u_{2 ， t}$ 为均值为0，标准差为1的状态扰动的高斯序列。
$ε_{1 ， t}$ 观测创新的高斯级数是否具有均值0和标准差 $σ_{1}$ ．
$ε_{2 ， t}$ 观测创新的高斯级数是否具有均值0和标准差 $σ_{2}$ ．

加载Nelson-Plosser数据集，其中包含失业率和nGNP系列。

负载Data_NelsonPlosser

对数据进行预处理，取nGNP系列的自然对数，以及每个的第一差值。另外，去掉起始部分南每个系列的值。

isNaN = any(ismissing(DataTable)，2);包含nan的标记句点gnpn = dattable . gnpn (~isNaN);u = dattable . ur (~isNaN);T = size(gnpn,1);%样本量y = 0 (t -1,2);% PreallocateY (:，1) = diff(u);Y (:，2) = diff(log(gnpn));

本例继续使用无序列南值。然而，使用卡尔曼滤波器框架，该软件可以容纳包含缺失值的系列。

为了确定模型对观测值的预测效果，可以去掉最后10个观测值进行比较。

numPeriods = 10;%预测范围isY = y(1:end-numPeriods，:);样本内观测值%oosY = y(end- numperiods +1:end，:);%样本外观测值

指定系数矩阵。

A = [NaN NaN 0;0 0 0 0;NaN 0 NaN NaN;0 0 0 0];B = [10 0;1 0;0 1;0 1];C = [1 0 0 0;0 0 10 0]; D = [NaN 0; 0 NaN];

使用指定状态空间模型舰导弹．验证模型规范是否与状态空间模型一致。

Mdl = ssm(A,B,C,D)

Mdl =状态空间模型类型:ssm状态向量长度:4观测向量长度:2状态扰动向量长度:2观测创新向量长度:2模型支持的样本量:无限估计未知参数:8状态变量:x1, x2，…金宝app状态扰动:u1, u2，…观测序列:y1, y2，…观测创新:e1, e2，…未知参数:c1, c2，…状态方程:x1(t) = (c1)x1(t-1) + (c3)x2(t-1) + (c4)x3(t-1) + u1(t) x2(t) = u1(t) x3(t) = (c2)x1(t-1) + (c5)x3(t-1) + (c6)x4(t) + u2(t) x4(t) = u2(t)观测方程:y1(t) = x1(t) + (c7)e1(t) y2(t) = x3(t) + (c8)e2(t)初态分布:未指定初态均值。初始状态协方差矩阵未指定。未指定状态类型。

估计模型参数，并使用一组随机的初始参数值进行优化。限制估计 $σ_{1}$ 而且 $σ_{2}$ 对所有正数，实数使用“磅”名称-值对参数。为了数值稳定性，在软件计算参数协方差矩阵时指定Hessian，使用“CovMethod”名称-值对参数。

rng (1);Params0 = rand(8,1);[EstMdl,estParams] =估计(Mdl,isY,params0，.．.“磅”，[-Inf -Inf -Inf -Inf -Inf -Inf -Inf 0 0]，“CovMethod”，“海赛”）;

方法:最大似然(fmincon)样本量:51对数似然:-170.92赤井信息准则:357.84贝叶斯信息准则:373.295 | t统计概率多项式系数Std犯错  ---------------------------------------------------- c (1) c(2) | | 0.06750 0.16548 0.40791 0.68334 -0.01372 0.05887 -0.23302 0.81575摄氏度(3)0 c(4) | | 2.71201 0.27039 10.03006 0.83815 2.84585 0.29452 0.76836摄氏度(5)c(6) | | 0.06273 2.83473 0.02213 0.98235 0.05197 2.56875 0.02023 0.98386摄氏度(7)c(8) | | 0.00273 2.40763 0.00113 0.99910 0.00016 0.13942 0.00113 0.99910 | |最终状态性病Dev t统计概率(1)| -0.00000 0.00273 -0.00033 0.99973 x (2) | 0.12237 - 0.929540.13164 0.89527 x(3) | 0.04049 0.00016 256.76330 0 x(4) | 0.01183 0.00016 72.51366 0

EstMdl是一个舰导弹模型，您可以使用点表示法访问它的属性。

过滤估计的状态空间模型，并从最终周期中提取过滤的状态及其方差。

[~，~，Output] = filter(EstMdl,isY);

修改估计的状态空间模型，使初始状态均值和协方差是过滤后的状态及其最终周期的协方差。这建立了预测范围内的模拟。

EstMdl1 = EstMdl;EstMdl1。Mean0 = Output(end).FilteredStates;EstMdl1。Cov0 = Output(end).FilteredStatesCov;

模拟5 e5从拟合的状态空间模型观测的路径EstMdl．指定模拟每个周期的观测值。

numPaths = 5e5;simmy =模拟(EstMdl1,10，“NumPaths”, numPaths);

SimY是一个10——- - - - - -2——- - - - - -numPaths包含模拟观测的数组。一排排的SimY对应于周期，列对应于模型中的观察结果，页对应于路径。

估计预测观测值及其在预测视界内的95%置信区间。

MCFY = mean(SimY,3);CIFY =分位数(SimY，[0.025 0.975]，3);

估计理论预测波段。

[Y,YMSE] = forecast(EstMdl,10,isY);Lb = Y -根号(YMSE)*1.96;Ub = Y +根号(YMSE)*1.96;

将预测观测值与它们的真实值和预测间隔绘制出来。

figure h = plot(dates(end- numperiods -9:end)，[isY(end-9:end,1);oosY(:，1)]，“- k”，.．.日期(end-numPeriods + 1:结束),MCFY (end-numPeriods + 1:结束,1),“。r”，.．.日期(end-numPeriods + 1:结束),CIFY (end-numPeriods + 1:最终,1,1),“- b”，.．.日期(end-numPeriods + 1:结束),CIFY (end-numPeriods + 1:结束,1,2),“- b”，.．.日期(end-numPeriods + 1:结束),Y (: 1),”:c '，.．.日期(end-numPeriods + 1:结束),磅(:1),“m”，.．.日期(end-numPeriods + 1:结束),乌兰巴托(:1),“m”，.．.“线宽”3);包含(“时间”) ylabel (“失业率的变化”)传说(h((1、2、6)){“观察”，“MC预测”，.．.“95%预测区间”，理论预测的，.．.“95%理论区间”}，“位置”，“最佳”)标题(“失业率的观察和预测变化”）

图中包含一个轴对象。标题为“失业率的观察和预测变化”的轴对象包含7个类型线对象。这些对象代表观测值、MC预测、95%预测区间、理论预测、95%理论区间。

figure h = plot(dates(end- numperiods -9:end)，[isY(end-9:end,2);oosY(:，2)]，，“- k”，.．.日期(end-numPeriods + 1:结束),MCFY (end-numPeriods + 1:结束,2),“。r”，.．.日期(end-numPeriods + 1:结束),CIFY (end-numPeriods + 1:结束,2,- 1),“- b”，.．.日期(end-numPeriods + 1:结束),CIFY (end-numPeriods + 1:最终,2,2),“- b”，.．.日期(end-numPeriods + 1:结束),Y (:, 2),”:c '，.．.日期(end-numPeriods + 1:结束),磅(:,2),“m”，.．.日期(end-numPeriods + 1:结束),乌兰巴托(:,2),“m”，.．.“线宽”3);包含(“时间”) ylabel (国民生产总值增长率)传说(h((1、2、6)){“观察”，“MC预测”，.．.'95% MC间隔'，理论预测的，“95%理论区间”}，.．.“位置”，“最佳”)标题(“观察和预测的国民生产总值增长率”）

图中包含一个轴对象。标题为“观察到的和预测的nGNP增长率”的轴对象包含7个类型为线的对象。这些对象代表观测值、MC预测、95% MC区间、理论预测、95%理论区间。

另请参阅

舰导弹|模拟|估计|完善|预测

用蒙特卡罗方法预测状态空间模型

另请参阅

相关的例子

更多关于