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

    %% Forecast a VAR Model Using Monte Carlo Simulation
% This example shows how to use Monte Carlo simulation via |vgxsim| to forecast a VAR
% model.  
%%
% |vgxsim| enables you to generate simulations of time series based on your
% model. If you have a trustworthy model structure, you can use these
% simulations as sample forecasts.
%%
% vgxsim requires:
% 
% * A model (|impSpec| in what follows)
% * The number of periods for the forecast (|FT| in what follows)
%%
% vgxsim optionally takes:
% 
% * An exogenous data series
% * A presample time series (|Y(end-3:end,:)| in what follows)
% * Presample innovations
% * The number of realizations to simulate (|2000| in what follows)
%%
% Load the |Data_USEconModel| data set.  This example uses two time series: the logarithm of
% real GDP, and the real 3-month T-bill rate, both differenced to be
% approximately stationary. For illustration, a VAR(4) model describes the
% time series.

% Copyright 2015 The MathWorks, Inc.

load Data_USEconModel
DEF = log(DataTable.CPIAUCSL);
GDP = log(DataTable.GDP);
rGDP = diff(GDP - DEF); % Real GDP is GDP - deflation
TB3 = 0.01*DataTable.TB3MS;
dDEF = 4*diff(DEF); % Scaling
rTB3 = TB3(2:end) - dDEF; % Real interest is deflated
Y = [rGDP,rTB3];
%%
% Fit a VAR(4) model specification:
Spec = vgxset('n',2,'nAR',4,'Constant',true);
impSpec = vgxvarx(Spec,Y(5:end,:),[],Y(1:4,:));
impSpec = vgxset(impSpec,'Series',...
  {'Transformed real GDP','Transformed real 3-mo T-bill rate'});
%%
% Define the forecast horizon.
FDates = datenum({'30-Jun-2009'; '30-Sep-2009'; '31-Dec-2009';
'31-Mar-2010'; '30-Jun-2010'; '30-Sep-2010'; '31-Dec-2010';
'31-Mar-2011'; '30-Jun-2011'; '30-Sep-2011'; '31-Dec-2011';
'31-Mar-2012'; '30-Jun-2012'; '30-Sep-2012'; '31-Dec-2012';
'31-Mar-2013'; '30-Jun-2013'; '30-Sep-2013'; '31-Dec-2013';
'31-Mar-2014'; '30-Jun-2014' });
FT = numel(FDates);
%%
% Simulate the model for 10 steps, replicated 2000 times:
rng(1); %For reproducibility
Ysim = vgxsim(impSpec,FT,[],Y(end-3:end,:),[],2000);
%%
% Calculate the mean and standard deviation of the simulated series:
Ymean = mean(Ysim,3); % Calculate means
Ystd = std(Ysim,0,3); % Calculate std deviations
%%
% Plot the means +/- 1 standard deviation for the simulated series:
subplot(2,1,1)
plot(dates(end-10:end),Y(end-10:end,1),'k')
hold('on')
plot([dates(end);FDates],[Y(end,1);Ymean(:,1)],'r')
plot([dates(end);FDates],[Y(end,1);Ymean(:,1)]+[0;Ystd(:,1)],'b')
plot([dates(end);FDates],[Y(end,1);Ymean(:,1)]-[0;Ystd(:,1)],'b')
datetick('x')
title('Transformed real GDP')
subplot(2,1,2)
plot(dates(end-10:end),Y(end-10:end,2),'k')
hold('on')
axis([dates(end-10),FDates(end),-.1,.1]);
plot([dates(end);FDates],[Y(end,2);Ymean(:,2)],'r')
plot([dates(end);FDates],[Y(end,2);Ymean(:,2)]+[0;Ystd(:,2)],'b')
plot([dates(end);FDates],[Y(end,2);Ymean(:,2)]-[0;Ystd(:,2)],'b')
datetick('x')
title('Transformed real 3-mo T-bill rate')