lunes, 10 de enero de 2011

ARIMA. Modelos no estacionarios estacionales


Muchas series temporales no son estacionarias, ya sea porque presentan tendencias o por efectos estacionales. En particular, las caminatas aleatorias, que caracterizan muchos tipos de series, no son estacionarias pero pueden transformarse en series estacionarias utilizando la diferenciación de primer orden. Como la serie diferenciada necesita ser agregada (o integrada) para recuperar la serie original, el proceso estocástico subyacente se llama media móvil integrada autorregersiva (ARIMA).
Los procesos ARIMA se pueden extender para incluir términos estacionales, dando paso a modelos ARIMA estacionales no-estacionarios (SARIMA).
Las series no-estacionarias porque la varianza está correlacionada en serie (heterocedástica condicional) resultan en períodos de volatilidad (un cambio claro en la varianza). Una aproximación para modelar este tipo de series es el modelo heterocedástico condicional autorregresivo (ARCH). Asimismo, podemos utilizar la generalización del modelo ARCH (GARCH).
  • Modelos ARIMA no-estacionales
Una serie temporal {X(t)} sigue un proceso ARIMA(p,d,q) si la diferencia dth de la serie es un proceso ARMA(p,q). Es decir, sea Y(t)=(1-B)^d*X(t) entonces theta_p(B)*Y(t)=phi_q(B)*W(t). Por tanto, podemos escribir el proceso ARIMA como_ thet_p(B)(1-B)^d*X(t)=phi_q(B)*W(t), donde theta_p y phi_q son polinomios de orden p y q, respectivamente..

Ejemplos:
  • ARIMA(0,1,1) o IMA(1,1): X(t)=X(t-1)+W(t)+beta*W(t-1), donde beta es el parámetro del modelo, o (1-B)*X(t)=(1+beta*B)*W(t).
  • ARIMA(1,1,0) o ARI(1,1): X(t)=alfa*X(t-1)+X(t-1)-alfa*X(t-2)+W(t) donde alfa es el parámetro modelado, o (1-alfa*B)(1-B)*X(t)=W(t).

Aplicación en R: ajuste ARIMA(1,1,1) para series simuladas
    #Ajuste del ARIMA(1,1,1) para series simuladas mediante el modelo xt = 0.5xt-1+xt-1-0.5xt-2+wt+0.3wt-1.
#manualmente
set.seed(1)
x <- w <- rnorm(1000); for (i in 3:1000) x[i] <- 0.5 * x[i - 1] + x[i - 1] - 0.5 *x[i - 2] + w[i] + 0.3 * w[i - 1]
arima(x, order = c(1, 1, 1))
#o automáticamente
x <- arima.sim(model = list(order = c(1, 1, 1), ar = 0.5,ma = 0.3), n = 1000)
arima(x, order = c(1, 1, 1))



Aplicación en R: ajuste ARIMA(0,1,1) para la producción de cerveza australiana
    #Ajuste del IMA(1,1) o ARIMA(0,1,1) para la producción de cerveza australiana
www <- "http://www.massey.ac.nz/~pscowper/ts/cbe.dat"
CBE <- read.table(www, he = T)
Beer.ts <- ts(CBE[, 2], start = 1958, freq = 12)
Beer.ima <- arima(Beer.ts, order = c(0, 1, 1)) #la serie tiene una tendencia creciente, por loque ajustamos un modelo IMA(1,1) (representa una tendencia linal con un ruido blanco)
Beer.ima #parámetros estimados xt=Xt-1+Wt-.33*wt-1.
acf(resid(Beer.ima)) #corelograma de los residuales. tiene picos en cada año y sugiere que se necesita un término estacional.
Beer.1991 <- predict(Beer.ima, n.ahead = 12) #predecimos los valores para los próximos 12 años
sum(Beer.1991$pred) #producción anual total para el año 1991



