You can make x an array and store the N different results. I do this in the code below and show the results for i = 1:20 for j = 1 and j = 100. Since the values of A1 are smaller for j = 1 than j = 100, the calculated values increase a little more slowly for j = 1 but the i = 20 value is about 1E+190, so the next step will be NaN even for j = 1
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 <- array(data = 0, c(k, N1, N))
myeps1 <- 0.25*rnorm(k)
for( i in (p+1):N1){
for (j in 1:N){
x[,i,j] <- A1[,,j]%*%x[,i-1,j] - A2%*%x[,i-2,j] - A3%*%x[,i-3,j] - A4%*%((x[,i-1,j])^3) + myeps1
}
}
x[,1:20,1]
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
#> [1,] 0 0 0 0.03903182 -0.1454364 -0.5779590 -1.0771888 -1.52723046
#> [2,] 0 0 0 -0.38873817 -0.5720410 -0.7220469 -0.6607307 -0.06567585
#> [,9] [,10] [,11] [,12] [,13] [,14] [,15]
#> [1,] -1.726565 -1.080278 1.190540 4.031453 6.4058654 8.372963 -26.71532
#> [2,] 1.168073 2.623233 3.098518 2.608855 -0.8293063 -10.414098 30.33007
#> [,16] [,17] [,18] [,19] [,20]
#> [1,] 645.6901 -29721583 3.189299e+21 -3.962298e+63 7.600320e+189
#> [2,] -1026.1268 48638138 -5.227981e+21 6.495696e+63 -1.245984e+190
x[,1:20, 100]
#> [,1] [,2] [,3] [,4] [,5] [,6] [,7] [,8]
#> [1,] 0 0 0 0.03903182 -0.4845456 -2.295505 -8.015300 -32.351066
#> [2,] 0 0 0 -0.38873817 -0.9111502 -2.397584 -5.627079 2.467767
#> [,9] [,10] [,11] [,12] [,13] [,14] [,15]
#> [1,] -386.5318 7492132 -5.152917e+19 1.671953e+58 -5.710498e+173 NaN NaN
#> [2,] 645.5113 -12293287 8.448013e+19 -2.740981e+58 9.361703e+173 NaN NaN
#> [,16] [,17] [,18] [,19] [,20]
#> [1,] NaN NaN NaN NaN NaN
#> [2,] NaN NaN NaN NaN NaN
Created on 2022-01-03 by the reprex package (v2.0.1)