ARMA. Modelos estacionarios no estacionales

Los modelos de regresión pueden dar cuenta de componentes no estacionarios, sin embargo, no consideran la correlación temporal de los residuales. Por ello, consideraremos modelos estacionarios que pueden ser utilizados para series de residuos que contengan tendencias o ciclos estacionales no obvios. El modelo estacionario ajustado puede combinarse con un modelo de regresión para mejorar los pronósticos.


1. Modelos de media móvil MA(q)

Definición: un proceso MA de orden 1 es una combinación lineal del término de ruido blanco actual y los q términos de ruido blanco del pasado más reciente.

X(t) = w(t) +beta_1*w(t−1) + . . . + beta_q*w(t−q)

donde {w(t)} es un ruido blanco con media cero y varianza sigma^2.
También podemos escribir esta ecuación en términos del operador de cambio de retardo B:

X(t)=(1+beta_1*B+…+beta_q*B^q)*w(t)=phi_q(B)*w(t)

donde phi_q es un polinomio de orden q.
Como el proceso MA consiste en la suma de términos de ruido blanco estacionarios, son estacionarios y tienen media y autocovarianza invariantes en el tiempo.

Función “arima”: El modelo MA(q) puede ser ajustado a los datos utilizando la función “arima” seleccionando el orden de los parámetros como c(0,0,q). Por defecto, la función “arima” no resta la media y estima un término de intercepto. Asimismo, minimiza la suma condicional de cuadrados para estimar los valores de los parámetros (para verlos debemos seleccionar “method”=c(“CSS”) o utilizarlos como entradas para estimación por máxima verosimilitud). Podemos seleccionar el valor de la media como cero, en logar de estimar el intercepto, utilizando “incluse.mean=FALSE” dentro de la función “arima”.

Aplicación en R: ajuste de modelos MA a series simuladas según el modelo MA(3) 
set.seed(1)

b <- c(0.8, 0.6, 0.4)

x <- w <- rnorm(1000)

