simulate

Simulate coefficients and innovations covariance matrix of Bayesian vector autoregression (VAR) model

Syntax

[Coeff,Sigma] = simulate(PriorMdl)

[Coeff,Sigma] = simulate(PriorMdl,Y)

[Coeff,Sigma] = simulate(___,Name,Value)

Description

[Coeff,Sigma] = simulate(PriorMdl) returns a random vector of coefficients Coeff and a random innovations covariance matrix Sigma drawn from the prior Bayesian VAR(p) model PriorMdl.

[Coeff,Sigma] = simulate(PriorMdl,Y) draws from the posterior distributions produced or updated by incorporating the response data Y.

NaNs in the data indicate missing values, which simulate removes by using list-wise deletion.

[Coeff,Sigma] = simulate(___,Name,Value) specifies options using one or more name-value pair arguments in addition to any of the input argument combinations in the previous syntaxes. For example, you can set the number of random draws from the distribution or specify the presample response data.

Examples

Draw Coefficients and Innovations Covariance Matrix from Prior Distribution

Consider the 3-D VAR(4) model for the US inflation (INFL), unemployment (UNRATE), and federal funds (FEDFUNDS) rates.

$$\left[\begin{array}{l}{\text{INFL}}_t\\ {\text{UNRATE}}_t\\ {\text{FEDFUNDS}}_t\end{array}\right]=c+\sum _{j=1}^4{\Phi }_j\left[\begin{array}{l}{\text{INFL}}_{t-j}\\ {\text{UNRATE}}_{t-j}\\ {\text{FEDFUNDS}}_{t-j}\end{array}\right]+\left[\begin{array}{c}{\epsilon }_{1,t}\\ {\epsilon }_{2,t}\\ {\epsilon }_{3,t}\end{array}\right].$$

For all $$t$$, $${\epsilon }_t$$ is a series of independent 3-D normal innovations with a mean of 0 and covariance $\Sigma$. Assume that a conjugate prior distribution $\pi \left({\left[{\Phi }_1,...,{\Phi }_4,\mathit{c}\right]}^{\prime },\Sigma \right)$ governs the behavior of the parameters.

Create a conjugate prior model. Specify the response series names. Obtain a summary of the prior distribution.

seriesnames = ["INFL" "UNRATE" "FEDFUNDS"];
numseries = numel(seriesnames);
numlags = 4;

PriorMdl = bayesvarm(numseries,numlags,'ModelType','conjugate',...
    'SeriesNames',seriesnames);
Summary = summarize(PriorMdl,'off');

Draw a set of coefficients and an innovations covariance matrix from the prior distribution.

rng(1) % For reproducibility
[Coeff,Sigma] = simulate(PriorMdl);

Display the selected coefficients with corresponding names and the innovations covariance matrix.

table(Coeff,'RowNames',Summary.CoeffMap)

ans=39×1 table
                     Coeff   
                   __________

    AR{1}(1,1)        0.44999
    AR{1}(1,2)       0.047463
    AR{1}(1,3)      -0.042106
    AR{2}(1,1)     -0.0086264
    AR{2}(1,2)       0.034049
    AR{2}(1,3)      -0.058092
    AR{3}(1,1)      -0.015698
    AR{3}(1,2)      -0.053203
    AR{3}(1,3)      -0.031138
    AR{4}(1,1)       0.036431
    AR{4}(1,2)      -0.058279
    AR{4}(1,3)       -0.02195
    Constant(1)        -1.001
    AR{1}(2,1)      -0.068182
    AR{1}(2,2)        0.51029
    AR{1}(2,3)      -0.094367
      ⋮

AR{r}(j,k) is the AR coefficient of response variable k (lagged r units) in response equation j.

Sigma

Sigma = 3×3

    0.1238   -0.0053   -0.0369
   -0.0053    0.0456    0.0160
   -0.0369    0.0160    0.1039

Rows and columns of Sigma correspond to the innovations in the response equations ordered by PriorMdl.SeriesNames.

Simulate Parameters from Analytically Tractable Posterior Distribution

Consider the 3-D VAR(4) model of Draw Coefficients and Innovations Covariance Matrix from Prior Distribution. In this case, assume that the prior distribution is diffuse.

Load and Preprocess Data

