Pronósticos de series temporales con reg-ARIMA y dummys temporales en R

En la siguiente entrada se simulan tres series autoregresivas (AR(1)) con coeficiente autoregresivo de 0.9 y una estructura de correlación dada (estimada a partir de Cholesky), para mostrar cómo realizar pronósticos ARIMA usando la función auto.arima de R, agregando regresores exógenos y realizando la transformación propuesta por Box-Cox. Se usan dos tipos de regresores: otras series de tiempo y dummy representando eventos o rangos de tiempo.

El primer bloque de códigos muestra como simular un conjunto de series temporales asumiendo una estructura de correlación entre estas [el código de esta simulación ha sido tomado de: stats.stackexchange.com].

library(tseries)

library(forecast)

n <- 100

corM <- rbind(c(1.0, 0.6, 0.8),

              c(0.6, 1.0, 0.5),

              c(0.8, 0.5, 1.0))

set.seed(100)

SigmaEV <- eigen(corM)

# eps <- rnorm(n * ncol(SigmaEV$vectors)) #Normal

eps  <- 10 + arima.sim(list(order=c(1,0,0), ar=.85), n=n * ncol(SigmaEV$vectors))

Meps <- matrix(eps, ncol = n, byrow = TRUE)   

Meps <- SigmaEV$vectors %*% diag(sqrt(SigmaEV$values)) %*% Meps

Meps <- ts(t(Meps), start = c(1991,1), frequency = 12)

cor(Meps)

autoplot(Meps)

Luego, se extrae las fechas del objeto ts (series de tiempo) que contiene nuestras series de tiempo o base de datos. Esto se realiza para poder agregar dummys temporales en función de fechas al modelo auto.arima identificado. Por tanto, luego de obtener el vector de fechas, podemos crear variables binarias a partir de operadores lógicos y relacionales, usando la función ifelse sobre el objeto fecha.arima.

fecha.xarima <- format(as.Date(Meps), "%Y-%m")

arimax_data <- cbind(Meps,

                     d1 = ifelse(fecha.xarima ==  "1991-06",1,0),

                     d2 = ifelse(fecha.xarima > "1998-03",1,0),

                     d3 = ifelse(fecha.xarima > "1995-01" & fecha.xarima < "1995-12",1,0))

head(arimax_data)

        Meps.Series 1 Meps.Series 2 Meps.Series 3 d1 d2 d3

Jan 1991     -5.380834     -16.21302     -8.985020  0  0  0

Feb 1991     -5.958240     -17.04476     -9.621183  0  0  0

Mar 1991     -5.914374     -16.31246     -9.252622  0  0  0

Apr 1991     -6.495822     -16.04300    -11.825374  0  0  0

May 1991     -5.014285     -14.42712     -9.175827  0  0  0

Jun 1991     -5.889311     -15.25092     -8.886918  1  0  0

Ahora, sobre la estructura de datos formada utilizamos el auto.arima de R para encontrar el modelo adecuado, tomando en consideración el ajuste de Box-Cox propuesto para estabilizar la varianza de los datos por medio de la función BoxCox.lambda que permite identificar el Lambda y de usar las variables endógenas. En este ejemplo pretendemos pronosticar la serie 1 y utilizar el resto de series como regresores del modelo ARIMA. La data se segmenta en dos periodos, uno de entrenamiento sobre el que se estima el modelo a partir de las primeras 95 observaciones y otro de evaluación del modelo sobre las 5 observaciones restantes. Dentro de la función auto.arima se incluye la variable a pronosticar trainData[,1] y el conjunto de regresores trainData[,-1] donde se incluyen las variables dummy temporales.

obsForcast <- tail(1:n, 5)

h <- length(obsForcast)

trainData <- arimax_data[-obsForcast,]

testData <- arimax_data[obsForcast,]

lamb.x <- BoxCox.lambda(trainData[,1])

lamb.x

[1] 1.999959

autoArimaReg <-auto.arima(trainData[,1],  xreg=trainData[,-1], seasonal=T, lambda = lamb.x, ic ="aic")

autoArimaReg

Series: trainData[, 1]

Regression with ARIMA(1,0,0) errors

Box Cox transformation: lambda= 1.999959

Coefficients:

         ar1  intercept  Meps.Series 2  Meps.Series 3       d1      d2      d3

      0.6397    25.8587         1.0771         2.4947  -3.5615  4.6658  0.7982

s.e.  0.0928     4.0908         0.3425         0.2760   1.9474  1.7414  1.6661

sigma^2 estimated as 5.575:  log likelihood=-213.04

AIC=442.09   AICc=443.76   BIC=462.52

Una vez estimado el modelo ARIMA es necesario verificar la no autocorrelación de su residuo como forma de validar los resultados del modelo y su posterior uso para generar pronósticos (verifique que en el ejemplo actuar claramente se rechaza la hipótesis nula de no autocorrelación en los residuos, por lo que, debe continuarse la búsqueda de un modelo para pronóstico):

# Evaluar la validez del modelo:

residuo <-autoArimaReg$residuals

Box.test(residuo, lag=20, type="Ljung-Box")

Box-Ljung test

data:  residuo

X-squared = 35.802, df = 20, p-value = 0.01623

Posteriormente se hace el pronóstico agregando los regresores exógenos para obtener los pronósticos. Es decir, deben otorgarse a la función los valores esperado de los regresores durante el horizonte de pronóstico. En base a lo anterior, note que se tenían los valores de los regresores en la ventana de prueba (horizonte de pronóstico), sin embargo, en la práctica por lo general no se cuenta con este dato, por lo que, se suelen realizar pronósticos ARIMA de las series individuales y usar reglas deterministas en el caso de las dummys temporales, dato que su valor indicador se conoce siempre.

fcast.cox <- forecast::forecast(autoArimaReg, xreg=testData[,-1], h=n.valid)

autoplot(fcast.cox)

La evaluación del pronóstico:

plot(arimax_data[obsForcast,1], type="l")

lines(as.vector(fcast.cox$mean))

accuracy(f=fcast.cox$mean, x=testData[,1])

                 ME      RMSE       MAE     MPE    MAPE

Test set -0.5159918 0.6478536 0.5159918 10.6098 10.6098