www.gusucode.com > econ 案例源码程序 matlab代码 > econ/EstimateMonteCarloForecastsOfAStateSpaceModelExample.m

    %% Estimate Monte-Carlo Forecasts of State-Space Model
% Suppose that the relationship between the change in the unemployment rate
% ($x_{1,t}$) and the nominal gross national product (nGNP) growth rate
% ($x_{3,t}$) can be expressed in the following, state-space model form.
%
% $$ {\left[ {\begin{array}{*{20}{c}}
% {{x_{1,t}}}\\
% {{x_{2,t}}}\\
% {{x_{3,t}}}\\
% {{x_{4,t}}}
% \end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
% {{\phi _1}}&{{\theta _1}}&{{\gamma _1}}&0\\
% 0&0&0&0\\
% {{\gamma _2}}&0&{{\phi _2}}&{{\theta _2}}\\
% 0&0&0&0
% \end{array}} \right]\left[ {\begin{array}{*{20}{c}}
% {{x_{1,t - 1}}}\\
% {{x_{2,t - 1}}}\\
% {{x_{3,t - 1}}}\\
% {{x_{4,t - 1}}}
% \end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
% 1&0\\
% 1&0\\
% 0&1\\
% 0&1
% \end{array}} \right]\left[ {\begin{array}{*{20}{c}}
% {{u_{1,t}}}\\
% {{u_{2,t}}}
% \end{array}} \right]}$$
%
% $${\left[ {\begin{array}{*{20}{c}}y_{1,t} \\ y_{2,t}\end{array}}\right]} = \left[\begin{array}{*{20}{c}}
% 1 & 0  & 0 & 0 \\
% 0 & 0  & 1 & 0 \\
% \end{array}\right]
% \left[ {\begin{array}{*{20}{c}}x_{1,t} \\ x_{2,t}\\x_{3,t}\\x_{4,t}\end{array}}\right]
% + \left[ {\begin{array}{*{20}{c}}\sigma_1 & 0\\ 0 &
% \sigma_2\end{array}}\right]\left[\begin{array}{*{20}{c}}\varepsilon_{1,t} \\ \varepsilon_{2,t}\end{array}\right],$$
%
% where:
%
% * $x_{1,t}$ is the change in the unemployment rate at time _t_.
% * $x_{2,t}$ is a dummy state for the MA(1) effect on $x_{1,t}$.
% * $x_{3,t}$ is the nGNP growth rate at time _t_.
% * $x_{4,t}$ is a dummy state for the MA(1) effect on $x_{3,t}$.
% * $y_{1,t}$ is the observed change in the unemployment rate.
% * $y_{2,t}$ is the observed nGNP growth rate.
% * $u_{1,t}$ and $u_{2,t}$ are Gaussian series of state disturbances having mean 0 and
% standard deviation 1.
% * $\varepsilon_{1,t}$ is the Gaussian series of observation innovations having
% mean 0 and standard deviation $\sigma_1$.
% * $\varepsilon_{2,t}$ is the Gaussian series of observation innovations having
% mean 0 and standard deviation $\sigma_2$.
%
%%
% Load the Nelson-Plosser data set, which contains the unemployment rate and
% nGNP series, among other things.

% Copyright 2015 The MathWorks, Inc.

load Data_NelsonPlosser
%%
% Preprocess the data by taking the natural logarithm of the nGNP series,
% and the first difference of each.   Also, remove the starting |NaN|
% values from each series.
isNaN = any(ismissing(DataTable),2);       % Flag periods containing NaNs
gnpn = DataTable.GNPN(~isNaN);
u = DataTable.UR(~isNaN);
T = size(gnpn,1);                          % Sample size
y = zeros(T-1,2);                          % Preallocate
y(:,1) = diff(u);
y(:,2) = diff(log(gnpn));
%%
% This example proceeds using series without |NaN| values.  However, using
% the Kalman filter framework, the software can accommodate series
% containing missing values.
%%
% To determine how well the model forecasts observations, remove the last
% 10 observations for comparison.
numPeriods = 10;                   % Forecast horizon
isY = y(1:end-numPeriods,:);       % In-sample observations
oosY = y(end-numPeriods+1:end,:);  % Out-of-sample observations
%%
% Specify the coefficient matrices.
A = [NaN NaN NaN 0; 0 0 0 0; NaN 0 NaN NaN; 0 0 0 0];
B = [1 0;1 0 ; 0 1; 0 1];
C = [1 0 0 0; 0 0 1 0];
D = [NaN 0; 0 NaN];
%%
% Specify the state-space model using |ssm|. Verify that the model
% specification is consistent with the state-space model.
Mdl = ssm(A,B,C,D)
%%
% Estimate the model parameters, and use a random set of initial parameter
% values for optimization. Restrict the estimate of $\sigma_1$ and
% $\sigma_2$ to all positive, real numbers using the |'lb'| name-value pair
% argument.  For numerical stability, specify the Hessian when the software
% computes the parameter covariance matrix, using the |'CovMethod'|
% name-value pair argument.
rng(1);
params0 = rand(8,1);
[EstMdl,estParams] = estimate(Mdl,isY,params0,...
    'lb',[-Inf -Inf -Inf -Inf -Inf -Inf 0 0],'CovMethod','hessian');
