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**hide

#### Best Answer

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.

### Similar Posts:

- Solved – Problem with the formulation of a gaussian copula likelihood function
- Solved – Efficient calculation of matrix inverse in R
- Solved – Prior selection for Gaussian Processes (GP)
- Solved – score function of bivariate/multivariate normal distribution
- Solved – score function of bivariate/multivariate normal distribution