I am trying to simulate data to generate data series for 2 variables to do some further analysis. I am writing a for loop to generate the data series according to the data generating process I have in mind. There are 1006 observations in each simulation (1000 observations for each variable plus 3 lags for each variable) and the matrix A1 can take 100 values. The code does generate data series for the 2 variables but after a few observations (around the 15th or so), I start getting NaNs instead of the values. For both the variables, the values keep getting bigger and then eventually hit NaN. Can someone help me with this? I want full numerical values for further analysis rather than NaNs. Thanks.

```
N <- 100
A1_sequence <- seq(from=0.5,to=1.5,length.out = N)
A1 <- array(data = A1_sequence, dim=c(2,2,100))
A2 <- matrix(c(0.2,1.1,-0.6,0.2), nrow=2, ncol=2)
A3 <- matrix(c(-0.8,1.2,-0.4,0.4), nrow=2, ncol=2)
A4 <- matrix(c(-0.01, 0.02,-0.03,0.05), nrow=2, ncol=2)
p <- 3 # Number of lags
N1 <- 1000+2*p # Number of observations in each simulation
k <- 2 #Number of endogenous variables
x <- matrix(0, k, N1)
myeps1 <- 0.25*rnorm(k)
for( i in (p+1):N1){
for (j in 1:N){
x[,i] <- A1[,,j]%*%x[,i-1] - A2%*%x[,i-2] - A3%*%x[,i-3] - A4%*%((x[,i-1])^3) + myeps1
}
}
```