BEKK model simulation in R

I have been working with a BEKK(1,1) model with dimension 3,4, and 5, for a time series analysis. I was given the feedback to include a simulation study. In order to trust the results that I obtain, I want to, via simulations, show that the estimation of the BEKK model parameters works also well for the sample sizes considered in the paper. I want to show that the distributional theory can be applied for my sample size.

I want to investigate if the sample size is enough to apply the asymptotic results?

Method:

I wish to generate data from the fitted model in the case of the dimension 3, with the sample size 3000. Estimate parameter by fitting BEKK model to the generated data sets and repeat this step, say 10000 times. Then I obtain 10000 estimators for each parameter for which the sampling distribution can be constructed and then it has to be compared with the asymptotic distribution.

Then repeat this procedure for dimension 4 and 5.

#I've been using the mgarchBEKK package when creating my BEKK models.
#The package provides the example below as help for simulation:

## Simulate series:
simulated <- simulateBEKK(2, 1000, c(1,1))

## Prepare the matrix:
simulated <- do.call(cbind, simulated$eps)

## Estimate with default arguments:
estimated <- BEKK(simulated)

I'm not a master in R by any means. So I'm not quite sure how to code the procedure that I describe above.

Any help is greatly apprecieated :)