2
$\begingroup$

enter image description here

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

$\endgroup$
3
$\begingroup$

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 
$\endgroup$
  • $\begingroup$ Thank You for your cooperation. Its means alot I am trying to understand and follow these steps. When I apply this command (arima) on my data, after applying it it gives an error "Error in optim(init[mask], armafn, method = optim.method, hessian = TRUE, : non-finite finite-difference value [1]" can you please explain wt does this error means? $\endgroup$ – user49456 Jul 5 '14 at 19:39
  • $\begingroup$ The error says that there was a problem in the optimization algorithm, apparently the numerical gradient was not a finite value. You can try: 1) pass different initial values through the argument init; 2) use method="CSS" in the function arima. If you post the data, the output of dput(data.object), I could have a look at it. $\endgroup$ – javlacalle Jul 5 '14 at 20:52
  • $\begingroup$ another query is in R these codes do not provide values of P-values and t-statistics, BIC, SIC, R-square etc $\endgroup$ – user49456 Jul 6 '14 at 17:36
  • 1
    $\begingroup$ You can obtain the AIC as AIC(fit, k = 2) and BIC as AIC(fit, k = log(length(x))); use logLik(fit) to extract the log-likelihood; t-statistics and p-values can be conveniently obtained using coeftest in package lmtest, require(lmtest); coeftest(fit); you may also be interested in confidence intervals for the parameter estimates, they can be obtained as confint(fit). The R-square is not that meaningful in the context of ARIMA models, see this post. $\endgroup$ – javlacalle Jul 6 '14 at 19:08

Your Answer

By clicking “Post Your Answer”, you agree to our terms of service, privacy policy and cookie policy

Not the answer you're looking for? Browse other questions tagged or ask your own question.