Modelización ARMA en R

Question

Tengo una serie temporal. Quiero modelarla utilizando ARMA, que se utilizará para la previsión.

En R estoy usando arima() para obtener los coeficientes. Pero arima() requiere orden(p,d,q) como entrada. ¿Cuál es la forma más sencilla en R de llegar a un buen valor para p y q (con d = 0) para no sobreajustar?

Answer 1

Una opción es ajustar una serie de modelos ARMA con combinaciones de $p$ y $q$ y trabajar con el modelo que tenga el mejor "ajuste". Aquí evalúo el "ajuste" utilizando el BIC para intentar penalizar los ajustes demasiado complejos. A continuación se muestra un ejemplo para el Mauna Loa incorporado $\mathrm{CO}_2$ conjunto de datos de concentración

## load the data
data(co2)
## take only data up to end of 1990 - predict for remaining data later
CO2 <- window(co2, end = c(1990, 12))

## Set up the parameter sets over which we want to operate
CO2.pars <- expand.grid(ar = 0:2, diff = 1, ma = 0:2, sar = 0:1,
                        sdiff = 1, sma = 0:1)
## As you are only wanting ARMA, then you would need something like
## pars <- expand.grid(ar = 0:4, diff = 0, ma = 0:4)
## and where you choose the upper and lower limits - here 0 and 4

## A vector to hold the BIC values for each combination of model
CO2.bic <- rep(0, nrow(CO2.pars))

## loop over the combinations, fitting an ARIMA model and recording the BIC
## for that model. Note we use AIC() with extra penalty given by `k`
for (i in seq(along = CO2.bic)) {
    CO2.bic[i] <- AIC(arima(CO2, unlist(CO2.pars[i, 1:3]), 
                            unlist(CO2.pars[i, 4:6])),
                      k = log(length(CO2)))
}

## identify the model with lowest BIC
CO2.pars[which.min(CO2.bic), ]

## Refit the model with lowest BIC
CO2.mod <- arima(CO2, order = c(0, 1, 1), seasonal = c(0, 1, 1))
CO2.mod
## Diagnostics plots
tsdiag(CO2.mod, gof.lag = 36)

## predict for the most recent data
pred <- predict(CO2.mod, n.ahead = 7 * 12)
upr <- pred$pred + (2 * pred$se) ## upper and lower confidence intervals
lwr <- pred$pred - (2 * pred$se) ## approximate 95% pointwise

## plot what we have done
ylim <- range(co2, upr, lwr)
plot(co2, ylab = ylab, main = expression(bold(Mauna ~ Loa ~ CO[2])),
     xlab = "Year", ylim = ylim)
lines(pred$pred, col = "red")
lines(upr, col = "red", lty = 2)
lines(lwr, col = "red", lty = 2)
legend("topleft", legend = c("Observed", "Predicted", "95% CI"),
       col = c("black", "red", "red"), lty = c(1, 1, 2), bty = "n")

Answer 2

La forma más sencilla de llegar a los valores de $p$ y $q$ está utilizando auto.arima del paquete previsión . No hay una forma más sencilla en ningún paquete estadístico para llegar a buenos valores. La razón principal es que no existe una definición universal de bueno.

Ya que menciona el sobreajuste, una forma posible es ajustar modelos arima para diferentes valores de $p$ y $q$ y, a continuación, elija el que sea el mejor según sus criterios de sobreajuste (rendimiento de previsión fuera de la muestra, por ejemplo). auto.arima hace básicamente lo mismo, puede elegir entre AIC, AICC y BIC para que auto.arima elija el mejor modelo.