This shows how to simulate a VAR with GARCH (DVECH) errors. A simpler example without the VAR mean model can be found at GARCHMVSIMULATE.RPF. An example with a univariate (AR-GARCH) process can be found at ARGARCHSIM.RPF.

This first simulates a mean zero GARCH error process to form the shocks, then simulates the VAR taking the GARCH shocks as input.

This sets the number of observations and number of variables:

compute nobs=500

compute n=2

This sets up the control parameters for the GARCH processâ€”as a DVECH process, the C, A and B are all $$n \times n$$SYMMETRIC matrices, so only the lower triangle needs to be specified. If you want a process larger than 2 variables, you'll need to expand these matrices to the proper dimension:

dec symm vc(n,n) va(n,n) vb(n,n)

compute vc=||.05|.01,.05||

compute va=||.08|.05,.08||

compute vb=||.90|.85,.90||

This computes the stationary covariance matrix. Because this is a DVECH, this can be done element by element.

dec symm h0(n,n)

ewise h0(i,j)=vc(i,j)/(1-va(i,j)-vb(i,j))

This creates the formulas which simulate the mean zero DVECH process. HF will compute the variance given the lagged outer product of residuals (UU) and lagged variance (H); HU (in order)

1. computes the covariance matrix
2. does mean zero normal draws with the just computed covariance matrix using %RANMVNORMAL
3. computes and saves the outer product of the residuals using %OUTERXX
4. saves the residuals into the elements of the VEC[SERIES] U using %PT
5. returns the covariance matrix as the final value of the FRML

frml hf   = vc+va.*uu{1}+vb.*h{1}

frml hu   = hx=hf,ux=%ranmvnormal(%decomp(hx)),uu=%outerxx(ux),%pt(u,t,ux),hx

This initializes the outer product and variance series to the stationary variance of the process, then simulates the process from entries 2 on.

gset uu   = h0

gset h    = h0

*

* Generate starting in entry 2. This creates the simulated GARCH process

* as a side effect of generating the sequence of covariance matrices.

*

gset h 2 nobs = hu(t)

If we just wanted a multivariate GARCH process, we would now be done, with U(1) and U(2) as the two series. The following creates a 2 lag VAR which YGARCH(1) and YGARCH(2) as the two series. The two processes will be

$$\begin{array}{*{20}{l}}{{y_{1,t}} = .80{y_{1,t - 1}} + .15{y_{1,t - 2}} + .30{y_{2,t - 1}} - .50{y_{2,t - 2}} + {\varepsilon _{1,t}}} \\ {{y_{2,t}} = - .20{y_{1,t - 1}} + .20{y_{1,t - 2}} + .85{y_{2,t - 1}} + .10{y_{2,t - 2}} + {\varepsilon _{2,t}}} \\\end{array}$$

dec vect[series] ygarch(2)

system(model=vargarch)

variables ygarch

lags 1 2

end(system)

*

dec rect varcoeffs(4,2)

input varcoeffs

.80 -.20

.15  .20

.30  .85

-.50  .10

compute %modelsetcoeffs(vargarch,varcoeffs)

Just in case, this checks whether the largest root isn't explosive (that it's less than 1):

disp "Largest root of VAR" %modellargestroot(vargarch)

This does the simulation of the process, using U as the input shocks (the PATHS option on FORECAST indicates that it will be forecast with full paths of shocks):

clear(zeros) ygarch

forecast(paths,from=3,to=nobs,model=vargarch,results=ygarch)

# u

For illustration, the program estimates the VAR by least squares and by the (properly specified) GARCH model and computes forecasts (of the mean of the data) for each. It uses @UFOREERRORS to compute statistics on the forecast errors. Note that this is a very bad idea. The standard forecast performance statistics are based upon an assumption of heteroscedasticityâ€”the GARCH model (intentionally) downweights the residuals in the high variance zones. As a result, least squares will be expected to look somewhat better by these measures than the properly specified GARCH model.

Full Program

