# simulate

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

## Syntax

## Description

`[`

returns a random vector of coefficients `Coeff`

,`Sigma`

]
= simulate(`PriorMdl`

)`Coeff`

and a random innovations covariance matrix `Sigma`

drawn from the prior Bayesian VAR(*p*) model
`PriorMdl`

.

`[`

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.`Coeff`

,`Sigma`

]
= simulate(___,`Name,Value`

)

## 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 \text{\hspace{0.17em}}\mathit{V}\right)$.

$$\Sigma \sim Inverse\phantom{\rule{0.16666666666666666em}{0ex}}Wishart(\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` |

`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.

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:

Determine the extent of the transient behavior in the sample by setting the

`BurnIn`

name-value argument to`0`

.Simulate a few thousand observations by using

`simulate`

.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*`

)-by-`k`

`NumDraws`

elements, where

, which is the number of coefficients in a response equation. For details on the structure of * k* = PriorMdl.NumSeries*PriorMdl.P + PriorMdl.IncludeIntercept + PriorMdl.IncludeTrend + PriorMdl.NumPredictors

`Coeff0`

, see the output `Coeff`

.By default, `Coeff0`

is the multivariate least-squares estimate.

**Tip**

To specify

`Coeff0`

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

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].$$

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*`

)-by-`k`

`NumDraws`

numeric matrix, where

, which is the number of coefficients in a response equation. Each column is a separate draw from the distribution.* k* = PriorMdl.NumSeries*PriorMdl.P + PriorMdl.IncludeIntercept + PriorMdl.IncludeTrend + PriorMdl.NumPredictors

For draw

, `j`

`Coeff(1:`

corresponds to all coefficients in the equation of response variable * k*,

*)*

`j`

`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

`(`

through– 1)*PriorMdl.NumSeries + 1`q`

correspond to the lag*PriorMdl.NumSeries`q`

AR coefficients of the response variables ordered by`q`

`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

compose the vector of regression coefficients of the exogenous variables.`k`

This figure shows the structure of `Coeff(L,`

for a 2-D VAR(3) model that contains a constant vector and four exogenous predictors.* j*)

$$[\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}}}\text{\hspace{0.17em}}\text{\hspace{0.17em}}\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*is the_{t}*m*-dimensional observed response vector, where*m*=`numseries`

.Φ

_{1},…,Φ_{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}\\ \vdots & \vdots & \ddots & \vdots \\ {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(Λ).*ε*is an_{t}*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*_{0}), 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*_{0}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.`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.`simulate`

draws a value of Σ from*π*(Σ|Λ,*Y*,*X*), the full conditional distribution of the innovations covariance matrix, by using the previously generated value of Λ.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**

## See Also

### Objects

### Functions

## Open Example

You have a modified version of this example. Do you want to open this example with your edits?

## MATLAB Command

You clicked a link that corresponds to this MATLAB command:

Run the command by entering it in the MATLAB Command Window. Web browsers do not support MATLAB commands.

# Select a Web Site

Choose a web site to get translated content where available and see local events and offers. Based on your location, we recommend that you select: .

You can also select a web site from the following list:

## How to Get Best Site Performance

Select the China site (in Chinese or English) for best site performance. Other MathWorks country sites are not optimized for visits from your location.

### Americas

- América Latina (Español)
- Canada (English)
- United States (English)

### Europe

- Belgium (English)
- Denmark (English)
- Deutschland (Deutsch)
- España (Español)
- Finland (English)
- France (Français)
- Ireland (English)
- Italia (Italiano)
- Luxembourg (English)

- Netherlands (English)
- Norway (English)
- Österreich (Deutsch)
- Portugal (English)
- Sweden (English)
- Switzerland
- United Kingdom (English)