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

Question

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

Answer 1

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