a. Use the following R code to generate data from an AR(1) model with y1 = 0.

y <- numeric(100)
e <- rnorm(100)
for(i in 2:100)
  y[i] <- 0.6*y[i-1] + e[i]
sim <- tsibble(idx = seq_len(100), y = y, index = idx)

b. Produce a time plot for the series. How does the plot change as you change ϕ1?

# Simulate AR(1) with ϕ₁ = 0.9 (closer to non-stationary)
y_b <- numeric(100)
e_b <- rnorm(100)

for(i in 2:100) {
  y_b[i] <- 0.9 * y_b[i - 1] + e_b[i]
}

sim_b <- tsibble(idx = 1:100, y = y_b, index = idx)

# Plot to observe the effect of higher ϕ₁
ggplot(sim_b, aes(x = idx, y = y)) +
  geom_line() +
  ggtitle("b. AR(1) Model: ϕ₁ = 0.9")

c. Write your own code to generate data from an MA(1) model with θ1 = 0.6 and σ2 = 1.

# Simulate MA(1): y[t] = e[t] + 0.6 * e[t-1]
e_c <- rnorm(101)  # One extra to cover lag
y_c <- numeric(100)

for(i in 1:100) {
  y_c[i] <- e_c[i + 1] + 0.6 * e_c[i]
}

sim_c <- tsibble(idx = 1:100, y = y_c, index = idx)

# Plot the MA(1) model
ggplot(sim_c, aes(x = idx, y = y)) +
  geom_line() +
  ggtitle("c. MA(1) Model: θ₁ = 0.6")

d. Produce a time plot for the series. How does the plot change as you change θ1?

# Simulate MA(1) with θ₁ = -0.9
e_d <- rnorm(101)
y_d <- numeric(100)

for(i in 1:100) {
  y_d[i] <- e_d[i + 1] - 0.9 * e_d[i]
}

sim_d <- tsibble(idx = 1:100, y = y_d, index = idx)

# Plot the MA(1) model with negative θ₁
ggplot(sim_d, aes(x = idx, y = y)) +
  geom_line() +
  ggtitle("d. MA(1) Model: θ₁ = -0.9")

e. Generate data from an ARMA(1,1) model with ϕ1 = 0.6, θ1 = 0.6 and σ2 = 1.

# Simulate ARMA(1,1): y[t] = ϕ₁ * y[t-1] + e[t] + θ₁ * e[t-1]
e_e <- rnorm(101)
y_e <- numeric(100)

for(i in 2:100) {
  y_e[i] <- 0.6 * y_e[i - 1] + e_e[i] + 0.6 * e_e[i - 1]
}

sim_e <- tsibble(idx = 1:100, y = y_e, index = idx)

# Plot the ARMA(1,1) model
ggplot(sim_e, aes(x = idx, y = y)) +
  geom_line() +
  ggtitle("e. ARMA(1,1) Model: ϕ₁ = 0.6, θ₁ = 0.6")

f. Generate data from an AR(2) model with ϕ1 = −0.8, ϕ2 = 0.3 and σ2 = 1. (Note that these parameters will give a non-stationary series.)

# Simulate AR(2): y[t] = ϕ₁ * y[t-1] + ϕ₂ * y[t-2] + e[t]
e_f <- rnorm(100)
y_f <- numeric(100)
y_f[1:2] <- 0  # Initial values

for(i in 3:100) {
  y_f[i] <- -0.8 * y_f[i - 1] + 0.3 * y_f[i - 2] + e_f[i]
}

sim_f <- tsibble(idx = 1:100, y = y_f, index = idx)

# Plot the AR(2) model
ggplot(sim_f, aes(x = idx, y = y)) +
  geom_line() +
  ggtitle("f. AR(2) Model: ϕ₁ = −0.8, ϕ₂ = 0.3 (Non-stationary)")

g. Graph the latter two series and compare them.

# Compare the ARMA(1,1) and AR(2) models
p_arma <- ggplot(sim_e, aes(x = idx, y = y)) +
  geom_line() +
  ggtitle("ARMA(1,1): ϕ₁ = 0.6, θ₁ = 0.6")

p_ar2 <- ggplot(sim_f, aes(x = idx, y = y)) +
  geom_line() +
  ggtitle("AR(2): ϕ₁ = −0.8, ϕ₂ = 0.3 (Non-stationary)")

# Stack vertically using patchwork
p_arma / p_ar2

## The ARMA(1,1) series (ϕ₁ = 0.6, θ₁ = 0.6) displays stationary behavior, with values oscillating around a stable mean and showing short-term dependence. The series remains bounded, with no clear trend or drift. In contrast, the AR(2) series (ϕ₁ = −0.8, ϕ₂ = 0.3) appears non-stationary. It exhibits increasing variance over time and sometimes resembles an expanding cone. This suggests the presence of unit root behavior or instability, especially due to the parameter combination, where the roots of the characteristic polynomial lie on or near the unit circle.