%% Estimate State-Space Model Containing Regression Component
% Suppose that the linear 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\\
% 1
% \end{array}} \right]u_{1,t}\\
% {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_{1,t} is the observed change in the unemployment 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, among other things.

load Data_NelsonPlosser

%
% 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 proceeds using series without NaN values. However, using the Kalman filter framework, the software can accommodate series containing missing values.
%
% Specify the state-transition coefficient matrix.

A = [NaN NaN; 0 0];

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

B = [1; 1];

%
% Specify the measurement-sensitivity coefficient matrix.

C = [1 0];

%
% Specify the observation-innovation coefficient matrix.

D = NaN;

%
% Specify the state-space model using ssm.

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

%
% Estimate the model parameters. Specify the regression component and its initial value for optimization using the 'Predictors' and 'Beta0' name-value pair arguments, respectively. Display the estimates and all optimization diagnostic information. Restrict the estimate of \sigma to all positive, real numbers.

params0 = [0.3 0.2 0.1]; % Chosen arbitrarily
EstMdl = estimate(Mdl,y,params0,'Predictors',Z,'Display','full',...
    'Beta0',[0.1 0.2],'lb',[-Inf,-Inf,0,-Inf,-Inf]);

% Optimization information and a table of estimates and statistics output to the Command Window.  EstMdl is an ssm model, and you can access its properties using dot notation.
% Copyright 2012 The MathWorks, Inc.