Load the US macroeconomic data set. Compute the inflation rate, stabilize the unemployment and federal funds rates, and remove missing values.

load Data_USEconModel
seriesnames = ["INFL" "UNRATE" "FEDFUNDS"];
DataTimeTable.INFL = 100*[NaN; price2ret(DataTimeTable.CPIAUCSL)];

DataTimeTable.DUNRATE = [NaN; diff(DataTimeTable.UNRATE)];
DataTimeTable.DFEDFUNDS = [NaN; diff(DataTimeTable.FEDFUNDS)];
seriesnames(2:3) = "D" + seriesnames(2:3);
rmDataTimeTable = rmmissing(DataTimeTable);

Create Prior Model

Create a diffuse Bayesian VAR(4) prior model for the three response series. Specify the response series names.

numseries = numel(seriesnames);
numlags = 4;

PriorMdl = bayesvarm(numseries,numlags,'SeriesNames',seriesnames);

Estimate Posterior Distribution

Estimate the posterior distribution. Return the estimation summary.

[PosteriorMdl,Summary] = estimate(PriorMdl,rmDataTimeTable{:,seriesnames});

Bayesian VAR under diffuse priors
Effective Sample Size:          197
Number of equations:            3
Number of estimated Parameters: 39
             |   Mean     Std  
-------------------------------
 Constant(1) |  0.1007  0.0832 
 Constant(2) | -0.0499  0.0450 
 Constant(3) | -0.4221  0.1781 
 AR{1}(1,1)  |  0.1241  0.0762 
 AR{1}(2,1)  | -0.0219  0.0413 
 AR{1}(3,1)  | -0.1586  0.1632 
 AR{1}(1,2)  | -0.4809  0.1536 
 AR{1}(2,2)  |  0.4716  0.0831 
 AR{1}(3,2)  | -1.4368  0.3287 
 AR{1}(1,3)  |  0.1005  0.0390 
 AR{1}(2,3)  |  0.0391  0.0211 
 AR{1}(3,3)  | -0.2905  0.0835 
 AR{2}(1,1)  |  0.3236  0.0868 
 AR{2}(2,1)  |  0.0913  0.0469 
 AR{2}(3,1)  |  0.3403  0.1857 
 AR{2}(1,2)  | -0.0503  0.1647 
 AR{2}(2,2)  |  0.2414  0.0891 
 AR{2}(3,2)  | -0.2968  0.3526 
 AR{2}(1,3)  |  0.0450  0.0413 
 AR{2}(2,3)  |  0.0536  0.0223 
 AR{2}(3,3)  | -0.3117  0.0883 
 AR{3}(1,1)  |  0.4272  0.0860 
 AR{3}(2,1)  | -0.0389  0.0465 
 AR{3}(3,1)  |  0.2848  0.1841 
 AR{3}(1,2)  |  0.2738  0.1620 
 AR{3}(2,2)  |  0.0552  0.0876 
 AR{3}(3,2)  | -0.7401  0.3466 
 AR{3}(1,3)  |  0.0523  0.0428 
 AR{3}(2,3)  |  0.0008  0.0232 
 AR{3}(3,3)  |  0.0028  0.0917 
 AR{4}(1,1)  |  0.0167  0.0901 
 AR{4}(2,1)  |  0.0285  0.0488 
 AR{4}(3,1)  | -0.0690  0.1928 
 AR{4}(1,2)  | -0.1830  0.1520 
 AR{4}(2,2)  | -0.1795  0.0822 
 AR{4}(3,2)  |  0.1494  0.3253 
 AR{4}(1,3)  |  0.0067  0.0395 
 AR{4}(2,3)  |  0.0088  0.0214 
 AR{4}(3,3)  | -0.1372  0.0845 
       Innovations Covariance Matrix       
           |   INFL     DUNRATE  DFEDFUNDS 
-------------------------------------------
 INFL      |  0.3028   -0.0217     0.1579  
           | (0.0321)  (0.0124)   (0.0499) 
 DUNRATE   | -0.0217    0.0887    -0.1435  
           | (0.0124)  (0.0094)   (0.0283) 
 DFEDFUNDS |  0.1579   -0.1435     1.3872  
           | (0.0499)  (0.0283)   (0.1470) 

PosteriorMdl is a conjugatebvarm model, which is analytically tractable.

Simulate Parameters from Posterior