%%
% |EstMdl| is an |ssm| model, and you can access its properties using dot
% notation.
%%
% Filter the estimated, state-space model, and extract the filtered states
% and their variances from the final period.
[~,~,Output] = filter(EstMdl,isY);
%%
% Modify the estimated, state-space model so that the initial state means
% and covariances are the filtered states and their covariances of the
% final period.  This sets up simulation over the forecast horizon.
EstMdl1 = EstMdl;
EstMdl1.Mean0 = Output(end).FilteredStates;
EstMdl1.Cov0 = Output(end).FilteredStatesCov;
%%
% Simulate |5e5| paths of observations from the fitted, state-space model
% |EstMdl|.  Specify to simulate observations for each period.
numPaths = 5e5;
SimY = simulate(EstMdl1,10,'NumPaths',numPaths);
%%
% |SimY| is a |10|-by- |2|-by- |numPaths| array containing the simulated
% observations. The rows of |SimY| correspond to periods, the columns
% correspond to an observation in the model, and the pages correspond to
% paths.
%%
% Estimate the forecasted observations and their 95% confidence intervals
% in the forecast horizon.
MCFY = mean(SimY,3);
CIFY = quantile(SimY,[0.025 0.975],3);
%%
% Estimate the theoretical forecast bands.
[Y,YMSE] = forecast(EstMdl,10,isY);
Lb = Y - sqrt(YMSE)*1.96;
Ub = Y + sqrt(YMSE)*1.96;
%%
% Plot the forecasted observations with their true values and the forecast
% intervals.
figure
h = plot(dates(end-numPeriods-9:end),[isY(end-9:end,1);oosY(:,1)],'-k',...
    dates(end-numPeriods+1:end),MCFY(end-numPeriods+1:end,1),'.-r',...
    dates(end-numPeriods+1:end),CIFY(end-numPeriods+1:end,1,1),'-b',...
    dates(end-numPeriods+1:end),CIFY(end-numPeriods+1:end,1,2),'-b',...
    dates(end-numPeriods+1:end),Y(:,1),':c',...
    dates(end-numPeriods+1:end),Lb(:,1),':m',...
    dates(end-numPeriods+1:end),Ub(:,1),':m',...
    'LineWidth',3);
xlabel('Period')
ylabel('Change in the unemployment rate')
legend(h([1,2,4:6]),{'Observations','MC forecasts',...
    '95% forecast intervals','Theoretical forecasts',...
    '95% theoretical intervals'},'Location','Best')
title('Observed and Forecasted Changes in the Unemployment Rate')

figure
h = plot(dates(end-numPeriods-9:end),[isY(end-9:end,2);oosY(:,2)],'-k',...
    dates(end-numPeriods+1:end),MCFY(end-numPeriods+1:end,2),'.-r',...
    dates(end-numPeriods+1:end),CIFY(end-numPeriods+1:end,2,1),'-b',...
    dates(end-numPeriods+1:end),CIFY(end-numPeriods+1:end,2,2),'-b',...
    dates(end-numPeriods+1:end),Y(:,2),':c',...
    dates(end-numPeriods+1:end),Lb(:,2),':m',...
    dates(end-numPeriods+1:end),Ub(:,2),':m',...
    'LineWidth',3);
xlabel('Period')
ylabel('nGNP growth rate')
legend(h([1,2,4:6]),{'Observations','MC forecasts',...
    '95% MC intervals','Theoretical forecasts','95% theoretical intervals'},...
    'Location','Best')
title('Observed and Forecasted nGNP Growth Rates')