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.