%% Compare Estimates from State-Space Model Filtering Methods
% The software implements the Kalman filter using the covariance filter by default, but you can specify to use the square-root filter instead. This example compares estimates from each method using simulated data.
% Suppose that a latent process is an AR(1). The state equation is
% x_t = 0.5x_{t-1} + u_t,
% where u_t is Gaussian with mean 0 and standard deviation 0.3.
% Generate a random series of 100 observations from x_t, assuming that the series starts at 1.5.

T = 100;
ARMdl = arima('AR',0.5,'Constant',0,'Variance',0.3^2);
x0 = 1.5;
rng(1); % For reproducibility
x = simulate(ARMdl,T,'Y0',x0);

% Suppose further that the latent process is subject to additive measurement error. The observation equation is
% y_t = x_t + \varepsilon_t,
% where \varepsilon_t is Gaussian with mean 0 and standard deviation 0.1.
%
% Use the random latent state process (x) and the observation equation to generate observations.

y = x + 0.1*randn(T,1);

% Together, the latent process and observation equations compose a state-space model. Supposing that the coefficients and variances are unknown parameter, the state-space model is
% \begin{array}{c}x_t = \phi x_{t-1} + \sigma_1u_t \\ y_t = x_t +
% \sigma_2\varepsilon_t\end{array}
%
% Specify the state-transition coefficient matrix. Use NaN values for unknown parameters.

A = NaN;

%
% Specify the state-disturbance-loading coefficient matrix.

B = NaN;

%
% Specify the measurement-sensitivity coefficient matrix.

C = 1;

%
% Specify the observation-innovation coefficient matrix.

D = NaN;

%
% Specify the state-space model using the coefficient matrices. Also, specify the initial state mean, variance, and distribution (which is stationary).

Mean0 = 0;
Cov0 = 10;
StateType = 0;
Mdl = ssm(A,B,C,D,'Mean0',Mean0,'Cov0',Cov0,'StateType',StateType);

% Mdl is an ssm model.
%
% Estimate the parameters using estimate two ways:
% Using the default, simple Kalman filter
% Using the square root filter variation
% In both cases, specify that no output should be returned to the Command Window. This is good practice if you plan on running estimate multiple times (such as a Monte Carlo simulation).

params0 = [10,10,10];
[~,estParamsSKF,EstParamCovSKF,logLSKF,OutputSKF] = estimate(Mdl,y,params0,...
    'Display','off');
[~,estParamsSR,EstParamCovSR,logLSR,OutputSR] = estimate(Mdl,y,params0,...
    'Squareroot',true,'Display','off');

%
% Check that the algorithms converged properly by printing the exit flag properties of OutputSKF and OutputSR.

exitFlagSKF = OutputSKF.ExitFlag
exitFlagSR = OutputSR.ExitFlag

% Both algorithms have an exit flag of 1, which indicates that the software met the convergence criteria.
%
% Compare the estimates from each algorithm.

fprintf('\n Parameter Estimates\n')
table(estParamsSKF',estParamsSR','VariableNames',...
    {'SimpleKalmanFilter','SquarerootFilter'})
fprintf('\nEstimated Parameter Covariance Matrix\n')
table(EstParamCovSKF,EstParamCovSR,'VariableNames',...
    {'SimpleKalmanFilter','SquarerootFilter'})

% In this case, the results are the same.
% If you use the default, covariance filter method, and you run into numerical problems during estimation, filtering, or smoothing, try using the squareroot method.
% Copyright 2015 The MathWorks, Inc.