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