compute nobs=500
compute n=2
*
all nobs
*
* u is the simulated GARCH process
*
dec vect[series] u(n)
*
* H is the series of variance matrices, UU is the series of lagged outer products of the residuals.
*
dec series[symm] h uu
dec symm vc(n,n) va(n,n) vb(n,n)
compute vc=||.05|.01,.05||
compute va=||.08|.05,.08||
compute vb=||.90|.85,.90||
*
* Compute the stationary variance
*
dec symm h0(n,n)
ewise h0(i,j)=vc(i,j)/(1-va(i,j)-vb(i,j))
*
dec frml[symm] hf hu
dec symm hx
dec vect ux
*
frml hf   = vc+va.*uu{1}+vb.*h{1}
frml hu   = hx=hf,ux=%ranmvnormal(%decomp(hx)),uu=%outerxx(ux),%pt(u,t,ux),hx
*
* Initialize the UU and H series for the pre-sample values
*
gset uu   = h0
gset h    = h0
*
* Generate starting in entry 2. This creates the simulated GARCH process
* as a side effect of generating the sequence of covariance matrices.
*
gset h 2 nobs = hu(t)
******************************************************************
*
* Create a VAR with the GARCH as an error process
*
dec vect[series] ygarch(2)
system(model=vargarch)
variables ygarch
lags 1 2
end(system)
*
dec rect varcoeffs(4,2)
input varcoeffs
.80 -.20
.15  .20
.30  .85
-.50  .10
compute %modelsetcoeffs(vargarch,varcoeffs)
*
* Make sure the dominant root isn't explosive
*
disp %modellargestroot(vargarch)
*
clear(zeros) ygarch
forecast(paths,from=3,to=nobs,model=vargarch,results=ygarch)
# u
*
* Discard the first 100 entries in doing the estimates to give the VAR process
* time to overcome the initial conditions.
*
* Estimate the model by OLS and do out-of-sample forecasts of the mean.
*
estimate(model=vargarch) 101 nobs
disp %logl
forecast(static,model=vargarch,from=%regstart(),to=%regend(),results=olsforecasts)
*
* Estimate the model by GARCH and do out-of-sample forecasts of the mean.
*
garch(p=1,q=1,model=vargarch,presample=%sigma) 101 nobs ygarch
forecast(static,model=vargarch,from=%regstart(),to=%regend(),results=garchforecasts)
*
* Do forecast error statistics on the two sets of forecasts. This is to
* demonstrate that this is a *really* bad idea for evaluating a GARCH
* model. Standard forecast statistics (RMSE, MAE) are based upon equal
* weights on the errors. A GARCH model (intentionally) downweights the
* errors on the high-variance parts of the sample, while OLS doesn't.
* Thus OLS would be expected to do no worse (and should theoretically be
* slightly better) even though, in this case, we know that the GARCH is
* correctly specified and OLS isn't.
*
@uforeerrors(title="Forecast Errors from GARCH Model, Variable 1") ygarch(1) garchforecasts(1)
@uforeerrors(title="Forecast Errors from Least Squares, Variable 1")  ygarch(1) olsforecasts(1)
@uforeerrors(title="Forecast Errors from GARCH Model, Variable 2") ygarch(2) garchforecasts(2)
@uforeerrors(title="Forecast Errors from Least Squares, Variable 2")  ygarch(2) olsforecasts(2)


Output

This depends upon random numbers so won't make exactly.

Largest root of VAR       0.95001

VAR/System - Estimation by Least Squares

Usable Observations                       400

Dependent Variable YGARCH(1)

Mean of Dependent Variable       13.242584903

Std Error of Dependent Variable  18.335985914

Standard Error of Estimate        1.436261677

Sum of Squared Residuals         816.88765143

Durbin-Watson Statistic                2.0319

Variable                        Coeff      Std Error      T-Stat      Signif

************************************************************************************

1.  YGARCH(1){1}                  0.793498054  0.045043430     17.61629  0.00000000