Draw 1000 samples from the posterior distribution.

rng(1)
[Coeff,Sigma] = simulate(PosteriorMdl,'NumDraws',1000);

Coeff is a 39-by-1000 matrix of randomly drawn coefficients. Each column is an individual draw, and each row is an individual coefficient. Sigma is a 3-by-3-by-1000 array of randomly drawn innovations covariance matrices. Each page is an individual draw.

Display the first coefficient drawn from the distribution with corresponding parameter names, and display the first drawn innovations covariance matrix.

Coeffs = table(Coeff(:,1),'RowNames',Summary.CoeffMap)

Coeffs=39×1 table
                     Var1   
                   _________

    AR{1}(1,1)       0.14994
    AR{1}(1,2)      -0.46927
    AR{1}(1,3)      0.088388
    AR{2}(1,1)       0.28139
    AR{2}(1,2)      -0.19597
    AR{2}(1,3)      0.049222
    AR{3}(1,1)        0.3946
    AR{3}(1,2)      0.081871
    AR{3}(1,3)      0.002117
    AR{4}(1,1)       0.13514
    AR{4}(1,2)      -0.23661
    AR{4}(1,3)      -0.01869
    Constant(1)     0.035787
    AR{1}(2,1)     0.0027895
    AR{1}(2,2)       0.62382
    AR{1}(2,3)      0.053232
      ⋮

Sigma(:,:,1)

ans = 3×3

    0.2653   -0.0075    0.1483
   -0.0075    0.1015   -0.1435
    0.1483   -0.1435    1.5042

Simulate Parameters from Analytically Intractable Posterior Distribution

Consider the 3-D VAR(4) model of Draw Coefficients and Innovations Covariance Matrix from Prior Distribution. In this case, assume that the prior distribution is semiconjugate.

Load and Preprocess Data

Load the US macroeconomic data set. Compute the inflation rate, stabilize the unemployment and federal funds rates, and remove missing values.

load Data_USEconModel
seriesnames = ["INFL" "UNRATE" "FEDFUNDS"];
DataTimeTable.INFL = 100*[NaN; price2ret(DataTimeTable.CPIAUCSL)];

DataTimeTable.DUNRATE = [NaN; diff(DataTimeTable.UNRATE)];
DataTimeTable.DFEDFUNDS = [NaN; diff(DataTimeTable.FEDFUNDS)];
seriesnames(2:3) = "D" + seriesnames(2:3);
rmDataTimeTable = rmmissing(DataTimeTable);

Create Prior Model

Create a semiconjugate Bayesian VAR(4) prior model for the three response series. Specify the response variable names.

numseries = numel(seriesnames);
numlags = 4;

PriorMdl = bayesvarm(numseries,numlags,'Model','semiconjugate',...
    'SeriesNames',seriesnames);

Simulate Parameters from Posterior

Because the joint posterior distribution of a semiconjugate prior model is analytically intractable, simulate sequentially draws from the full conditional distributions.

Draw 1000 samples from the posterior distribution. Specify a burn-in period of 10,000, and a thinning factor of 5. Start the Gibbs sampler by assuming the posterior mean of $\Sigma$ is the 3-D identity matrix.

rng(1)
[Coeff,Sigma] = simulate(PriorMdl,rmDataTimeTable{:,seriesnames},...
    'NumDraws',1000,'BurnIn',1e4,'Thin',5,'Sigma0',eye(3));

Coeff is a 39-by-1000 matrix of randomly drawn coefficients. Each column is an individual draw, and each row is an individual coefficient. Sigma is a 3-by-3-by-1000 array of randomly drawn innovations covariance matrices. Each page is an individual draw.

Simulate Coefficients from VARX Model

Consider the 2-D VARX(1) model for the US real GDP (RGDP) and investment (GCE) rates that treats the personal consumption (PCEC) rate as exogenous:

$$\left[\begin{array}{l}{\text{RGDP}}_t\\ {\text{GCE}}_t\end{array}\right]=c+\Phi \left[\begin{array}{l}{\text{RGDP}}_{t-1}\\ {\text{GCE}}_{t-1}\end{array}\right]+{\text{PCEC}}_t\beta +{\epsilon }_t.$$

For all $$t$$, $${\epsilon }_t$$ is a series of independent 2-D normal innovations with a mean of 0 and covariance $\Sigma$. Assume the following prior distributions:

