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

    %% Simulate and Forecast a VEC Model
% This example shows how to generate forecasts from a VEC(_q_) model several
% ways:
%
% * Monte Carlo forecast paths using the VEC(_q_) model directly
% * Minimum mean square error (MMSE) forecasts using the VAR model
% representation of the VEC(_q_) model.
% * Monte Carlo forecast paths using the VAR(_p_) model representation
% of the VEC(_q_) model.
% 
% This example follows from <docid:econ_ug.bsxb_m3-1>.
%% Load Data and Preprocess
% Load the |Data_Canada| data set. Extract the interest rate data.
load Data_Canada
Y = Data(:,3:end);
%% Estimate VEC(_q_) Model
% Assuming that the interest rate data follows a VEC(2) model, fit the
% model to the data.
[~,~,~,~,reg] = egcitest(Y,'test','t2');
c0 = reg.coeff(1);
b = reg.coeff(2:3);
beta = [1; -b];
q = 2;
[numObs,numDims] = size(Y);
tBase = (q+2):numObs;                % Commensurate time base, all lags
T = length(tBase);                   % Effective sample size
DeltaYLags = zeros(T,(q+1)*numDims);
YLags = lagmatrix(Y,0:(q+1));        % Y(t-k) on observed time base
LY = YLags(tBase,(numDims+1):2*numDims);
for k = 1:(q+1)
    DeltaYLags(:,((k-1)*numDims+1):k*numDims) = ...
               YLags(tBase,((k-1)*numDims+1):k*numDims) ...
             - YLags(tBase,(k*numDims+1):(k+1)*numDims);
end
DY = DeltaYLags(:,1:numDims);        % (1-L)Y(t)
DLY = DeltaYLags(:,(numDims+1):end); % [(1-L)Y(t-1),...,(1-L)Y(t-q)]
X = [(LY*beta-c0),DLY,ones(T,1)];
P = (X\DY)';                         % [alpha,B1,...,Bq,c1]
alpha = P(:,1);
C = alpha*beta';                     % Error-correction coefficient matrix
B1 = P(:,2:4);                       % VEC(2) model coefficient
B2 = P(:,5:7);                       % VEC(2) model coefficient
c1 = P(:,end);
b = (alpha*c0 + c1)';                % VEC(2) model constant offset
res = DY-X*P';
EstCov = cov(res);
%% Monte Carlo Forecasts Using VEC Model
% Specify a 10-period forecast horizon. Set |numPaths| to generate 1000
% paths. Because _q_, the degree of the VEC model, is 2, reserve three
% presample observations for $y_t$.
forecastHorizon = 10;
numPaths        = 1000;
psSize          = 3;