2.  YGARCH(1){2}                  0.155062033  0.043628417      3.55415  0.00042497

3.  YGARCH(2){1}                  0.224214420  0.045841859      4.89104  0.00000146

4.  YGARCH(2){2}                 -0.446576268  0.047950196     -9.31334  0.00000000

Dependent Variable YGARCH(2)

Mean of Dependent Variable       -2.884887297

Std Error of Dependent Variable   5.956208544

Standard Error of Estimate        1.522227477

Sum of Squared Residuals         917.60189102

Durbin-Watson Statistic                2.0186

Variable                        Coeff      Std Error      T-Stat      Signif

************************************************************************************

1.  YGARCH(1){1}                 -0.218659572  0.047739453     -4.58027  0.00000623

2.  YGARCH(1){2}                  0.216467336  0.046239746      4.68141  0.00000392

3.  YGARCH(2){1}                  0.859444947  0.048585671     17.68927  0.00000000

4.  YGARCH(2){2}                  0.093811273  0.050820200      1.84594  0.06564643

Log likelihood of OLS VAR   -1443.64784

MV-GARCH - Estimation by BFGS

Convergence in   136 Iterations. Final criterion was  0.0000000 <=  0.0000100

Usable Observations                       400

Log Likelihood                     -1422.3435

Variable                        Coeff      Std Error      T-Stat      Signif

************************************************************************************

Mean Model(YGARCH(1))

1.  YGARCH(1){1}                  0.794186519  0.046840609     16.95509  0.00000000

2.  YGARCH(1){2}                  0.155187638  0.045416674      3.41697  0.00063321

3.  YGARCH(2){1}                  0.211788446  0.044428845      4.76691  0.00000187

4.  YGARCH(2){2}                 -0.432950739  0.046639755     -9.28287  0.00000000

Mean Model(YGARCH(2))

5.  YGARCH(1){1}                 -0.243085415  0.037981481     -6.40010  0.00000000

6.  YGARCH(1){2}                  0.239604031  0.036832921      6.50516  0.00000000

7.  YGARCH(2){1}                  0.867813537  0.047891438     18.12043  0.00000000

8.  YGARCH(2){2}                  0.077213525  0.049434715      1.56193  0.11830466

9.  C(1,1)                        0.106172309  0.156402252      0.67884  0.49723842

10. C(2,1)                        0.002298371  0.011572317      0.19861  0.84256833

11. C(2,2)                        0.071408205  0.052429489      1.36199  0.17320245

12. A(1,1)                        0.047823722  0.047332958      1.01037  0.31231885

13. A(2,1)                       -0.007751368  0.028395583     -0.27298  0.78487018

14. A(2,2)                        0.123587284  0.040303313      3.06643  0.00216632

15. B(1,1)                        0.900149358  0.116376873      7.73478  0.00000000

16. B(2,1)                        0.930626910  0.257317756      3.61664  0.00029845

17. B(2,2)                        0.852228747  0.048781865     17.47020  0.00000000

Forecast Errors from GARCH Model, Variable 1

From 101 to 500

Mean Error             0.02839049

Mean Absolute Error    1.10026909

Root Mean Square Error 1.42925210

Mean Square Error        2.042762

Theil's U                0.812777

Forecast Errors from Least Squares, Variable 1

From 101 to 500

Mean Error             0.03584847

Mean Absolute Error    1.09951033

Root Mean Square Error 1.42906232

Mean Square Error        2.042219

Theil's U                0.812669

Forecast Errors from GARCH Model, Variable 2

From 101 to 500

Mean Error             -0.1202795

Mean Absolute Error     1.2029680

Root Mean Square Error  1.5153372

Mean Square Error        2.296247

Theil's U                0.959166

Forecast Errors from Least Squares, Variable 2

From 101 to 500

Mean Error             -0.1134268

Mean Absolute Error     1.2000395

Root Mean Square Error  1.5145972

Mean Square Error        2.294005

Theil's U                0.958698