• ${\left[\begin{array}{ccc}\Phi & \mathit{c}& \beta \end{array}\right]}^{\prime }|\Sigma \sim {{\rm N}}_{4\times 2}\left({\rm M},\mathit{V},\Sigma \right)$, where M is a 4-by-2 matrix of means and $\mathit{V}$ is the 4-by-4 among-coefficient scale matrix. Equivalently, $\mathrm{vec}\left({\left[\begin{array}{ccc}\Phi & \mathit{c}& \beta \end{array}\right]}^{\prime }\right)|\Sigma \sim {{\rm N}}_8\left(\mathrm{vec}\left({\rm M}\right),\Sigma \otimes \mathit{V}\right)$.

• $$\Sigma \sim InverseWishart(\Omega ,\nu ),$$ where Ω is the 2-by-2 scale matrix and $\nu$ is the degrees of freedom.

Load the US macroeconomic data set. Compute the real GDP, investment, and personal consumption rate series. Remove all missing values from the resulting series.

load Data_USEconModel
DataTimeTable.RGDP = DataTimeTable.GDP./DataTimeTable.GDPDEF;
seriesnames = ["PCEC"; "RGDP"; "GCE"];
rates = varfun(@price2ret,DataTimeTable,'InputVariables',seriesnames);
rates = rmmissing(rates);
rates.Properties.VariableNames = seriesnames;

Create a conjugate prior model for the 2-D VARX(1) model parameters.

numseries = 2;
numlags = 1;
numpredictors = 1;
PriorMdl = conjugatebvarm(numseries,numlags,'NumPredictors',numpredictors,...
    'SeriesNames',seriesnames(2:end));

Simulate directly from the posterior distribution. Specify the exogenous predictor data.

[Coeff,Sigma] = simulate(PriorMdl,rates{:,PriorMdl.SeriesNames},...
    'X',rates{:,seriesnames(1)});

By default, simulate uses the first p = 1 observations of the response data to initialize the dynamic component of the model, and removes the corresponding observations from the predictor data.

Input Arguments

PriorMdl — Prior Bayesian VAR model

conjugatebvarm model object | semiconjugatebvarm model object | diffusebvarm model object | normalbvarm model object

Prior Bayesian VAR model, specified as a model object in this table.

Model Object	Description
conjugatebvarm	Dependent, matrix-normal-inverse-Wishart conjugate model returned by bayesvarm or conjugatebvarm
semiconjugatebvarm	Independent, normal-inverse-Wishart semiconjugate prior model returned by bayesvarm or semiconjugatebvarm
diffusebvarm	Diffuse prior model returned by bayesvarm or diffusebvarm
normalbvarm	Normal conjugate model with a fixed innovations covariance matrix, returned by bayesvarm or normalbvarm

Note

If PriorMdl is a diffusebvarm model, then you must also supply Y because simulate cannot draw from an improper prior distribution. Consequently, Coeff and Sigma represent draws from the posterior distribution.

Y — Observed multivariate response series

numeric matrix

Observed multivariate response series to which simulate fits the model, specified as a numobs-by-numseries numeric matrix.

numobs is the sample size. numseries is the number of response variables (PriorMdl.NumSeries).

Rows correspond to observations, and the last row contains the latest observation. Columns correspond to individual response variables.

Y represents the continuation of the presample response series in Y0.

Data Types: double

Name-Value Arguments

Specify optional pairs of arguments as Name1=Value1,...,NameN=ValueN, where Name is the argument name and Value is the corresponding value. Name-value arguments must appear after other arguments, but the order of the pairs does not matter.

Before R2021a, use commas to separate each name and value, and enclose Name in quotes.

Example: 'Y0',Y0,'X',X specifies the presample response data Y0 to initialize the VAR model for posterior estimation, and the predictor data X for the exogenous regression component.

Options for All Prior Distributions

NumDraws — Number of random draws

1 (default) | positive integer

Number of random draws from the distributions, specified as the comma-separated pair consisting of 'NumDraws' and a positive integer.

Example: 'NumDraws',1e7

Data Types: double

Y0 — Presample response data

numeric matrix

Presample response data to initialize the VAR model for estimation, specified as the comma-separated pair consisting of 'Y0' and a numpreobs-by-numseries numeric matrix. numpreobs is the number of presample observations.

