Construcción de modelos ARIMA estacionales (SARIMA) con R

#Construcción de modelos ARIMA estacionales multiplicativos:
#ARMA(P,Q)s estacional puro, causal e invertible. cuando contamos con diferenciación estacional.
#SARIMA o ARIMA(p,d,q)x(P,D,Q)s cuando contamos con diferenciación estacional y diferenciación ordinaria.

prod=ts(scan("c:/data/prod.dat"), start=1948, frequency=12) #analizamos los datos del BPI de la Reserva Federal
par(mfrow=c(2,1)); acf(prod, 48); pacf(prod, 48) # (P)ACF de los datos muestran un decaimiento lento del ACF y un pico en lag=1 del PACF, indicando comportamiento no.estacionario y la necesidad de diferenciar la serie (diferencia primera diff(Xt)).
par(mfrow=c(2,1)); acf(diff(prod), 48);pacf(diff(prod), 48) # (P)ACF de los datos diferenciados d1. vemos picos en lag=12, 24, 36, 48 con decaimiento lento, indicando la necesidad de diferenciación estacional de los datos (diff(diff(Xt))).
par(mfrow=c(2,1));acf(diff(diff(prod),12), 48);pacf(diff(diff(prod),12), 48) # (P)ACF de los datos d1-d12. Vemos un pico en lag=12 del ACF con menores picos para lag=24,36 y picos en lag=12, 24, 36, 48 del PACF. Por lo tanto, se sugieren los modelos MA estacional de orden Q=1 y el AR estacional de orden P=2 o ambos juntos. Asimismo, para lags no estacionales (h=1,...,11) ACF y PACf decaen, indicando que debemos ajustar modelos con p>0 y q>0 para el comportamiento no-estacional de la serie, en principio seleccionamos p=1 y q=1.
#los 3 modelos sugeridos son: ARIMA(1,1,1)x(0,1,1)12, ARIMA(1,1,1)x(2,1,0)12 y ARIMA(1,1,1)x(2,1,1)12

#ajustamos el tercer modelo sugerido ARIMA(1,1,1)x(2,1,1)12 y calculamos sus AIC
prod.fit3 = arima(prod, order=c(1,1,1),seasonal=list(order=c(2,1,1), period=12)); prod.fit3
tsdiag(prod.fit3, gof.lag=48) # diagnostics
prod.pr = predict(prod.fit3, n.ahead=12) ### forecasts for the final model
U = prod.pr$pred + 2*prod.pr$se; L = prod.pr$pred - 2*prod.pr$se
month=337:372; plot(month,prod[month],type="o",xlim=c(337,384),ylim=c(100,180),ylab="Production"); lines(prod.pr$pred, col="red",type="o");lines(U, col="blue", lty="dashed");lines(L, col="blue", lty="dashed"); abline(v=372.5,lty="dotted")
ts.plot(prod,prodpr$pred, col=1:2, type="o",ylim=c(105,175), xlim=c(1975,1980));lines(U, col="blue", lty="dashed");lines(L, col="blue", lty="dashed")

Referencias:

Comentarios