Aplicación en R: diferenciación de la serie de electricidad
#Ajuste de ARIMA no estacional a la serie de electricidad. Diferenciación para eliminar tendencia.
www <- "http://www.massey.ac.nz/~pscowper/ts/cbe.dat"
CBE <- read.table(www, he = T)
Elec.ts <- ts(CBE[, 3], start = 1958, freq = 12)
layout(c(1, 1, 2, 3))
plot(Elec.ts)
plot(diff(Elec.ts))
plot(diff(log(Elec.ts)) #eliminamos la tendencia creciente




  • Modelos ARIMA estacionales
El modelo ARIMA estacional usa la diferenciación en lag igual al número de estaciones (s) para elimina el efecto estacional aditivo. Al igual que la diferenciación a lag=1 que remueve la tendencia, la diferenciación a lag=s introduce un término de medias móviles, además del término autorregresivo en lag=s. El ARIMA estacional ARIMA(p,d,q)(P,D,Q)s puede expresarse usando el operador de retardo B:

Theta_P(B^s)*theta_p(B)(1-B^s)^D(1-B)^d*x(t)=Phi_Q(B^s)*phi_q(B)*W(t)

donde Theta_P, theta_p, Phi_Q y phi_q son polinomios de orden P, p, Q y q, respectivamente. El modelo es no-estacionario, pero si D=d=0 y las raíces de la ecuación característica son mayores a la unidad, el modelo resultante será estacionario.

Ejemplos:
  • AR(1) con período estacional de 12 unidades o ARIMA(0,0,0)(1,0,0)12: es el modelo xt = alfa*xt−12 + wt. Este modelo es apropiado para datos mensuales donde mes actual solo está influenciado por el valor del mes del año anterior.
  • ARIMA(0,1,0)(1,0,0)12 o xt = xt−1 + alfa*xt−12 – alfa*xt−13 + wt, es un modelo con tendencia estocástica e influencias estacionales. El cambio en el tiempo t depende del cambio en el mismo tiempo (i.e. mes) del año anterior.
  • MA(1) con período estacional de 4 unidades o ARIMA(0,0,0)(0,0,1)4: es el modelo xt=wt-beta*wt-4, estacionario y solo apto para datos sin tendencia.
  • ARIMA(0,1,0)(0,0,1)4 o xt=xt-1+wt-beta*wt-4: para datos con tendencia estocástica e influencias estacionales cuatrimestrales.
  • ARIMA(0,0,0)(0,1,1)4 o xt = xt−4 + wt – beta*wt−4: para datos donde el término estacional contiene una tendencia estocástica (la diferenciación se aplica al período estacional). con término el modelo se puede extender

Criterios de selección de modelos ARIMA:
Los modelos SARIMA pueden tener un gran número de parámetros y combinaciones de ellos. Por tanto, será apropiado considerar un amplio rango de modelos posibles y elegir entre ellos según un criterio apropiado (e.g. AIC). Para ello deberemos tener cuidado con el sobre-ajuste (el cual produce mejores AIC). Una vez que seleccionamos el mejor modelo a ajustar, realizamos el correlograma de los residuos para verificar que sigue un comportamiento de ruido blanco.

Aplicación en R: selección entre ajuste ARIMA(1,1,0)(1,0,0)12 y ARIMA(0,1,1)(0,0,1)12 para series de producción eléctrica, luego buscamos el mejor modelo según un procedimiento automático y lo analizamos.

    #Ajuste de ARIMA estacional a la serie de electricidad.
#ajustamos dos modelos ajuste ARIMA(1,1,0)(1,0,0)12 y ARIMA(0,1,1)(0,0,1)12 y evaluamos cuál se ajusta mejor a los datos mediante el criterio AIC
AIC (arima(log(Elec.ts), order = c(1,1,0),seas = list(order = c(1,0,0), 12)))
AIC (arima(log(Elec.ts), order = c(0,1,1),seas = list(order = c(0,0,1), 12)))
#el modelo ARI estacional provee el mejor ajuste, con el menor AIC
#para chequear varios modelos podemos utilizar la siguiente función con el criterio CAIC (que permite eliminar la sobre-parametrización)
get.best.arima <- function(x.ts, maxord = c(1,1,1,1,1,1)) #construimos una función para obtener el mejor modelo arima según el criterio CAIC
{
best.aic <- 1e8
n <- length(x.ts)
for (p in 0:maxord[1]) for(d in 0:maxord[2]) for(q in 0:maxord[3])
for (P in 0:maxord[4]) for(D in 0:maxord[5]) for(Q in 0:maxord[6])
{
fit <- arima(x.ts, order = c(p,d,q),
seas = list(order = c(P,D,Q),
frequency(x.ts)), method = "CSS")
fit.aic <- -2 * fit$loglik + (log(n) + 1) * length(fit$coef)
if (fit.aic < best.aic)
{
best.aic <- fit.aic
best.fit <- fit
best.model <- c(p,d,q,P,D,Q)
}
}
list(best.aic, best.fit, best.model)
}
best.arima.elec <- get.best.arima( log(Elec.ts),maxord = c(2,2,2,2,2,2)) #observamos que el mejor modelo ajustado usando términos de orden mayor a 2, es ARIMA(0,1,1)(2,0,2)12.
best.fit.elec <- best.arima.elec[[2]]
acf( resid(best.fit.elec) ) #los residuales parecen aproximarse a un ruido blanco
best.arima.elec [[3]]
ts.plot( cbind( window(Elec.ts,start = 1981),exp(predict(best.fit.elec,12)$pred) ), lty = 1:2) #se utiliza un factor de correción sesgado para los valores predichos. Sin embargo, esa corrección parece innecesaria en este caso porque la desviación estandar de los residuos es pequeña en comparación con las predicciones.



RESUMEN DE PROCESOS ARIMA ESTACIONALES (SARIMA)

  • Procesos no estacionarios, donde analizaremos la falta de estacionariedad en la media (comportamiento estacional).
  • Generalmente podemos incorporar la estacionalidad dentro de un modelo ARIMA de una forma multiplicativa, resultando en un modelo ARIMA estacional multiplicativo.
  • Concepto de estacionalidad y sus tipos:
  1. Una serie estacional presenta valores que no son constantes pero varían con una pauta cíclica. Cuando E[Xt]=E[Xt+s] decimos que s es el período de estacionalidad, y éste define el número de observaciones que forman el ciclo estacional: s=12 (serie mensual), s=4 (trimestral), s=7 (semanal).
  2. El modelo más simple trata un efecto constante que se suma a los valores de la serie: Xt=St+nt. La serie se escribe como suma de un componente estacional St y un proceso estacionario nt (con mu su media).
  3. El componente estacional St puede tener comportamiento: 1) determinista (función constante para el mismo mes en distintos años) St=St+ks, 2) estacionario (evoluciona en el tiempo y su evolución es estacionaria, oscilando alrededor de un valor medio) St=mu+vt, donde vt es un proceso estacionario de media cero que introduce variabilidad en el año, o 3= no-estacionario (es cambiante sin ningún valor medio fijo) St=St-s+vt.