Rows correspond to presample observations, and the last row contains the latest observation. Y0 must have at least PriorMdl.P rows. If you supply more rows than necessary, simulate uses the latest PriorMdl.P observations only.

Columns must correspond to the response series in Y.

By default, simulate uses Y(1:PriorMdl.P,:) as presample observations, and then estimates the posterior using Y((PriorMdl.P + 1):end,:). This action reduces the effective sample size.

Data Types: double

X — Predictor data

numeric matrix

Predictor data for the exogenous regression component in the model, specified as the comma-separated pair consisting of 'X' and a numobs-by-PriorMdl.NumPredictors numeric matrix.

Rows correspond to observations, and the last row contains the latest observation. simulate does not use the regression component in the presample period. X must have at least as many observations as the observations used after the presample period.

• If you specify Y0, then X must have at least numobs rows (see Y).

• Otherwise, X must have at least numobs – PriorMdl.P observations to account for the presample removal.

In either case, if you supply more rows than necessary, simulate uses the latest observations only.

Columns correspond to individual predictor variables. All predictor variables are present in the regression component of each response equation.

Data Types: double

Options for Semiconjugate Prior Distributions

BurnIn — Number of draws to remove from beginning of sample

0 (default) | nonnegative scalar

Number of draws to remove from the beginning of the sample to reduce transient effects, specified as the comma-separated pair consisting of 'BurnIn' and a nonnegative scalar. For details on how simulate reduces the full sample, see Algorithms.

Tip

To help you specify the appropriate burn-in period size:

1. Determine the extent of the transient behavior in the sample by setting the BurnIn name-value argument to 0.

2. Simulate a few thousand observations by using simulate.

3. Create trace plots.

Example: 'BurnIn',0

Data Types: double

Thin — Adjusted sample size multiplier

1 (default) | positive integer

Adjusted sample size multiplier, specified as the comma-separated pair consisting of 'Thin' and a positive integer.

The actual sample size is BurnIn + NumDraws*Thin. After discarding the burn-in, simulate discards every Thin – 1 draws, and then retains the next draw. For more details on how simulate reduces the full sample, see Algorithms.

Tip

To reduce potential large serial correlation in the sample, or to reduce the memory consumption of the draws stored in Coeff and Sigma, specify a large value for Thin.

Example: 'Thin',5

Data Types: double

Coeff0 — Starting value of VAR model coefficients for Gibbs sampler

numeric column vector

Starting value of the VAR model coefficients for the Gibbs sampler, specified as the comma-separated pair consisting of 'Coeff0' and a numeric column vector with (PriorMdl.NumSeries*k)-by-NumDraws elements, where k = PriorMdl.NumSeries*PriorMdl.P + PriorMdl.IncludeIntercept + PriorMdl.IncludeTrend + PriorMdl.NumPredictors, which is the number of coefficients in a response equation. For details on the structure of Coeff0, see the output Coeff.

By default, Coeff0 is the multivariate least-squares estimate.

Tip

• To specify Coeff0:

  1. Set separate variables for the initial values each coefficient matrix and vector.

  2. Horizontally concatenate all coefficient means in this order:

     $$Coeff=\left[\begin{array}{ccccccc}{\Phi }_1& {\Phi }_2& \cdots & {\Phi }_p& c& \delta & {\rm B}\end{array}\right].$$

  3. Vectorize the transpose of the coefficient mean matrix.

     Coeff0 = Coeff.';
     Coeff0 = Coeff0(:);

• A good practice is to run simulate multiple times with different parameter starting values. Verify that the estimates from each run converge to similar values.

Data Types: double

Sigma0 — Starting value of innovations covariance matrix for Gibbs sampler

positive definite numeric matrix

Starting value of the innovations covariance matrix for the Gibbs sampler, specified as the comma-separated pair consisting of 'Sigma0' and a PriorMdl.NumSeries-by-PriorMdl.NumSeries positive definite numeric matrix. By default, Sigma0 is the residual mean squared error from multivariate least-squares. Rows and columns correspond to innovations in the equations of the response variables ordered by PriorMdl.SeriesNames.

Tip

A good practice is to run simulate multiple times with different parameter starting values. Verify that the estimates from each run converge to similar values.

