Simulate Responses of Estimated VARX Model

This example shows how to estimate a multivariate time series model that contains lagged endogenous and exogenous variables, and how to simulate responses. The response series are the quarterly:

• Changes in real gross domestic product (rGDP) rates (${y}_{1t}$)

• Real money supply (rM1SL) rates (${y}_{2t}$)

• Short-term interest rates (i.e., three-month treasury bill yield, ${y}_{3t}$)

from March, 1959 through March, 2009 . The exogenous series is the quarterly changes in the unemployment rates (${x}_{t}$).

Suppose that a model for the responses is this VARX(4,3) model

$\left[\begin{array}{c}{y}_{1t}\\ {y}_{2t}\\ {y}_{3t}\end{array}\right]=\left[\begin{array}{c}{a}_{1}\\ {a}_{2}\\ {a}_{3}\end{array}\right]+\left[\begin{array}{ccc}{x}_{t}& 0& 0\\ 0& {x}_{t}& 0\\ 0& 0& {x}_{t}\end{array}\right]\left[\begin{array}{c}{b}_{1}\\ {b}_{2}\\ {b}_{3}\end{array}\right]+\left[\begin{array}{ccc}{\varphi }_{11}^{\left(1\right)}& {\varphi }_{12}^{\left(1\right)}& {\varphi }_{13}^{\left(1\right)}\\ {\varphi }_{21}^{\left(1\right)}& {\varphi }_{22}^{\left(1\right)}& {\varphi }_{23}^{\left(1\right)}\\ {\varphi }_{31}^{\left(1\right)}& {\varphi }_{32}^{\left(1\right)}& {\varphi }_{33}^{\left(1\right)}\end{array}\right]\left[\begin{array}{c}{y}_{1,t-1}\\ {y}_{2,t-1}\\ {y}_{3,t-1}\end{array}\right]+...+\left[\begin{array}{ccc}{\varphi }_{11}^{\left(4\right)}& {\varphi }_{12}^{\left(4\right)}& {\varphi }_{13}^{\left(4\right)}\\ {\varphi }_{21}^{\left(4\right)}& {\varphi }_{22}^{\left(4\right)}& {\varphi }_{23}^{\left(4\right)}\\ {\varphi }_{31}^{\left(4\right)}& {\varphi }_{32}^{\left(4\right)}& {\varphi }_{33}^{\left(4\right)}\end{array}\right]\left[\begin{array}{c}{y}_{1,t-1}\\ {y}_{2,t-1}\\ {y}_{3,t-1}\end{array}\right]+\left[\begin{array}{c}{\epsilon }_{1t}\\ {\epsilon }_{2t}\\ {\epsilon }_{3t}\end{array}\right]$

Preprocess the Data

Load the U.S. macroeconomic data set. Flag the series and their periods that contain missing values (indicated by NaN values).

varNaN = any(ismissing(DataTable),1); % Variables containing NaN values
seriesWithNaNs = series(varNaN)
seriesWithNaNs = 1x5 cell
Columns 1 through 3

{'(FEDFUNDS) Eff...'}    {'(GS10) Ten-yea...'}    {'(M1SL) M1 mone...'}

Columns 4 through 5

{'(M2SL) M2 mone...'}    {'(UNRATE) Unemp...'}

In this data set, the variables that contain missing values entered the sample later than the other variables. There are no missing values after sampling started for a particular variable.

Flag all periods corresponding to missing values in the model variables.

idx = all(~ismissing(DataTable(:,{'UNRATE' 'M1SL'})),2);

For the rest of the example, consider only those values that of the series indicated by a true in idx.

Compute rGDP and rM1SL, and the growth rates of rGDP, rM1SL, short-term interest rates, and the unemployment rate. Description contains a description of the data and the variable names. Reserve the last three years of data to investigate the out-of-sample performance of the estimated model.

rGDP = DataTable.GDP(idx)./(DataTable.GDPDEF(idx)/100);
rM1SL = DataTable.M1SL(idx)./(DataTable.GDPDEF(idx)/100);

dLRGDP = diff(log(rGDP));              % rGDP growth rate
dLRM1SL = diff(log(rM1SL));            % rM1SL growth rate
d3MTB = diff(DataTable.TB3MS(idx));    % Change in short-term interest rate (3MTB)
dUNRATE = diff(DataTable.UNRATE(idx)); % Change in unemployment rate

T = numel(d3MTB);    % Total sample size
oosT = 12;           % Out-of-sample size
estT = T - oosT;     % Estimation sample size
estIdx = 1:estT;     % Estimation sample indices
oosIdx = (T - 11):T; % Out-of-sample indices
dates = dates((end - T + 1):end);

