www.gusucode.com > econ 案例源码程序 matlab代码 > econ/RefineEstimationOfAStateSpaceModelWithARegressionComponeExample.m

    %% 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.