%% Refine Estimation of State-Space Model Containing Regression Component
% This example shows how to choose a better set of initial parameter values to pass to estimate when you estimate ssm models containing a regression component.
% 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 the relationship between the change in the unemployment rate and the nominal gross national product (nGNP) growth rate is of interest. Suppose further that the first difference of the unemployment rate is an ARMA(1,1) series. Symbolically, and in state-space form, the model is
% \begin{array}{l}
% \left[ {\begin{array}{*{20}{c}}
% {{x_{1,t}}}\\
% {{x_{2,t}}}
% \end{array}} \right] = \left[ {\begin{array}{*{20}{c}}
% \phi &\theta \\
% 0&0
% \end{array}} \right]\left[ {\begin{array}{*{20}{c}}
% {{x_{1,t - 1}}}\\
% {{x_{2,t - 1}}}
% \end{array}} \right] + \left[ {\begin{array}{*{20}{c}}
% 1&0\\
% 1&0
% \end{array}} \right]\left[ {\begin{array}{*{20}{c}}
% {{u_{1,t}}}\\
% {{u_{2,t}}}
% \end{array}} \right]\\
% {y_t} - \beta {Z_t} = {x_{1,t}} + \sigma\varepsilon_t,
% \end{array}
% 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.
% y_{t} is the observed change in the unemployment rate being deflated by the growth rate of nGNP (Z_t).
% u_{1,t} is the Gaussian series of state disturbances having mean 0 and standard deviation 1.
% \varepsilon_t is the Gaussian series of observation innovations having mean 0 and standard deviation \sigma.
% Load the Nelson-Plosser data set, which contains the unemployment rate and nGNP series data.

load Data_NelsonPlosser
rng(1); % For reproducibility

%
% 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
Z = [ones(T-1,1) diff(log(gnpn))];
y = diff(u);

% This example continues excluding NaN values from the analysis. However, using the Kalman filter framework, the software can accommodate series containing missing values.
%
% Specify the coefficient matrices.

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

%
% Specify the state-space model using ssm.

Mdl = ssm(A,B,C,D);

%
% To help you choose a good set of initial parameter values for estimation, pass Mdl and an initial set of parameter values to refine.

params0 = [0 0 0]; % Chosen arbitrarily
Beta0 = [0 0];
Output = refine(Mdl,y,params0,'Predictors',Z,'Beta0',Beta0);

% Output is a 1-by-5 structure array containing the recommended initial parameter values.
%
% Choose the initial parameter values corresponding to the largest loglikelihood.

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

% The algorithm that yields the highest loglikelihood value is Starting value perturbation, which is the fourth struct in the structure array Output.
%
% Estimate Mdl using the refined initial parameter values refinedParams0.

pBeta = numel(Beta0);
EstMdl = estimate(Mdl,y,refinedParams0(1:(end - pBeta)),'Predictors',Z,...
    'Beta0',refinedParams0((end - pBeta + 1):end),...
    'lb',[-Inf,-Inf,0,-Inf,-Inf]);

% estimate returns reasonable parameter estimates and their corresponding standard errors.
% Copyright 2012 The MathWorks, Inc.