%% Compare Simulation Smoother to Smoothed States
% This example shows how the results of the state-space model simulation smoother (simsmooth) compare to the smoothed states (smooth).
% 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.
%
% 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,y,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.
%
% Simulate 1e4 paths of observations from the fitted, state-space model EstMdl using the simulation smoother. Specify to simulate observations for each period.

numPaths = 1e4;
SimX = simsmooth(EstMdl,y,'NumPaths',numPaths);

% SimX is a T - 1-by- 4-by- numPaths matrix containing the simulated states. The rows of SimX correspond to periods, the columns correspond to a state in the model, and the pages correspond to paths.
%
% Estimate the smoothed state means, standard deviations, and 95% confidence intervals.

SmoothBar = mean(SimX,3);
SmoothSTD = std(SimX,0,3);
SmoothCIL = SmoothBar - 1.96*SmoothSTD;
SmoothCIU = SmoothBar + 1.96*SmoothSTD;

%
% Estimate smooth states using smooth.

SmoothX = smooth(EstMdl,y);

%
% Plot the smoothed states, and the means of the simulated states and their 95% confidence intervals.

figure
h = plot(dates(2:T),SmoothBar(:,1),'-r',...
    dates(2:T),SmoothCIL(:,1),':b',...
    dates(2:T),SmoothCIU(:,1),':b',...
    dates(2:T),SmoothX(:,1),':k',...
    'LineWidth',3);
xlabel 'Period';
ylabel 'Unemployment rate';
legend(h([1,2,4]),{'Simulated, smoothed state mean','95% confidence interval',...
    'Smoothed states'},'Location','Best');
title 'Smoothed Unemployment Rate';
axis tight

figure
h = plot(dates(2:T),SmoothBar(:,3),'-r',...
    dates(2:T),SmoothCIL(:,3),':b',...
    dates(2:T),SmoothCIU(:,3),':b',...
    dates(2:T),SmoothX(:,3),':k',...
    'LineWidth',3);
xlabel 'Period';
ylabel 'nGNP';
legend(h([1,2,4]),{'Simulated, smoothed state mean','95% confidence interval',...
    'Smoothed states'},'Location','Best');
title 'Smoothed nGNP';
axis tight

% The simulated state means are practically identical to the smoothed states.
% Copyright 2015 The MathWorks, Inc.