Simulate 100 AR(2) time series

library(FitAR)
set.seed(54321)
n=50
phi <- c(0.1,0.5)

count <- 0

for(i in 1:100){
  yt <- unclass(arima.sim(n=n,list(ar=phi),innov=rnorm(n,0,1)))
  p=SelectModel(as.ts(yt), lag.max = 20, Criterion = "BIC", Best=1)
  fit.monthly <- arima(yt, order = c(p, 0, 0))
  my_coefficients =fit.monthly$coef
  my_coefficients=my_coefficients[!names(my_coefficients) == 'intercept']
  print(my_coefficients)

  if(length(my_coefficients) == 2){
    count <- count + 1
  }
}

print(paste0("AR(2) model count is: ", count))

Use the rGARMA function in the ts.extend package

You can generate random vectors from any stationary Gaussian ARMA model using the ts.extend package. This package generates random vectors directly form the multivariate normal distribution using the computed autocorrelation matrix for the random vector, so it gives random vectors from the exact distribution and does not require "burn-in" iterations. Here is an example from an AR(2) model.

#Load the package
library(ts.extend)

#Set parameters
AR <- c(0.9, -0.2)
m  <- 50

#Generate n = 100 random vectors from this model
set.seed(1)
SERIES <- rGARMA(n = 100, m = m, ar = AR, errorvar = 1)

#Plot the series using ggplot2 graphics
library(ggplot2)
plot(SERIES)