www.gusucode.com > econ 案例源码程序 matlab代码 > econ/EstimateMonteCarloForecastsOfAStateSpaceModelExample.m
%% Estimate Monte-Carlo Forecasts of State-Space Model % Suppose that the relationship between the change in the unemployment rate % ($x_{1,t}$) and the nominal gross national product (nGNP) growth rate % ($x_{3,t}$) can be expressed in the following, state-space model form. % % $$ {\left[ {\begin{array}{*{20}{c}} % {{x_{1,t}}}\\ % {{x_{2,t}}}\\ % {{x_{3,t}}}\\ % {{x_{4,t}}} % \end{array}} \right] = \left[ {\begin{array}{*{20}{c}} % {{\phi _1}}&{{\theta _1}}&{{\gamma _1}}&0\\ % 0&0&0&0\\ % {{\gamma _2}}&0&{{\phi _2}}&{{\theta _2}}\\ % 0&0&0&0 % \end{array}} \right]\left[ {\begin{array}{*{20}{c}} % {{x_{1,t - 1}}}\\ % {{x_{2,t - 1}}}\\ % {{x_{3,t - 1}}}\\ % {{x_{4,t - 1}}} % \end{array}} \right] + \left[ {\begin{array}{*{20}{c}} % 1&0\\ % 1&0\\ % 0&1\\ % 0&1 % \end{array}} \right]\left[ {\begin{array}{*{20}{c}} % {{u_{1,t}}}\\ % {{u_{2,t}}} % \end{array}} \right]}$$ % % $${\left[ {\begin{array}{*{20}{c}}y_{1,t} \\ y_{2,t}\end{array}}\right]} = \left[\begin{array}{*{20}{c}} % 1 & 0 & 0 & 0 \\ % 0 & 0 & 1 & 0 \\ % \end{array}\right] % \left[ {\begin{array}{*{20}{c}}x_{1,t} \\ x_{2,t}\\x_{3,t}\\x_{4,t}\end{array}}\right] % + \left[ {\begin{array}{*{20}{c}}\sigma_1 & 0\\ 0 & % \sigma_2\end{array}}\right]\left[\begin{array}{*{20}{c}}\varepsilon_{1,t} \\ \varepsilon_{2,t}\end{array}\right],$$ % % 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 on $x_{1,t}$. % * $x_{3,t}$ is the nGNP growth rate at time _t_. % * $x_{4,t}$ is a dummy state for the MA(1) effect on $x_{3,t}$. % * $y_{1,t}$ is the observed change in the unemployment rate. % * $y_{2,t}$ is the observed nGNP growth rate. % * $u_{1,t}$ and $u_{2,t}$ are Gaussian series of state disturbances having mean 0 and % standard deviation 1. % * $\varepsilon_{1,t}$ is the Gaussian series of observation innovations having % mean 0 and standard deviation $\sigma_1$. % * $\varepsilon_{2,t}$ is the Gaussian series of observation innovations having % mean 0 and standard deviation $\sigma_2$. % %% % 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. 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 y = zeros(T-1,2); % Preallocate y(:,1) = diff(u); y(:,2) = diff(log(gnpn)); %% % This example proceeds using series without |NaN| values. However, using % the Kalman filter framework, the software can accommodate series % containing missing values. %% % 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 %% % Specify the coefficient matrices. A = [NaN NaN NaN 0; 0 0 0 0; NaN 0 NaN NaN; 0 0 0 0]; B = [1 0;1 0 ; 0 1; 0 1]; C = [1 0 0 0; 0 0 1 0]; D = [NaN 0; 0 NaN]; %% % Specify the state-space model using |ssm|. Verify that the model % specification is consistent with the state-space model. Mdl = ssm(A,B,C,D) %% % Estimate the model parameters, and use a random set of initial parameter % values for optimization. Restrict the estimate of $\sigma_1$ and % $\sigma_2$ to all positive, real numbers using the |'lb'| name-value pair % argument. For numerical stability, specify the Hessian when the software % computes the parameter covariance matrix, using the |'CovMethod'| % name-value pair argument. rng(1); params0 = rand(8,1); [EstMdl,estParams] = estimate(Mdl,isY,params0,... 'lb',[-Inf -Inf -Inf -Inf -Inf -Inf 0 0],'CovMethod','hessian'); %% % |EstMdl| is an |ssm| model, and you can access its properties using dot % notation. %% % Filter the estimated, state-space model, and extract the filtered states % and their variances from the final period. [~,~,Output] = filter(EstMdl,isY); %% % Modify the estimated, state-space model so that the initial state means % and covariances are the filtered states and their covariances of the % final period. This sets up simulation over the forecast horizon. EstMdl1 = EstMdl; EstMdl1.Mean0 = Output(end).FilteredStates; EstMdl1.Cov0 = Output(end).FilteredStatesCov; %% % Simulate |5e5| paths of observations from the fitted, state-space model % |EstMdl|. Specify to simulate observations for each period. numPaths = 5e5; SimY = simulate(EstMdl1,10,'NumPaths',numPaths); %% % |SimY| is a |10|-by- |2|-by- |numPaths| array containing the simulated % observations. The rows of |SimY| correspond to periods, the columns % correspond to an observation in the model, and the pages correspond to % paths. %% % Estimate the forecasted observations and their 95% confidence intervals % in the forecast horizon. MCFY = mean(SimY,3); CIFY = quantile(SimY,[0.025 0.975],3); %% % Estimate the theoretical forecast bands. [Y,YMSE] = forecast(EstMdl,10,isY); Lb = Y - sqrt(YMSE)*1.96; Ub = Y + sqrt(YMSE)*1.96; %% % Plot the forecasted observations with their true values and the forecast % intervals. figure h = plot(dates(end-numPeriods-9:end),[isY(end-9:end,1);oosY(:,1)],'-k',... dates(end-numPeriods+1:end),MCFY(end-numPeriods+1:end,1),'.-r',... dates(end-numPeriods+1:end),CIFY(end-numPeriods+1:end,1,1),'-b',... dates(end-numPeriods+1:end),CIFY(end-numPeriods+1:end,1,2),'-b',... dates(end-numPeriods+1:end),Y(:,1),':c',... dates(end-numPeriods+1:end),Lb(:,1),':m',... dates(end-numPeriods+1:end),Ub(:,1),':m',... 'LineWidth',3); xlabel('Period') ylabel('Change in the unemployment rate') legend(h([1,2,4:6]),{'Observations','MC forecasts',... '95% forecast intervals','Theoretical forecasts',... '95% theoretical intervals'},'Location','Best') title('Observed and Forecasted Changes in the Unemployment Rate') figure h = plot(dates(end-numPeriods-9:end),[isY(end-9:end,2);oosY(:,2)],'-k',... dates(end-numPeriods+1:end),MCFY(end-numPeriods+1:end,2),'.-r',... dates(end-numPeriods+1:end),CIFY(end-numPeriods+1:end,2,1),'-b',... dates(end-numPeriods+1:end),CIFY(end-numPeriods+1:end,2,2),'-b',... dates(end-numPeriods+1:end),Y(:,2),':c',... dates(end-numPeriods+1:end),Lb(:,2),':m',... dates(end-numPeriods+1:end),Ub(:,2),':m',... 'LineWidth',3); xlabel('Period') ylabel('nGNP growth rate') legend(h([1,2,4:6]),{'Observations','MC forecasts',... '95% MC intervals','Theoretical forecasts','95% theoretical intervals'},... 'Location','Best') title('Observed and Forecasted nGNP Growth Rates')