Data Types: double

Output Arguments

Coeff — Simulated VAR model coefficients

numeric matrix

Simulated VAR model coefficients, returned as a (PriorMdl.NumSeries*k)-by-NumDraws numeric matrix, where k = PriorMdl.NumSeries*PriorMdl.P + PriorMdl.IncludeIntercept + PriorMdl.IncludeTrend + PriorMdl.NumPredictors, which is the number of coefficients in a response equation. Each column is a separate draw from the distribution.

For draw j, Coeff(1:k,j) corresponds to all coefficients in the equation of response variable PriorMdl.SeriesNames(1), Coeff((k + 1):(2*k),j) corresponds to all coefficients in the equation of response variable PriorMdl.SeriesNames(2), and so on. For a set of indices corresponding to an equation:

• Elements 1 through PriorMdl.NumSeries correspond to the lag 1 AR coefficients of the response variables ordered by PriorMdl.SeriesNames.

• Elements PriorMdl.NumSeries + 1 through 2*PriorMdl.NumSeries correspond to the lag 2 AR coefficients of the response variables ordered by PriorMdl.SeriesNames.

• In general, elements (q – 1)*PriorMdl.NumSeries + 1 through q*PriorMdl.NumSeries correspond to the lag q AR coefficients of the response variables ordered by PriorMdl.SeriesNames.

• If PriorMdl.IncludeConstant is true, element PriorMdl.NumSeries*PriorMdl.P + 1 is the model constant.

• If PriorMdl.IncludeTrend is true, element PriorMdl.NumSeries*PriorMdl.P + 2 is the linear time trend coefficient.

• If PriorMdl.NumPredictors > 0, elements PriorMdl.NumSeries*PriorMdl.P + 3 through k compose the vector of regression coefficients of the exogenous variables.

This figure shows the structure of Coeff(L,j) for a 2-D VAR(3) model that contains a constant vector and four exogenous predictors.

$$[\stackrel{y_{1,t}}{\overbrace{\begin{array}{ccccccccccc}{\varphi }_{1,11}& {\varphi }_{1,12}& {\varphi }_{2,11}& {\varphi }_{2,12}& {\varphi }_{3,11}& {\varphi }_{3,12}& c_1& {\beta }_{11}& {\beta }_{12}& {\beta }_{13}& {\beta }_{14}\end{array}}}\stackrel{y_{2,t}}{\overbrace{\begin{array}{ccccccccccc}{\varphi }_{1,21}& {\varphi }_{1,22}& {\varphi }_{2,21}& {\varphi }_{2,22}& {\varphi }_{3,21}& {\varphi }_{3,22}& c_2& {\beta }_{21}& {\beta }_{22}& {\beta }_{23}& {\beta }_{24}\end{array}}}],$$

where

• ϕ_q,jk is element (j,k) of the lag q AR coefficient matrix.

• c_j is the model constant in the equation of response variable j.

• B_ju is the regression coefficient of exogenous variable u in the equation of response variable j.

Sigma — Simulated innovations covariance matrices

array of positive definite numeric matrices

Simulated innovations covariance matrices, returned as a PriorMdl.NumSeries-by-PriorMdl.NumSeries-by-NumDraws array of positive definite numeric matrices.

Each page is a separate draw (covariance) from the distribution. Rows and columns correspond to innovations in the equations of the response variables ordered by PriorMdl.SeriesNames.

If PriorMdl is a normalbvarm object, all covariances in Sigma are equal to PriorMdl.Covariance.

Limitations

• simulate cannot draw values from an improper distribution, which is a distribution whose density does not integrate to 1.

More About

Bayesian Vector Autoregression (VAR) Model

A Bayesian VAR model treats all coefficients and the innovations covariance matrix as random variables in the m-dimensional, stationary VARX(p) model. The model has one of the three forms described in this table.

Model	Equation
Reduced-form VAR(p) in difference-equation notation	$$y_t={\Phi }_1{y}_{t-1}+\mathrm{...}+{\Phi }_p{y}_{t-p}+c+\delta t+{\rm B}{x}_t+{\epsilon }_t.$$
Multivariate regression	$$y_t=Z_t\lambda +{\epsilon }_t.$$
Matrix regression	$$y_t={\Lambda }^{\prime }{z}_t^{\prime }+{\epsilon }_t.$$

