How can I implement my own AR(p) simulation algorithm?

Recall that an AR(p) process can be defined as $$X_t = \phi_1 X_{t-1} + \cdots + \phi_p X_{t-p} + Z_t$$ where $Z_t$ is I.I.D. white noise. I want to write an algorithm simulating these processes, but I'm not sure if I'm taking the right approach. Here's my idea:

INPUT:
n = length of time series
[\phi_1,...,\phi_p] the parameters of the process

1. Construct a white noise vector Z of length n + p
2. Have the first rolling sum of the X_i's be set as the noise
3. Construct the rest of the time series from there


In this case \begin{align} X_0 &= \phi_1Z_{-1} + \cdots + \phi_pZ_{-p} + Z_0 \\ X_1 &= \phi_1X_0 + \phi_2Z_{-1} + \cdots + \phi_pZ_{-(p-1)} \\ &\cdots \end{align}

Another option would be to start by constructing $X_1,\ldots, X_p$ by having the lower order terms in the recursion equation be set to 0. So, for the first few \begin{align} X_1 &= Z_1 \\ X_2 &= \phi_1 X_1 + Z_2 \\ X_3 &= \phi_1 X_2 + \phi_2 X_1 + Z_3 \\ &\cdots \end{align} Which approach is the "right" approach for constructing the simulation?

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).