%% Refine Estimation of State-Space Model Containing Regression Component
% Suppose that the relationship between the unemployment rate and the
% nominal gross national product (nGNP) is linear. Suppose further that the
% unemployment rate is an AR(1) series. Symbolically, and in state-space
% form, the model is
%
% $$\begin{array}{l}
% x_{t} = \phi x_{t - 1} + \sigma u_{t}\\
% {y_t} - \beta {Z_t} = {x_{t}}, 
% \end{array}$$
%
% where:
%
% * $x_{t}$ is the unemployment rate at time _t_. 
% * $y_{t}$ is the observed unemployment rate being deflated by the log of
% nGNP ($Z_t$).
% * $u_{t}$ is the Gaussian series of state disturbances having mean 0 and
% unknown standard deviation $\sigma$.
%
%%
% Load the Nelson-Plosser data set, which contains the unemployment rate
% and nGNP series data.
load Data_NelsonPlosser
%%
% Preprocess the data by taking the first difference of the unemployment
% rate and converting nGNP to a series of returns. Remove the observations
% corresponding to the sequence of |NaN| values at the beginning of the
% unemployment rate series.
isNaN = any(ismissing(DataTable),2);  % Flag periods containing NaNs
gnpn = DataTable.GNPN(~isNaN);
y = DataTable.UR(~isNaN);
y = diff(y); 
T = size(y,1);
Z = [ones(T,1) price2ret(gnpn)];
%%
% This example continues using the series without |NaN| values.  However, using
% the Kalman filter framework, the software can accommodate series
% containing missing values.
%%
% Specify the coefficient matrices.
A = NaN;
B = NaN;
C = 1;
%%
% Specify the state-space model using |ssm|.
Mdl = ssm(A,B,C);
%%
% Find a good set of starting parameters to use for estimation.
params0 = [150 1000]; % Initial values chosen arbitrarily
Beta0 = [1 -100];
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
%%
% Estimate |Mdl| using the refined initial parameter values
% |refinedParams0|.
EstMdl = estimate(Mdl,y,refinedParams0(1:(end - 2)),'Predictors',Z,...
    'Beta0',refinedParams0((end - 1):end));
%%
% |estimate| returns reasonable parameter estimates and their corresponding
% standard errors.

%% 
% Copyright 2012 The MathWorks, Inc.