En los 3 casos podemos convertir una serie estacional en estacionaria aplicando una diferencia estacional. Definimos el operador diferencia estacional de período s: diffs=1-B^s, es decir diffsXt=(1-B^s)Xt=Xt-Xt-s p diffsXt=diffsSt+diffsnt para nuestro modelo más simple de estacionalidad.
En el modelo ARIMA estacional podemos convertir series no estacionarias en estacionarias tomando diferencias regulares (entre períodos consecutivos) o eliminar la estacionalidad mediante diferencias estacionales, ambos mediante la transformación: Wt=diffsDdiffdXt, donde D es el número de diferencias estacionales (generalmente D=1) y del número de diferencias regulares (d<=3). Cuando existe dependencia estacional podemos generalizar el modelo ARMA para series estacionarias incorporando además la dependencia regular (asociada a intervalos de medida de la serie) y la dependencia estacional (asociada a observaciones separadas por s períodos). Debemos modelar 2 tipos de dependencias:
  1. Incorporar la dependencia estacional a la regular, añadiendo a los operadores AR o MA el operador B, términos B^s para representar la dependencia entre observaciones separadas por s períodos.
  2. Modelar de forma separada la dependencia regular y la estacional, y construir el modelo incorporando ambas de forma multiplicativa: modelo ARIMA estacional multiplicativo. Esta opción es la más simple.
Esta clase de modelos se escriben de forma simplificada como el modelo ARIMA(P,D,Q)sx(p,d,q).
Nota: el modelo ARIMA estacional multiplicativo se basa en la hipótesis central de que la relación de dependencia estacional (modelo estacional) es la misma para todos los períodos.

En el ACF:
fijarse en los retardos iniciales 1,2,3,4,… para identificar la estructura regular, y en los retardos estacionales s, 2s, 3s,… para identificar la estructura estacional. La interacción alrededor de los coeficientes estacionales puede entonces utilizarse como confirmación de la identificación realizada.
En el PACF: fijarse en los retardos iniciales 1,2,3,4… para identificar la estructura regular y en los retardos estacionales s,Ds,3s,… para la identificación de la estructura estacional.


IDENTIFICACIÓN DE LOS POSIBLES MODELOS ARIMA

  1. Selección del modelo de la serie estacionaria: identificar la estructura no-estacionaria (si existe) y después la estructura ARMA estacionaria.
  • a. Decidir qué transformaciones aplicar para convertir la serie observada en una serie estacionaria. Transformar la serie para que tenga varianza constante y/o diferenciar la serie para que tenga media constante.
  • -Si la impresión visual de la serie es que la variabilidad no aumenta con el tiempo sino con el nivel de la serie, nos fijaremos en el gráfico de desviación típica (medida de la variabilidad) vs. Media local (medida del nivel) para confirmar o no nuestra sospecha visual. Cuando la variabilidad aumenta linealmente con el nivel de la serie, tomaremos logaritmos para que la variabilidad sea constante. En general, aplicamos la transformación Box-Cox Yt={xt^(1-alfa)-1}/[1-alfa], donde alfa lo estimamos como la pendiente del gráfico log(si) vs. Log(mean(xi)), donde si es la desviación típica y mean(xi) es la media. Excepciones: si la varianza cambia sin relación con el nivel o si es constante.
  • b. Determinar un modelo ARMA para la serie estacionaria (p y q de su estructura AR y MA) y si el proceso es estacional, determinar los órdenes P, Q de la estructura ARMA estacional.
  • Estimaremos todos los modelos que consideremos posibles para explicar la serie y después elegir entre ellos.

  1. Estimación: los parámetros AR y MA del modelo se estiman por máxima verosimilitud y se obtienen sus errores estándar y los residuos del modelo.
  2. Diagnosis: se comprueba que los residuos no tienen estructura de dependencia y siguen un proceso de ruido blanco. Si los residuos tienen estructura modificaremos el modelo para incorporarla y repetiremos las 3 etapas anteriores hasta obtener un modelo adecuado.

Referencias
  • Cowpertwait, Paul S.P. & Metcalfe, Andrew V. (2009). Introductory Time Series with R. Springer-Verlag.

No hay comentarios:

Publicar un comentario en la entrada

Libros para descargar sobre Análisis de Series Temporales