Time series estimation on specific lags in ARMA model

Question

How can I fit this ARMA model having specific lags using R?

Answer 1

You can use the argument fixed in function arima and fix some coefficients to zero. Only those elements in fixed that are NA values will be updated by the optimization algorithm. An example with your model:

x <- diff(diff(log(AirPassengers)), 12)
fit <- arima(x, order = c(10,0,13), include.mean = FALSE, 
  fixed = c(NA, 0, NA, rep(0, 6), NA, NA, rep(0, 11), NA),
  transform.pars = FALSE)
fit
#Coefficients:
#          ar1  ar2      ar3  ar4  ar5  ar6  ar7  ar8  ar9     ar10      ma1
#      -0.0070    0  -0.2016    0    0    0    0    0    0  -0.0514  -0.3583
#s.e.   0.2402    0   0.0902    0    0    0    0    0    0   0.1012   0.2297
#      ma2  ma3  ma4  ma5  ma6  ma7  ma8  ma9  ma10  ma11  ma12     ma13
#        0    0    0    0    0    0    0    0     0     0     0  -0.0219
#s.e.    0    0    0    0    0    0    0    0     0     0     0   0.0974
#sigma^2 estimated as 0.001751:  log likelihood = 229.74,  aic = -447.48

Be aware that transform.pars=FALSE is required when some of the parameters are fixed. You may want to check that the roots of the estimated AR polynomial lie in the region of stationarity and to enforce an invertible MA polynomial. For that, you can use the following functions that are defined inside the function arima.

arCheck <- function(ar)
{
    p <- max(which(c(1, -ar) != 0)) - 1
    if(!p) return(TRUE)
    all(Mod(polyroot(c(1, -ar[1L:p]))) > 1)
}
arcoefs <- coef(fit)[grepl("^ar\\d+$", names(coef(fit)))]
arCheck(arcoefs)
#[1] TRUE

maInvert <- function(ma)
{
    ## polyroot can't cope with leading zero.
    q <- length(ma)
    q0 <- max(which(c(1,ma) != 0)) - 1L
    if(!q0) return(ma)
    roots <- polyroot(c(1, ma[1L:q0]))
    ind <- Mod(roots) < 1
    if(all(!ind)) return(ma)
    if(q0 == 1) return(c(1/ma[1L], rep(0, q - q0)))
    roots[ind] <- 1/roots[ind]
    x <- 1
    for (r in roots) x <- c(x, 0) - c(0, x)/r
    c(Re(x[-1L]), rep(0, q - q0))
}
macoefs <- coef(fit)[grepl("^ma\\d+$", names(coef(fit)))]
maInvert(macoefs) # no transformation required in this case
#        ma1         ma2         ma3         ma4         ma5         ma6 
#-0.35830514  0.00000000  0.00000000  0.00000000  0.00000000  0.00000000 
#        ma7         ma8         ma9        ma10        ma11        ma12 
# 0.00000000  0.00000000  0.00000000  0.00000000  0.00000000  0.00000000 
#       ma13 
#-0.02190729