\[ y_t = \phi_1 y_{t-1} + \phi_2 y_{t-2} + \phi_3 y_{t-3} + \phi_4 y_{t-4} + a_t \]

set.seed(123)
at <- rnorm(1000, mean = 0, sd = 2)
phi1 <- 0.5
phi2 <- -0.4
phi3 <- 0.3
phi4 <- -0.2
n <- length(at)
y <- numeric(n)

y[1] <- 0
y[2] <- 0
y[3] <- 0
y[4] <- 0

for (t in 5:n) {
  y[t] <- phi1 * y[t-1] +
          phi2 * y[t-2] +
          phi3 * y[t-3] +
          phi4 * y[t-4] +
          at[t]
}

y <- y[-(1:100)]
ts.plot(y)

par(mfrow=c(1,2))
acf(y)
pacf(y)