EstY = [dLRGDP(estIdx) dLRM1SL(estIdx) d3MTB(estIdx)]; % In-sample responses
estX = dUNRATE(estIdx);                                % In-sample exogenous data
n = size(EstY,2);

OOSY = [dLRGDP(oosIdx) dLRM1SL(oosIdx) d3MTB(oosIdx)]; % Out-of-sample responses
oosX = dUNRATE(oosIdx);                                % Out-of-sample exogenous data

Create VARX Model

Create a VARX(4) model using varm.

Mdl = varm(n,4);

Mdl is a varm model object serving as a template for estimation. Currently, Mdl does know have the structure in place for the regression component. However, MATLAB® creates the required structure during estimation.

Estimate the VAR(4) Model

Estimate the parameters of the VARX(4) model using estimate. Display the parameter estimates.

EstMdl = estimate(Mdl,EstY,'X',estX);
summarize(EstMdl)

AR-Stationary 3-Dimensional VARX(4) Model with 1 Predictor

Effective Sample Size: 184
Number of Estimated Parameters: 42
LogLikelihood: 1037.52
AIC: -1991.04
BIC: -1856.01

Value       StandardError    TStatistic      PValue
___________    _____________    __________    __________

Constant(1)      0.0080266     0.00097087         8.2674     1.3688e-16
Constant(2)     0.00063838      0.0015942        0.40044        0.68883
Constant(3)       0.068361          0.143        0.47803        0.63263
AR{1}(1,1)       -0.034045        0.06633       -0.51327        0.60776
AR{1}(2,1)      -0.0024555        0.10891      -0.022546        0.98201
AR{1}(3,1)         -1.7163           9.77       -0.17567        0.86056
AR{1}(1,2)       -0.013882       0.046481       -0.29867        0.76519
AR{1}(2,2)         0.17753       0.076323          2.326       0.020017
AR{1}(3,2)         -6.7572         6.8464       -0.98697        0.32366
AR{1}(1,3)       0.0010682     0.00048092         2.2212       0.026337
AR{1}(2,3)      -0.0050252     0.00078967        -6.3636     1.9705e-10
AR{1}(3,3)        -0.16256       0.070837        -2.2948       0.021744
AR{2}(1,1)        0.077748       0.064014         1.2145        0.22454
AR{2}(2,1)       0.0047257        0.10511       0.044959        0.96414
AR{2}(3,1)          3.4244         9.4289        0.36318        0.71647
AR{2}(1,2)        0.077867       0.046954         1.6584       0.097245
AR{2}(2,2)         0.29087       0.077099         3.7727     0.00016148
AR{2}(3,2)         0.39284         6.9161         0.0568         0.9547
AR{2}(1,3)      -0.0010719     0.00056413        -1.9001       0.057423
AR{2}(2,3)      -0.0016135     0.00092631        -1.7419       0.081533
AR{2}(3,3)        -0.21556       0.083094        -2.5942      0.0094802
AR{3}(1,1)       -0.090881       0.062563        -1.4526        0.14633
AR{3}(2,1)        0.064249        0.10273        0.62542        0.53169
AR{3}(3,1)         -7.9727         9.2152       -0.86517        0.38695
AR{3}(1,2)       -0.024092        0.04631       -0.52024        0.60289
AR{3}(2,2)        0.068565       0.076041        0.90168        0.36723
AR{3}(3,2)          10.263         6.8212         1.5046        0.13242
AR{3}(1,3)     -0.00055981     0.00056073       -0.99836        0.31811
AR{3}(2,3)      -0.0021302     0.00092073        -2.3136        0.02069
AR{3}(3,3)         0.22969       0.082593         2.7809      0.0054203
AR{4}(1,1)        0.066151       0.056841         1.1638        0.24451
AR{4}(2,1)        0.028826       0.093334        0.30885        0.75744
AR{4}(3,1)          1.0379         8.3724        0.12397        0.90134
AR{4}(1,2)       -0.078735       0.043804        -1.7975       0.072263
AR{4}(2,2)       0.0096425       0.071926        0.13406        0.89335
AR{4}(3,2)         -12.007         6.4521        -1.8609       0.062761
AR{4}(1,3)     -0.00018454     0.00053356       -0.34586        0.72945
AR{4}(2,3)     -0.00019036     0.00087611       -0.21728        0.82799
AR{4}(3,3)        0.053812       0.078591        0.68471        0.49353
Beta(1,1)        -0.016084      0.0016037        -10.029     1.1365e-23
Beta(2,1)         -0.00154      0.0026333       -0.58482        0.55867
Beta(3,1)          -1.5317        0.23622        -6.4841     8.9252e-11

Innovations Covariance Matrix:
0.0000    0.0000    0.0000
0.0000    0.0001   -0.0019
0.0000   -0.0019    0.7790

