set.seed(47)
n = 200

a = rnorm(n,0,1)
a0 = rnorm(1,0,1)
an1 = rnorm(1,0,1)

theta = c(0.2, 0.5)

y = NULL
y[1] = a[1] + theta[1]*a0 + theta[2]*an1
y[2] = a[2] + theta[1]*a[1] + theta[2]*a0
for (t in 3:n) {
  y[t] = a[t] + theta[1]*a[t-1] + theta[2]*a[t-2]
}

MAfit = arima(y,order = c(0,0,2))

residuals(MAfit)[1]      #  (returns 2.745251)
coef(MAfit)[3]     # (returns 3.122668)

plot(y)

hist(residuals(MAfit))
coef(MAfit)

###

n <- 50 
y2 <- 1:50
a = c(0, 0.5, 0.3)
y2[1] <- 0
y2[2] <- 1

#y = NULL
#a[1] + theta[1]*a0 + theta[2]*an1
#y[2] = a[2] + theta[1]*a[1] + theta[2]*a0
for (t in 3:n) { 
 epsilon <- rnorm(1,0,2)
  y2[t] = a[1] + a[2]*y2[t-1] + a[3]*y2[t-2] + epsilon
}

acf(y2)
pacf(y2)

ARfit = arima(y2,order = c(2,0,0))

residuals(ARfit)[1]      #  (returns 2.745251)
coef(ARfit)[3]     # (returns 3.122668)

plot(y2)

hist(residuals(ARfit))
coef(ARfit)

summary(ARfit)
ARfit$sigma2