For each time t = 1,...,T:

• y_t is the m-dimensional observed response vector, where m = numseries.

• Φ₁,…,Φ_p are the m-by-m AR coefficient matrices of lags 1 through p, where p = numlags.

• c is the m-by-1 vector of model constants if IncludeConstant is true.

• δ is the m-by-1 vector of linear time trend coefficients if IncludeTrend is true.

• Β is the m-by-r matrix of regression coefficients of the r-by-1 vector of observed exogenous predictors x_t, where r = NumPredictors. All predictor variables appear in each equation.

• $$z_t=\left[\begin{array}{ccccccc}{y}_{t-1}^{\prime }& y_{t-2}^{\prime }& \cdots & y_{t-p}^{\prime }& 1& t& x_t^{\prime }\end{array}\right],$$ which is a 1-by-(mp + r + 2) vector, and Z_t is the m-by-m(mp + r + 2) block diagonal matrix

  $$\left[\begin{array}{cccc}{z}_t& 0_z& \cdots & 0_z\\ 0_z& z_t& \cdots & 0_z\\ ⋮& ⋮& \ddots & ⋮\\ 0_z& 0_z& 0_z& z_t\end{array}\right],$$

  where 0_z is a 1-by-(mp + r + 2) vector of zeros.

• $$\Lambda ={\left[\begin{array}{ccccccc}{\Phi }_1& {\Phi }_2& \cdots & {\Phi }_p& c& \delta & {\rm B}\end{array}\right]}^{\prime }$$, which is an (mp + r + 2)-by-m random matrix of the coefficients, and the m(mp + r + 2)-by-1 vector λ = vec(Λ).

• ε_t is an m-by-1 vector of random, serially uncorrelated, multivariate normal innovations with the zero vector for the mean and the m-by-m matrix Σ for the covariance. This assumption implies that the data likelihood is

  $$\ell \left(\Lambda ,\Sigma |y,x\right)={\displaystyle \prod _{t=1}^{T}f\left(y_t;\Lambda ,\Sigma ,z_t\right)},$$

  where f is the m-dimensional multivariate normal density with mean z_tΛ and covariance Σ, evaluated at y_t.

Before considering the data, you impose a joint prior distribution assumption on (Λ,Σ), which is governed by the distribution π(Λ,Σ). In a Bayesian analysis, the distribution of the parameters is updated with information about the parameters obtained from the data likelihood. The result is the joint posterior distribution π(Λ,Σ|Y,X,Y₀), where:

• Y is a T-by-m matrix containing the entire response series {y_t}, t = 1,…,T.

• X is a T-by-m matrix containing the entire exogenous series {x_t}, t = 1,…,T.

• Y₀ is a p-by-m matrix of presample data used to initialize the VAR model for estimation.

Tips

• Monte Carlo simulation is subject to variation. If simulate uses Monte Carlo simulation, then estimates and inferences might vary when you call simulate multiple times under seemingly equivalent conditions. To reproduce estimation results, set a random number seed by using rng before calling simulate.

Algorithms

• If simulate estimates a posterior distribution (when you supply Y) and the posterior is analytically tractable, simulate simulates directly from the posterior. Otherwise, simulate uses the Gibbs sampler to estimate the posterior.

• This figure shows how simulate reduces the sample by using the values of NumDraws, Thin, and BurnIn. Rectangles represent successive draws from the distribution. simulate removes the white rectangles from the sample. The remaining NumDraws black rectangles compose the sample.

  

• If PriorMdl is a semiconjugatebvarm object and you do not specify starting values (Coeff0 and Sigma0), simulate samples from the posterior distribution by applying the Gibbs sampler.

  1. simulate uses the default value of Sigma0 for Σ and draws a value of Λ from π(Λ|Σ,Y,X), the full conditional distribution of the VAR model coefficients.

  2. simulate draws a value of Σ from π(Σ|Λ,Y,X), the full conditional distribution of the innovations covariance matrix, by using the previously generated value of Λ.

  3. The function repeats steps 1 and 2 until convergence. To assess convergence, draw a trace plot of the sample.

  If you specify Coeff0, simulate draws a value of Σ from π(Σ|Λ,Y,X) to start the Gibbs sampler.

• simulate does not return default starting values that it generates.

Version History

Introduced in R2020a