%% Refine Initial Parameters After Fitting State-Space Models
% This example shows how to choose a better set of initial parameter values to pass to estimate when you estimate ssm models.
% State-space models tend to induce likelihood surfaces that have multiple maxima. Also, the software uses numerical optimization to find the parameter estimates. Therefore, it is important to use a good initial value for the parameters to attain the global maximum of the likelihood function.
% Suppose that a latent process is a random walk. The state equation is
% x_t = x_{t-1} + u_t,
% where u_t is Gaussian with mean 0 and standard deviation 1.
% Generate a random series of 100 observations from x_t, assuming that the series starts at 1.5.

T = 100;
rng(1); % For reproducibility
u = randn(T,1);
x = cumsum([1.5;u]);
x = x(2:end);

% 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 1.
%
% Use the random latent state process (x) and the observation equation to generate observations.

y = x + randn(T,1);

% Together, the latent process and observation equations compose a state-space model. Assume that the state is a stationary AR(1) process. Then, 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 coefficient matrices. Use NaN values for unknown parameters.

A = NaN;
B = NaN;
C = 1;
D = NaN;

%
% Specify the state-space model using the coefficient matrices. Specify that the initial state distribution is stationary using the StateType name-value pair argument.

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

% Mdl is an ssm model. The software sets values for the initial state mean and variance. Verify that the model is correctly specified using the display in the Command Window.
%
% Pass the observations to estimate to estimate the parameters. Set starting values for the parameters to params0 that are likely not their corresponding true values. Also, specify lower bound constraints of 0 for the standard deviations.

params0 = [-1e7 1e-6 2000];
EstMdl = estimate(Mdl,y,params0,'lb',[-Inf,0,0]);

% estimate failed to converge, and so the results are undesirable.
%
% Refine params0 using refine.

Output = refine(Mdl,y,params0);
logL = cell2mat({Output.LogLikelihood})';
[~,maxLogLIndx] = max(logL)
refinedParams0 = Output(maxLogLIndx).Parameters
Description = Output(maxLogLIndx).Description

% The algorithm that yields the highest loglikelihood value is Loose bound interior point, which is the third struct in the structure array Output.
%
% Estimate Mdl using refinedParams0, which is the vector of refined initial parameter values.

EstMdl = estimate(Mdl,y,refinedParams0,'lb',[-Inf,0,0]);

% estimate converged, making the parameter estimates much more desirable.
% Copyright 2015 The MathWorks, Inc.