%% Extract Other Estimates from Output
% Estimate a diffuse state-space model, filter the states, and then extract other estimates from the Output output argument.
% Consider the diffuse state-space model
% \begin{array}{l}
% \left[ {\begin{array}{*{20}{c}}
% {{x_{1,t}}}\\
% {{x_{2,t}}}
% \end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
% \phi &0\\
% 0&1
% \end{array}} \right]\left[ {\begin{array}{*{20}{c}}
% {{x_{1,t - 1}}}\\
% {{x_{2,t - 1}}}
% \end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
% {{\sigma _1}}&0\\
% 0&{{\sigma _2}}
% \end{array}} \right]\left[ {\begin{array}{*{20}{c}}
% {{u_{1,t}}}\\
% {{u_{2,t}}}
% \end{array}} \right]\\
% y_t = \left[ {\begin{array}{*{20}{c}}
% 1&1
% \end{array}} \right]\left[ {\begin{array}{*{20}{c}}
% {{x_{1,t}}}\\
% {{x_{2,t}}}
% \end{array}} \right]
% \end{array}
% The state variable x_{1,t} is an AR(1) model with autoregressive coefficient \phi.  x_{2,t} is a random walk. The disturbances u_{1,t} and u_{2,t} are independent Gaussian random variables with mean 0 and standard deviations \sigma_1 and \sigma_2, respectively. The observation y_{t} is the error-free sum of x_{1,t} and x_{2,t}.
% Generate data from the state-space model. To simulate the data, suppose that the sample size T = 100, \phi = 0.6, \sigma_1 = 0.2, \sigma_2 = 0.1, and x_{1,0} = x_{2,0} = 2.

rng(1); % For reproducibility
T = 100;
ARMdl = arima('AR',0.6,'Constant',0,'Variance',0.2^2);
x1 = simulate(ARMdl,T,'Y0',2);
u3 = 0.1*randn(T,1);
x3 = cumsum([2;u3]);
x3 = x3(2:end);
y = x1 + x3;

%
% Specify the coefficient matrices of the state-space model. To indicate unknown parameters, use NaN values.

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

%
% Create a diffuse state-space model that describes the model above. Specify that x_{1,t} and x_{2,t} have diffuse initial state distributions.

StateType = [2 2];
Mdl = dssm(A,B,C,'StateType',StateType);

%
% Estimate the unknown parameters of Mdl. Choose initial parameter values for optimization. Specify that the standard deviations are constrained to be positive, but all other parameters are unconstrained using the 'lb' name-value pair argument.

params0 = [0.01 0.1 0.01]; % Initial values chosen arbitrarily
EstMdl = estimate(Mdl,y,params0,'lb',[-Inf 0 0]);

% The parameters are close to their true values.
%
% Filter the states of EstMdl, and request all other available output.

[X,logL,Output] = filter(EstMdl,y);

%
% X is a T-by-2 matrix of filtered states, logL is the final optimized log-likelihood value, and Output is a structure array containing various estimates that the Kalman filter requires. List the fields of output using fields.

fields(Output)

%
% Convert Output to a table.

OutputTbl = struct2table(Output);
OutputTbl(1:10,1:5) % Display first ten rows of first five variables

%
% The first two rows of the table contain empty cells or zeros. These correspond to the observations required to initialize the diffuse Kalman filter. That is, SwitchTime is 2.

SwitchTime = 2;

%
% Plot the filtered and forecasted states.

ForeX = cell2mat(OutputTbl.ForecastedStates')'; % Orient forecasted states
ForeX = [zeros(SwitchTime,2);ForeX]; % Include zeros for initialization

figure;
plot(1:T,X(:,1),'r',1:T,ForeX(:,1),'b');
xlabel('Period');
ylabel('State estimate');
title('State 1 Estimates')
legend('Filtered','Forecasted');
grid on;

figure;
plot(1:T,X(:,2),'r',1:T,ForeX(:,2),'b');
xlabel('Period');
ylabel('State estimate');
title('State 2 Estimates')
legend('Filtered','Forecasted');
grid on;

% Copyright 2012 The MathWorks, Inc.