Generating Error Vectors (White Noise) for Simulation of Vector Autoregressive Model (VAR)

I just estimated a Vector Autoregressive Model with 6 lags and 10 variables in R. My goal is to simulate the given original time series (on which the model parameters were estimated) to see how the model fits.

As the simulation in R doesn't work (but this is not the problem) I'd like to do the simulation with excel. I know how to use the estimated parameters and the previous values of the time series (respective the number of lags) but my problem is how to generate the right white noise for the single time series.

My first attempt was generating individual white noise vectors with mean=0 and sigma=squared Standard error of residual (got that from the R-Output of estimation) and it seems to work for some of the variables. But I have the feeling that this is not the right way to do it.

So, can you tell me how I generate the error terms for the simulation of the time series? Do I have to use the covariance matrix of residuals?

As the estimation and simulation of multivariate time series is just step 1 of 3 in my masters thesis (and time is running) I appreciate every hint I can get. All my books and papers couldn't help me with this question so I really hope some of you might have an idea.

EDIT: Thanks for the quick answer. Regarding the simulation with R: I tried the VAR.sim() command from the package tsDyn.

simulation<-VAR.sim(h, innov=rmnorm(n, mean=0, varcov=m), n=1095, lag=6, include="none", starting=v)


Explanation: I want to simulate a VAR(6) model and determined the matrix h ($10 \times 60$ matrix) with the estimated coefficients for all variables and lags. m is the covariance matrix of the residuals (also output from the estimation). I want to simulate 1095 days (n) without constant/trend (include="none") and give in v starting values for days 1 to 6 for all the 10 variables.

The hint for using the Command is:

VAR.sim(B, n=200, lag=1, include = c("const", "trend","none", "both"), starting=NULL, innov=rmnorm(n, mean=0, varcov=varcov), varcov=diag(1,nrow(B)), show.parMat=FALSE)


But the output doesn't fit to any of my original time series and I don't get what I'm doing wrong.

Do you have any idea what could be wrong? Or if there's another easy way to simulate a VAR(p)-model?

By the way, why doesn't the simulation in R work? You can simulate multivariate normal errors using the function mvrnorm from the MASS package.