You can estimate an MA model using OLS, but you need to do it iteratively.
Consider an MA(1) without intercept:
$$
y_t = \epsilon_t + \theta\epsilon_{t-1}.
$$
The algorithm is:
Set an initial value, $\theta^{(0)}=0$. Fix $\hat{\epsilon}_0=0$. For $j=1, \dots, n_{iter}$ do:
- Compute the current estimate of the error terms,
$$
\hat{\epsilon}_t^{(j-1)}=y_t-\hat{\theta}^{(j-1)}\hat{\epsilon}_{t-1}^{(j-1)}
$$
- Run a regression of $y_t$ on $\hat{\epsilon}_{t-1}^{(j-1)}$ and update $\hat{\theta}^{(j)}$
- Repeat
Here is some R code and a figure showing how the estimate evolves in 10 iterations:
set.seed(1)
N <- 1000
nIter <- 10
eps <- rnorm(N+1)
theta <- 0.5
y <- eps[2:(N+1)] + theta * eps[1:N]
epsEst <- numeric(N+1)
epsEst[1] <- 0
thetaEst <- 0
thetaEstStorage <- numeric(nIter)
thetaEstStorage[1] <- thetaEst
for (j in 2:nIter) {
for (i in 1:N) {
epsEst[i+1] <- y[i] - thetaEst * epsEst[i]
}
mod <- lm(y ~ 0 + epsEst[1:N])
thetaEst <- coef(mod)
thetaEstStorage[j] <- thetaEst
}
ts.plot(thetaEstStorage, ylim = c(-1, 1))
abline(h = theta, col = "red")