# Maximum Likelihood estimation for ARMA(1,1) in R

I have written the following code in R to estimate the parameters of a ARMA(1,1) process.

armacoeff <- function(x) {

l     = length(x)
param = c(mu=0, phi=0, theta=0)

SSE <- function(param) {
mu    = param[1]
phi   = param[2]
theta = param[3]

res    = vector()
res[1] = 0
for(i in (2:l)) {
res[i] = z[i] - (mu+z[i-1]*phi) - (res[i-1]*theta)
}

return(sum(res*res))
}

return(nlminb(objective=SSE, start= param))
}


Now, as far as I understand this code should give me the Maximum Likelihood Estimates for $$\mu, \phi$$ and $$\theta$$ but they do not align with the estimates given from the arima function.

Namely, the AR1 estimate from arima corresponds to $$\theta$$ and the MA1 estimate corresponds to $$\phi$$. According to my derived likelihood function this should not be the case. Can anyone point out my errors?

I have attached the following results for a example time series called "test"

ARIMA estimate

arima(test, order=c(1,0,1))

Call:
arima(x = test, order = c(1, 0, 1))

Coefficients:
ar1     ma1  intercept
-0.0735  0.1030      1e-04
s.e.   0.2799  0.2815      4e-04

sigma^2 estimated as 0.0005476:  log likelihood = 8311.68,  aic = -16615.36


And now the result for armacoeff

 h=armacoeff(test)
> h
$par mu phi theta 1.944046e-05 9.743507e-02 -7.261513e-02$objective
[1] 1.943927

$convergence [1] 0$iterations
[1] 11

• You are copy-pasting your code and asking us to fnd the bug. You need to show that you have done some research effort before asking for help here. Besides this, the question is off-topic as it is about debagging. See stackoverflow.com for this. – utobi Nov 28 '16 at 5:58
• @utobi, Code Review Stack Exchange could be a better fit than Stack Overflow for debugging problems. However, this one needs quite some expertise in statistics (rather than programming) to be answered, so I don't think better help could be found elsewhere. But in present shape it could indeed be considered off topic here, IMHO. – Richard Hardy Nov 28 '16 at 9:45
• I appreciate your comments. I got the answer I was looking for. My problem was realizing the difference between my calculated MLE and the one realized by arima. I just used the code as an example. As @RobHyndmann pointed out this was due to the unconditional vs conditional likelihood. – T-at-R Nov 28 '16 at 10:47

First, you are optimizing the conditional likelihood as you are conditioning on the first observation. arima optimizes the full likelihood by default.

Second, your parameterization of the intercept is different.

Here is some code and output which match. I've also corrected your armacoeff function as you had z in place of x inside the function.

armacoeff <- function(x) {

l = length(x)
param=c(mu=0, phi=0, theta=0)

SSE <- function(param){
mu=param[1]
phi=param[2]
theta=param[3]

res = vector()
res[1] = 0
for(i in (2:l)){
res[i] = (x[i]-mu) - phi*(x[i-1]-mu) - theta*res[i-1]
}
return(sum(res*res))
}

bla =nlminb(objective=SSE, start= param)
return(bla)

}

set.seed(30)
test <- rnorm(20)

arima(test, order=c(1,0,1), method="CSS")
armacoeff(test)


> arima(test, order=c(1,0,1), method="CSS")

Call:
arima(x = test, order = c(1, 0, 1), method = "CSS")

Coefficients:
ar1     ma1  intercept
-0.2324  0.5670    -0.3746
s.e.   0.3153  0.2815     0.2571

sigma^2 estimated as 0.8492:  part log likelihood = -26.74
> armacoeff(test)
$par mu phi theta -0.3746246 -0.2324023 0.5669316$objective
[1] 16.13464

$convergence [1] 0$iterations
[1] 13

$evaluations function gradient 18 50$message
[1] "relative convergence (4)"

• Thank you very much. So there is nothing fundamentally wrong with my code but with my understanding of the ML estimation of ARMA models. – T-at-R Nov 28 '16 at 7:35
• From your answer I assume that one is more interested in the unconditional likelihood, am I right? – T-at-R Nov 28 '16 at 7:36