# Solved – Exact maximum likelihood estimation of MA(1)

I want to calculate the MLEs of the MA(1) model and for this purpose I have written the exact likelihood for the same. I built a programme in R for the log-likelihood, but it seems some problem in it whenever I use `nlm` function to optimize this. Here is my code in R.

``rm(list = ls()) y <- c(2.33,2.84,6.09,4.25,2.56,5.69,7.59,2.19,7.08,3.1) y=as.matrix(y) T=length(y) #psi=c(0.0054,0.045) #psi=as.matrix(psi) #q=length(psi) #sigma=9.6 q=1  loglike=function(parm){  psi=parm[1];sigma=parm[2];  L=diag(rep(1,T)); for(i in 1:T){ if(i>q){psi[i]=0} for(j in 1:T-1){ if((i+j)>T)break(); L[i+j,j]=psi[i] } } L  B=diag(rep(psi[q]),q); for(i in 1:q-1){ for(j in 1:q){ if((i+j)>q)break(); B[i,i+j]=psi[q-j] } } B  O=matrix(0,T-q,q); F=t(cbind(t(B),t(O)))  I=diag(rep(1),q) D=I+(t(F))%*%(t(solve(L)))%*%(solve(L))%*%F ehat=(solve(D))%*%(t(F))%*%(t(solve(L)))%*%(solve(L))%*%y S=(t(y-F%*%ehat))%*%(t(solve(L)))%*%(solve(L))%*%(y-F%*%ehat)+(t(ehat))%*%ehat  loglike=-(T/2)*(log(6.28)+log(sigma))-(1/2)*solve(det(D))-(1/2*sigma)*S  return(-loglike)  }  nlm(loglike,c(0.0054,9.6),hessian=T) ``
Contents

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.

Rate this post