for (t in 4:1000) {

for (j in 1:3)

{

x[t] <- x[t] + b[j] * w[t - j]

}

plot(x, type = "l")

acf(x)

x.ma <- arima(x, order = c(0, 0, 3))

x.ma #observamos los parámetros estimados y los IC al 95%, que contienen los valores de los parámetros subyacentes.



Aplicación en R: ajuste del modelo MA(1) para la serie de tasas de intercambio 
#Ajuste de modelo MA(1). Serie de tasa de intercambio.

www <- "http://www.massey.ac.nz/~pscowper/ts/pounds_nz.dat"

x <- read.table(www, header = T)

x.ts <- ts(x, st = 1991, fr = 4)

x.ma <- arima(x.ts, order = c(0, 0, 1)) #ajuste MA(1)

x.ma

acf(x.ma$res[-1]) #el modelo no provee un ajuste satisfactorio, los residuos no son una realización realista del ruido blanco



2. Modelos mixtos ARMA

Definición: cuando juntamos los términos AR y MA en una única expresión obtenemos un modelo de media móvil autorregresivo (ARMA) de orden (p,q):

X(t)=alfa_1*X(t-1)+…+ alfa_p*X(t-p)+W(t)+beta_1*W(t-1)+…+beta_q*W(t-q)

donde {W(t)} es un ruido blanco. En términos del operador de retardos B, tenemos:

theta_p(B)*X(t)=phi_q(B)*W(t)

El modelo ARMA(p,q) tiene las siguientes características:
• Es estacionario cuando las raíces de theta son mayores a la unidad, en valor absoluto.
• Es invertible cuando las raíces de phi son mayores a la unidad, en valor absoluto.
• El modelo AR(p) es un caso especial del ARMA(p,0)
• El modelo MA(q) es un caso especial del ARMA(0,q)
• Parsimonia sobre los parámetros: el modelo ARMA es más eficiente (i.e. requiere menores parámetros) que los modelos AR o MA solos.
• Redundancia de los parámetros: cuando theta y phi tienen un factor común, el modelo estacionario puede simplificarse.

Función “arima.sim”: los procesos ARMA y ARIMA pueden simularse con esta función.

Aplicación en R: ajuste una series simulada según el modelo ARMA(1,1) con alfa=-.6 y beta=.5. 
   #Ajuste de modelo ARMA(1,1) para una serie simulada con el modelo ARMA(1,1) y parámetros alfa=-.6 t beta=.5.

set.seed(1)

x <- arima.sim(n = 10000, list(ar = -0.6, ma = 0.5))

coef(arima(x, order = c(1, 0, 1))) #como era de esperar las estimaciones muestrales para alfa y beta son cercanos a los parámetros del modelo subyacente


Aplicación en R: ajuste de modelos ARMA para las series de tasas de intercambio. 
#Ajuste de modelos ARMA para las series de tasa de intercambio.

x.ma <- arima(x.ts, order = c(0, 0, 1)) #ajuste MA(1) o ARIMA(0,0,1)

x.ar <- arima(x.ts, order = c(1, 0, 0)) #ajuste AR(1) o ARIMA(1,0,0)

x.arma <- arima(x.ts, order = c(1, 0, 1)) #ajuste ARMA(1,1) o ARIMA(1,0,1)

AIC(x.ma);   AIC(x.ar) ;   AIC(x.arma)  #comparación de los modelos utilizando el criterio AIC. El ARMA(1,1) da mejor ajuste a los datos, seguido por el AR(1) y finalmente el MA(1).

x.arma

acf(resid(x.arma)) #los residuales del ARMA(1,1) tienen pequeñas autocorrelaciones, lo cual es consistente con la realización de un ruido blanco.


Aplicación en R: Ajuste de modelos ARMA para las series de producción eléctrica. 
#Ajuste de modelo ARMA para la producción de electricidad.

www <- "http://www.massey.ac.nz/~pscowper/ts/cbe.dat"

CBE <- read.table(www, header = T)

Elec.ts <- ts(CBE[, 3], start = 1958, freq = 12)

tsdisplay(Elec.ts) #existe una tendencia positiva y un ciclo estacional regular. La varianza aumenta con el tiempo sugiriendo una transformación log.

Time <- 1:length(Elec.ts)

Imth <- cycle(Elec.ts)

Elec.lm <- lm(log(Elec.ts) ~ Time + I(Time^2) + factor(Imth)) #se ajuste un modelo de regersión a los logaritmos de la serie original.

acf(resid(Elec.lm))  #el correlograma de los residuales parece ciclar con período de 12 meses sugiriendo que las variables indicadoras mensuales no son suficientes para dar cuenta de la estacionalidad en la serie.

#podemos dar cuenta de esto mediante un modelo no-estacionario con un componente estacional estocástico

best.order <- c(0, 0, 0)

best.aic <- Inf

for (i in 0:2)

{

for (j in 0:2)

 {

 fit.aic <- AIC(arima(resid(Elec.lm), order = c(i, 0,j)))

 if (fit.aic < best.aic)

   {

   best.order <- c(i, 0, j)

   best.arma <- arima(resid(Elec.lm), order = best.order)

   best.aic <- fit.aic

   }

 }

}

best.order

acf(resid(best.arma)) #correlograma del modelo ARMA de los residuales con mayor AIC



#predecimos la producción eléctrica para cada mes de los siguientes 3 años

new.time <- seq(length(Elec.ts), length = 36)

new.data <- data.frame(Time = new.time, Imth = rep(1:12,3))

predict.lm <- predict(Elec.lm, new.data)

predict.arma <- predict(best.arma, n.ahead = 36)

elec.pred <- ts(exp(predict.lm + predict.arma$pred), start = 1991,freq = 12)

ts.plot(cbind(Elec.ts, elec.pred), lty = 1:2)


Aplicación en R: ajuste de modelos ARMA para los datos de ondas en experimentos de tanques/baldes. 
 #Ajuste ARMA. Datos del tanque de ondas.

www <- "http://www.massey.ac.nz/~pscowper/ts/wave.dat"

wave.dat <- read.table(www, header = T)

attach (wave.dat)

layout(1:3)

plot (as.ts(waveht), ylab = 'Wave height')

acf (waveht);   pacf (waveht) #sugiere que p debe ser al menos 2

wave.arma <- arima(waveht, order = c(4,0,4))  #el mejor modelo aRMA(p,q) obtenido según la varianza mínima de los residuales, se obtiene para p=q=4. Ajustamos un ARMA(4,4)

acf (wave.arma$res[-(1:4)]);   pacf (wave.arma$res[-(1:4)]) #se adecua a un proceso de ruido blanco

hist(wave.arma$res[-(1:4)], xlab='height / mm', main='')



RESUMEN DE PROCESOS ARMA
  • Son combinación de las propiedades de los procesos AR y MA y permiten representar de forma “escueta” (utilizan pocos parámetros) procesos cuyas primeros q coeficientes son cualesquiera, mientras que los siguientes decrecerán según leyes simples.
ARMA(1,1)
Xt=phi*Xt-1+epsilont-phi1*epsilont-1 o (1-phi1B)Xt=(1-theta1B)epsilon, donde |phi1|<1>
ARMA(p,q)
(1-phi1B-…-phipB^p)Xt=(1-theta1B-…-thetaqB^q)epsilont o phip(B)Xt=thetaq(B)epsilont.



Referencias
  • Cowpertwait, Paul S.P. & Metcalfe, Andrew V. (2009). Introductory Time Series with R. Springer-Verlag.
  • Cowpertwait, Paul S.P. & Metcalfe, Andrew V. (2009). http://elena.aut.ac.nz/~pcowpert/ts/

Comentarios

Publicar un comentario