PSY = repmat(Y(end-(psSize-1):end,:),1,1,numPaths); % Presample
YSimVEC = zeros(forecastHorizon,numDims,numPaths);  % Preallocate forecasts
YSimVEC = [PSY; YSimVEC];
%%
% Generate Monte Carlo forecasts by adding the simulated innovations to the
% estimated VEC(2) model.
rng('default'); % For reproducibility
for p = 1:numPaths
    for t = (1:forecastHorizon) + psSize;
        eps = mvnrnd([0 0 0],EstCov);
        YSimVEC(t,:,p) = YSimVEC(t-1,:,p) + (C*YSimVEC(t-1,:,p)')'...
                    + (YSimVEC(t-1,:,p) - YSimVEC(t - 2,:,p))*B1'...
                    + (YSimVEC(t-2,:,p) - YSimVEC(t - 3,:,p))*B2'...
                    + b + eps;
    end
end
%%
% |YSimVEC| is a 13-by-3-by-1000 numeric array.  Its rows correspond to
% presample and forecast periods, columns correspond to the time series,
% and the pages correspond to a draw.
%%
% Compute the mean of the forecasts for each period
% and time series over all paths.  Construct 95% percentile forecast
% intervals for each period and time series.
FMCVEC   = mean(YSimVEC((psSize + 1):end,:,:),3);
CIMCVEC = quantile(YSimVEC((psSize + 1):end,:,:),[0.25,0.975],3);
%%
% |FMCVEC| is a 10-by-3 numeric matrix containing the Monte Carlo forecasts
% for each period (row) and time series (column). |CIMCVEC| is a
% 10-by-3-by-2 numeric array containing the 2.5% and 97.5% percentiles
% (pages) of the draws for each period (row) and time series (column).
%%
% Plot the effective-sample observations, the mean forecasts, and the 95%
% percentile confidence intervals.
fDates = dates(end) + (0:forecastHorizon)';
figure;
h1 = plot([dates; fDates(2:end)],[Y; FMCVEC],'LineWidth',2);
h2 = gca;
hold on
h3 = plot(repmat(fDates,1,3),[Y(end,:,:); CIMCVEC(:,:,1)],'--',...
    'LineWidth',2);
h3(1).Color = h1(1).Color;
h3(2).Color = h1(2).Color;
h3(3).Color = h1(3).Color;
h4 = plot(repmat(fDates,1,3),[Y(end,:,:); CIMCVEC(:,:,2)],'--',...
    'LineWidth',2);
h4(1).Color = h1(1).Color;
h4(2).Color = h1(2).Color; 
h4(3).Color = h1(3).Color;
patch([fDates(1) fDates(1) fDates(end) fDates(end)],...
    [h2.YLim(1) h2.YLim(2) h2.YLim(2) h2.YLim(1)],'b','FaceAlpha',0.1)
xlabel('Year')
ylabel('Percent')
title('{\bf VEC Model Monte Carlo Forecasts}')
axis tight
grid on
legend(h1,DataTable.Properties.VariableNames(3:end),'Location','Best');
%%
% The plot suggests that |INT_S| has lower forecast accuracy that the other
% two series because its confidence bounds are the widest.
%% MMSE Forecasts Using VAR(p) Representation  
% Compute the autoregressive coefficient matrices of the equivalent VAR(3)
% model to the VEC(2) model (i.e., the coefficients of $y_{t-1}$,
% $y_{t-2}$, and $y_{t-3}$).
A = vec2var({B1 B2},C);
%%
% |A| is a 1-by-3 row cell vector.  |A{j}| is the autoregressive
% coefficient matrix for lag term |j|. The constant offset (|b|) of the
% VEC(2) model and the constant offset of the equivalent VAR(3) model are
% equal.
%%
% Create a VAR(3) model object.
VAR3 = vgxset('AR',A,'a',b,'Q',EstCov);
%%
% Forecast over a 10 period horizon.  Compute 95% individual, Wald-type
% confidence intervals for each series.
[MMSEF,CovF] = vgxpred(VAR3,forecastHorizon,[],Y);
var = cellfun(@diag,CovF,'UniformOutput',false);
CIF = zeros(forecastHorizon,numDims,2); % Preallocation
for j = 1:forecastHorizon
    stdev = sqrt(var{j});
    CIF(j,:,1) = MMSEF(j,:) - 1.96*stdev';
    CIF(j,:,2) = MMSEF(j,:) + 1.96*stdev';
end
%%
% The confidence intervals do not account for the correlation between
% the forecasted series at a particular time.
%%
% |MMSEF| is a 10-by-3 numeric matrix of the MMSE forecasts for the VAR(3)
% model.  Rows correspond to forecast periods and columns to time series.
% |CovF| is a 10-by-1 cell vector of forecast covariance matrices, in which
% each row corresponds to a forecast period.
%%
% Plot the effective-sample observations, the MMSE forecasts, and the 95%
% Wald confidence intervals.
figure;
h1 = plot([dates; fDates(2:end)],[Y; MMSEF],'LineWidth',2);
h2 = gca;
hold on
h3 = plot(repmat(fDates,1,3),[Y(end,:,:); CIF(:,:,1)],'--',...
    'LineWidth',2);
h3(1).Color = h1(1).Color;
h3(2).Color = h1(2).Color;
h3(3).Color = h1(3).Color;
h4 = plot(repmat(fDates,1,3),[Y(end,:,:); CIF(:,:,2)],'--',...
    'LineWidth',2);
h4(1).Color = h1(1).Color;
h4(2).Color = h1(2).Color; 
h4(3).Color = h1(3).Color;
patch([fDates(1) fDates(1) fDates(end) fDates(end)],...
    [h2.YLim(1) h2.YLim(2) h2.YLim(2) h2.YLim(1)],'b','FaceAlpha',0.1)
xlabel('Year')
ylabel('Percent')
title('{\bf VAR Model MMSE Forecasts}')
axis tight
grid on
legend(h1,DataTable.Properties.VariableNames(3:end),'Location','Best');
%%
% The MMSE forecasts are very close to the VEC(2) model Monte Carlo forecast
% means.  However, the confidence bounds are wider.  
%% Monte Carlo Forecasts Using VAR(p) Representation
% Simulate 1000 paths of the VAR(3) model into the forecast horizon.  Use
% the observations as presample data.
YSimVAR = vgxsim(VAR3,forecastHorizon,[],Y,[],numPaths);
%%
% |YSimVAR| is a 10-by-3-by-1000 numeric array similar to |YSimVEC|.
%%
% Compute the mean of the forecasts for each period
% and time series over all paths.  Construct 95% percentile forecast
% intervals for each period and time series.
FMCVAR   = mean(YSimVAR,3);
CIMCVAR  = quantile(YSimVAR,[0.25,0.975],3);
%%
% |FMCVAR| is a 10-by-3 numeric matrix containing the Monte Carlo forecasts
% for each period and time series. |CIMCVAR| is a 10-by-3-by-2
% numeric array containing the 2.5% and 97.5% percentiles of the
% draws for each period and time series.
%%
% Plot the effective-sample observations, the mean forecasts, and the 95%
% percentile confidence intervals.
figure;
h1 = plot([dates; fDates(2:end)],[Y; FMCVAR],'LineWidth',2);
h2 = gca;
hold on
h3 = plot(repmat(fDates,1,3),[Y(end,:,:); CIMCVAR(:,:,1)],'--',...
    'LineWidth',2);
h3(1).Color = h1(1).Color;
h3(2).Color = h1(2).Color;
h3(3).Color = h1(3).Color;
h4 = plot(repmat(fDates,1,3),[Y(end,:,:); CIMCVAR(:,:,2)],'--',...
    'LineWidth',2);
h4(1).Color = h1(1).Color;
h4(2).Color = h1(2).Color; 
h4(3).Color = h1(3).Color;
patch([fDates(1) fDates(1) fDates(end) fDates(end)],...
    [h2.YLim(1) h2.YLim(2) h2.YLim(2) h2.YLim(1)],'b','FaceAlpha',0.1)
xlabel('Year')
ylabel('Percent')
title('{\bf VAR Representation Monte Carlo Forecasts}')
axis tight
grid on
legend(h1,DataTable.Properties.VariableNames(3:end),'Location','Best');
%%
% All sets of mean forecasts are very close.  Both the VAR(3) and VEC(2)
% sets of Monte Carlo forecasts are almost equivalent.