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 &amp;= \phi_1Z_{-1} + \cdots + \phi_pZ_{-p} + Z_0 \\
X_1 &amp;= \phi_1X_0 + \phi_2Z_{-1} + \cdots + \phi_pZ_{-(p-1)} \\
&amp;\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 &amp;= Z_1 \\
X_2 &amp;= \phi_1 X_1 + Z_2 \\
X_3 &amp;= \phi_1 X_2 + \phi_2 X_1 + Z_3 \\
&amp;\cdots
\end{align}
$$
Which approach is the "right" approach for constructing the simulation?
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).