set.seed(666)

N <- 1000

noise <- rnorm(N, 0, 1)

phi.1 <- .4
theta.1 <- -.3
theta.2 <- -.8

Z <- noise[1]

Z[2] <- phi.1*Z[1] + noise[2] + theta.1*noise[1]

# t >= 3
for (t in 3:N) {
  Z[t] = phi.1*Z[t-1] + noise[t] + theta.1*noise[t-1] + theta.2*noise[t-2]
}