Fitting an ARIMA model with lme() (package nlme)

An ARIMA model is one way to model time-series with two components, autoregressive and mobile average. See this post Modèle ARMA et ARIMA.

Let create a time series using both AR1 and MA1:

set.seed(101)
d <- data.frame(y=10,x=1:300)
epsilon <- arima.sim(model = list(ar=0.9, ma=+0.2279), n = 300)

d$y <- d$y + epsilon
layout(matrix(1))
plot(d$x, d$y, bty="n", las=1, pch=19, cex=0.5)

The parameters describing this time series can be obtained using different functions. The order parameter is c(number of AR, number of diff, number of MA) for arima() and FitARMA() and c(number of AR, number of MA) for arma():

library(stats)
ai_ <- arima(d$y, order=c(1, 0, 1), include.mean = TRUE)
library(tseries)
ai2_ <- arma(d$y, order = c(1, 1), include.intercept = TRUE)
library(FitARMA)
ai3_ <- FitARMA(d$y, order = c(1, 0, 1), demean = TRUE)

coef(ai_)
coef(ai2_)
coef(ai3_)

1/ Note that the intercept for arma() and FitArma() and arima() are different because the fitted model is different:
tseries::arma "y[t] = a[0] + a[1]y[t-1] + …"  vs stats::arima "(X[t]-m) = a[1](X[t-1]-m) + …" (from an answer of Adrian Trapletti <adrian@trapletti.org>)
2/ Note also that estimates of parameters are wrong if intercept is set to FALSE.
3/ The sign of MA can be positive or negative depending how the model is parametrized.

This arima model can be used also as an auto correlative structure for fitting. All these codes give the same result:


set.seed(101)

d <- data.frame(y=10,x=rep(1:30, 10),f=factor(sort(rep(1:10,30))))
epsilon <- arima.sim(model = list(ar=0.9, ma=+0.2279), n = 300)


d$y <- d$y + epsilon

m1 <- lme(y~x,random=~1|f,correlation=corARMA(form= ~ 1 | f, p=1, q=1),data=d)
summary(m1)
m1$modelStruct$corStruct
coef(m1$modelStruct$corStruct, unconstrained=FALSE)

m2 <- lme(y~x,random=~1|f,correlation=corARMA(p=1, q=1),data=d)
summary(m2)
coef(m2)
m2$modelStruct$corStruct
coef(m2$modelStruct$corStruct, unconstrained=FALSE)

cs1ARMA <- corARMA(c(0.95, 0.20), form = ~1 | f, p = 1, q = 1)
cs1ARMA. <- Initialize(cs1ARMA, data = d)
m3 <- lme(y~x,random=~1|f,correlation=cs1ARMA.,data=d)
summary(m3)
m3$modelStruct$corStruct
coef(m3$modelStruct$corStruct, unconstrained=FALSE)

cs1ARMA <- corARMA(c(0.95, 0.20), p = 1, q = 1)
cs1ARMA. <- Initialize(cs1ARMA, data = d)
m4 <- lme(y~x,random=~1|f,correlation=cs1ARMA.,data=d)
summary(m4)
m4$modelStruct$corStruct
coef(m4$modelStruct$corStruct, unconstrained=FALSE)