Innovations Correlation Matrix:
1.0000    0.1198    0.0011
0.1198    1.0000   -0.2177
0.0011   -0.2177    1.0000

EstMdl is a varm model object containing the estimated parameters.

Simulate Out-Of-Sample Response Paths

Simulate 1000, 3 year response series paths from the estimated model assuming that the exogenous unemployment rate is a fixed series. Since the model contains 4 lags per endogenous variable, specify the last 4 observations in the estimation sample as presample data.

numPaths = 1000;
Y0 = EstY((end-3):end,:);
rng(1); % For reproducibility
YSim = simulate(EstMdl,oosT,'X',oosX,'Y0',Y0,'NumPaths',numPaths);

YSim is a 12-by-3-by-1000 numeric array of simulated responses. The rows of YSim correspond to out-of-sample periods, the columns correspond to the response series, and the pages correspond to paths.

Plot the response data and the simulated responses. Identify the 5%, 25%, 75% and 95% percentiles, and the mean and median of the simulated series at each out-of-sample period.

YSimBar = mean(YSim,3);
YSimQrtl = quantile(YSim,[0.05 0.25 0.5 0.75 0.95],3);
RepDates = repmat(dates(oosIdx),1,1000);
respNames = {'dLRGDP' 'dLRM1SL' 'd3MTB'};

figure;
for j = 1:n
subplot(3,1,j);
h1 = plot(dates(oosIdx),squeeze(YSim(:,j,:)),'Color',0.75*ones(3,1));
hold on;
h2 = plot(dates(oosIdx),YSimBar(:,j),'.-k','LineWidth',2);
h3 = plot(dates(oosIdx),squeeze(YSimQrtl(:,j,:)),':r','LineWidth',1.5);
h4 = plot(dates((end - 30):end),[EstY((end - 18):end,j);OOSY(:,j)],...
'b','LineWidth',2);
title(sprintf('%s',respNames{j}));
datetick;
axis tight;
hold off;
end
legend([h1(1) h2(1) h3(1) h4],{'Simulated Series','Simulation Mean',...
'Simulation Quartile','Data'},'Location',[0.4 0.1 0.01 0.01],...
'FontSize',8); Simulate Out-Of-Sample Response Paths Using Random Exogenous Data

Suppose that the change in the unemployment rate is an AR(4) model, and fit the model to the estimation sample data.

MdlUNRATE = arima('ARLags',1:4);
EstMdlUNRATE = estimate(MdlUNRATE,estX,'Display','off');

EstMdlUNRATE is an arima model object containing the parameter estimates.

Simulate 1000, 3 year paths from the estimated AR(4) model for the change in unemployment rate. Since the model contains 4 lags, specify the last 4 observations in the estimation sample as presample data.

XSim = simulate(EstMdlUNRATE,oosT,'Y0',estX(end-3:end),...
'NumPaths',numPaths);

XSim is a 12-by-1000 numeric matrix of simulated exogenous paths. The rows correspond to periods and the columns correspond to paths.

Simulate 3 years of 1000 future response series paths from the estimated model using the simulated exogenous data. simulate does not accept multiple paths of predictor data, so you must simulate responses in a loop. Since the model contains 4 lags per endogenous variable, specify the last 4 observations in the estimation sample as presample data.

YSimRX = zeros(oosT,n,numPaths); % Preallocate

for j = 1:numPaths
YSimRX(:,:,j) = simulate(EstMdl,oosT,'X',XSim(:,j),'Y0',Y0);
end

YSimRX is a 12-by-3-by-1000 numeric array of simulated responses.

Plot the response data and the simulated responses. Identify the 5%, 25%, 75% and 95% percentiles, and the mean and median of the simulated series at each out-of-sample period.

YSimBarRX = mean(YSimRX,3);
YSimQrtlRX = quantile(YSimRX,[0.05 0.25 0.5 0.75 0.95],3);

figure;
for j = 1:n;
subplot(3,1,j);
h1 = plot(dates(oosIdx),squeeze(YSimRX(:,j,:)),'Color',0.75*ones(3,1));
hold on;
h2 = plot(dates(oosIdx),YSimBarRX(:,j),'.-k','LineWidth',2);
h3 = plot(dates(oosIdx),squeeze(YSimQrtlRX(:,j,:)),':r','LineWidth',1.5);
h4 = plot(dates((end - 30):end),[EstY((end - 18):end,j);OOSY(:,j)],...
'b','LineWidth',2);
title(sprintf('%s with Simulated Unemployment Rate',respNames{j}));
datetick;
axis tight;
hold off;
end
legend([h1(1) h2(1) h3(1) h4],{'Simulated Series','Simulation Mean',...
'Simulation Quartile','Data'},'Location',[0.4 0.1 0.01 0.01],...
'FontSize',8) 