# 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)


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(2\pi) - 0.5 \log(\sigma^2) - 0.5 \log(|Q_n|) - \frac{y^\top Q_n^{-1}y}{2\sigma^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
# [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.

• thanx for reference.... – Praveen Apr 26 '16 at 6:39