Solved – Exact maximum likelihood estimation of MA(1)

It seems that the matrix L is not properly defined since the error solve.default(L) is returned. Looking at the elements that you are using, the approach that it seems you are seeking can be sketched and implemented as follows (the following is based on Pollock (1999) [1] Chapter 22):

The logarithm of the likelihood function of the following moving average (MA) process $$ y_t = varepsilon_t + theta varepsilon_{t-1} ,, quad t=1,dots,n ,, quad varepsilon sim NID(0, sigma^2) $$

is given by:

$$ L(y,sigma^2) = -0.5 n log(2pi) – 0.5 log(sigma^2) – 0.5 log(|Q_n|) – frac{y^top Q_n^{-1}y}{2sigma^2} ,, $$

where $Q_n$ is the matrix of theoretical autocovariances of the MA process. By concentrating $sigma^2$ out of the above expression, it can be seen that maximizing the log-likelihood function is equivalent to minimizing the following expression:

$$ S = frac{y^top Q_n^{-1}y}{n}|Q_n|^{1/n} ,. $$

Now, the expression for $S$ can be evaluated as follows:

Q <- toeplitz(ARMAacf(ma = x, lag.max = n - 1) * (1 + theta^2)) Qinv <- chol2inv(chol(Q)) tmp <- drop(y %*% Qinv %*% cbind(y)) / n Sval <- tmp * (det(Q))^(1/n) 

First, we take advantage of the function ARMAacf, which returns the theoretical autocorrelations. Since the variance of a MA(1) process is $1 + theta^2$, we can recover the autocovariances by multiplying those autocrrelations by $1 + theta^2$. Then, as $Q$ is a symmetric matrix, it is more efficient to compute its inverse upon the Cholesky decomposition. Finally, we put all the elements together to evaluate the value of objective function $S$.

You can use this code to compare with the elements that you are using in your implementation.

Let's test this code in a small simulation exercise:

# objective function introduced above Sfnc <- function(x, y, n) {   Q <- toeplitz(ARMAacf(ma = x, lag.max = n - 1) * (1 + x^2))   tmp <- drop(y %*% chol2inv(chol(Q)) %*% cbind(y)) / n   Sval <- tmp * (det(Q))^(1/n)   Sval } # parameters of the simulation niter <- 1000 # number of iterations n <- 120 # sample size ma <- 0.6 # MA coefficient of the data generating process mres <- matrix(nrow = niter, ncol = 2) # matrix for storage of results colnames(mres) <- c("Sfnc", "arima") # main loop set.seed(123) for (i in seq_len(niter)) {   # simulate a MA(1) of length 'n' and coefficient 'ma'   y <- arima.sim(n = n, model = list(ma = ma))   # minimize the objective function Sfnc   mres[i,1] <- optimize(f = Sfnc, interval = c(-1, 1), y = y, n = n)$min   # fit the same model using stats::arima   mres[i,2] <- coef(arima(y, order = c(0,0,1), include.mean = FALSE)) } 

As shown below, the results obtained with the method described above are in line with those estimates obtained by stats::arima:

# first parameter estimates for each method head(mres) # [1,] 0.7255807 0.7255885 # [2,] 0.6235342 0.6235284 # [3,] 0.6087424 0.6087536 # [4,] 0.4811465 0.4811588 # [5,] 0.5823097 0.5823167 # [6,] 0.6665756 0.6665669 # summary statistics for the estimates obtained with each method apply(mres, MARGIN = 2, summary) #           Sfnc  arima # Min.    0.3656 0.3656 # 1st Qu. 0.5484 0.5484 # Median  0.6026 0.6026 # Mean    0.6027 0.6027 # 3rd Qu. 0.6568 0.6568 # Max.    0.8664 0.8663 

I don't recommend this implementation as a substitute to stats::arima or other software packages, among other things because its reliability relies on the searching interval (-1,1) set in optimize; but it can be helpful for pedagogical purposes and some experimentation.

References

[1] D.S.G. Pollock (1999) A handbook of time series analysis, signal processing and dynamics, Academic Press.