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]
}