%% 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.

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')

% Copyright 2015 The MathWorks, Inc.