How can I implement my own AR(p) simulation algorithm?
You question boils down to:


should you try to do the simulation in a vectorised manner?
should you iterate over a vector?


In general we should always try to vectorise (option 1) as it's much more efficient, and it is easy to lag vectors in R (embed) and Python Pandas (shift).

However, I can't see any efficient vectorisation-based solution. The autoregressive nature of the simulated values force us to iterate over a vector to simulate each of $X_1, X_2, \ldots, X_n$ individually (option 2).