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

    %% Forecast Observations of 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 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, among other things.

% Copyright 2015 The MathWorks, Inc.

load Data_NelsonPlosser
%%
% Preprocess the data by taking the natural logarithm of the nGNP series,
% and the first difference of each series.   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);
%%
% Though this example removes missing values, the software can accommodate
% series containing missing values in the Kalman filter framework.
%%
% To determine how well the model forecasts observations, remove the last
% 10 observations for comparison.
numPeriods = 10;                   % Forecast horizon
isY = y(1:end-numPeriods);         % In-sample observations
oosY = y(end-numPeriods+1:end);    % Out-of-sample observations
ISZ = Z(1:end-numPeriods,:);       % In-sample predictors
OOSZ = Z(end-numPeriods+1:end,:);  % Out-of-sample predictors
%%
% Specify the coefficient matrices.
A = [NaN NaN; 0 0];
B = [1; 1];
C = [1 0];
D = NaN;
%%
% Specify the state-space model using |ssm|.
Mdl = ssm(A,B,C,D);
%%
% Estimate the model parameters, and use a random set of initial parameter
% values for optimization. Specify the regression component and its initial
% value for optimization using the |'Predictors'| and |'Beta0'| name-value
% pair arguments, respectively. Restrict the estimate of $\sigma$ to all
% positive, real numbers.  For numerical stability, specify the Hessian
% when the software computes the parameter covariance matrix, using the
% |'CovMethod'| name-value pair argument.
params0 = [0.3 0.2 0.1]; % Chosen arbitrarily
[EstMdl,estParams] = estimate(Mdl,isY,params0,'Predictors',ISZ,...
    'Beta0',[0.1 0.2],'lb',[-Inf,-Inf,0,-Inf,-Inf],'CovMethod','hessian');
%%
% |EstMdl| is an |ssm| model, and you can access its properties using dot
% notation.
%%
% Forecast observations over the forecast horizon. |EstMdl| does not store
% the data set, so you must pass it in appropriate name-value pair
% arguments. 
[fY,yMSE] = forecast(EstMdl,numPeriods,isY,'Predictors0',ISZ,...
    'PredictorsF',OOSZ,'Beta',estParams(end-1:end));
%%
% |fY| is a 10-by-1 vector containing the forecasted observations, and
% |yMSE| is a 10-by-1 vector containing the variances of the forecasted
% observations.
%%
% Obtain 95% Wald-type forecast intervals. Plot the forecasted observations
% with their true values and the forecast intervals.
ForecastIntervals(:,1) = fY - 1.96*sqrt(yMSE);
ForecastIntervals(:,2) = fY + 1.96*sqrt(yMSE);

figure
h = plot(dates(end-numPeriods-9:end-numPeriods),isY(end-9:end),'-k',...
    dates(end-numPeriods+1:end),oosY,'-k',...
    dates(end-numPeriods+1:end),fY,'--r',...
    dates(end-numPeriods+1:end),ForecastIntervals,':b',...
    dates(end-numPeriods:end-numPeriods+1),...
    [isY(end)*ones(3,1),[oosY(1);ForecastIntervals(1,:)']],':k',...
    'LineWidth',2);
xlabel('Period')
ylabel('Change in the unemployment rate')
legend(h([1,3,4]),{'Observations','Forecasted responses',...
    '95% forecast intervals'})
title('Observed and Forecasted Changes in the Unemployment Rate')
%%
% This model seems to forecast the changes in the unemployment rate well.