Modelos de Series de Tiempo
Modelos Estacionarios
Las series de tiempo frecuentemente tendrán componentes bien definidos, como la tendencia o el patrón estacional. Una regresión lineal bien elegida puede servir para estos componentes no estacionarios. Sin embargo, los residuos están usualmente correlacionados en el tiempo, y esto no se toma en cuenta a la hora de hacer la regresión. En esta sesión se verán modelos estacionarios que pueden ser adecuados para series residuales que no contienen tendencias o ciclos estacionales. Estos modelos estacionarios pueden ser combinados con modelos de regresión para obtener pronósticos. Los modelos autorregresivos introducidos en sesiones pasadas ofrecen muchas veces modelos satisfactorios para las series de tiempo residuales.
Series Estrictamente Estacionarias
Una serie de tiempo \(x_t\) es estrictamente estacionaria si la distribución estadística conjunta de \(x_{t_1}, ..., x_{t_m}\) es la misma que la distribución conjunta de \(x_{t1+m},...,x_{tn+m}\) para todo \(t_1,...,t_n\) y \(m\), por lo tanto, la distribución no varía después de un cambio arbitrario de tiempo. Se debe notar que la estacionariedad estricta implica que la media y la varianza son constantes en el tiempo y la autocovarianza \(Cov(x_t,x_s)\) solo depende del rezago \(k = |t-s|\) y puede ser escrita como \(y(K)\). Si una serie no estrictamente estacionaria pero la media y la varianza son constantes en el tiempo y la autocovarianza solo depende del rezago, la serie tenrdrá estacionariedad de segundo orden. La estacionariedad es una idealización que propiedad en los modelos. Si se ajusta un modelo estacionario a los datos, se asume que nuestra data es una realización de un proceso estacionario. Es por eso que el primer paso para realizar este tipo de ajuste es observar si existe alguna evidencia de tendencia o afecto estacional, y si es el caso, removerla.
Modelos de Mevias Móviles
Proceso MA(q)
Un proceso de medias móviles (MA) de orden \(q\) es una combinación lineal del término de ruido blanco y de los más recientes \(q\) términos de ruido blanco pasados, y está definida así:
\[ x_t = w_t+\beta_1w_{t-1}+...+\beta_q w_{t-q} \] Donde \(w_t\) es un proceso de ruido blanco con media 0 y varianza \(\sigma_w^2\). La ecuación anterior puede ser reescrita en términos del operador de rezago.
\[ x_t = (1+\beta_1B+...+\beta_qb^q)w_t = \phi_q(B)w_t \] Donde \(\phi_q\) es un polinomio de orden \(q\). Ya que los procesos MA consisten en sumas finitas de términos de ruido blanco estacionarios, estos son estacionarios y por lo tanto tienen una media y una autocovarianza que no varía en el tiempo.
Simulación en R
La simulación de un proceso MA(q) puede ser implementada en R, el siguiente código muestra la simulación de un proceso de Medias Móviles en R con parámetros \(\beta_1=0.8, \beta_2=0.6 \text{ y } \beta_3=0.4\).
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")
Las dos últimas líneas del código muestran os gráficos del proceso MA(3) y su correlograma. Se puede observar que este proceso tiene en sus primeras 3 autocorrelaciones significancia diferente de 0.
Ajuste de un proceso de Medias Móviles
Un modelo MA(q) puede ser ajustado en R usando la función arima con la función order fijada a c(0,0,q). De esta forma diferente a la función ar , la función arima minimiza la suma condicional de cuadrados de los valores estimados de los parámetros y devolverá estos si method=c(“CSS”) es especificado o usados como valores iniciales para la estimación por máxima verosimilitud. En el siguiente código se ajusta un modelo de medias móviles para la serie simulada anteriormente.
##
## Call:
## arima(x = x, order = c(0, 0, 3))
##
## Coefficients:
## ma1 ma2 ma3 intercept
## 0.7898 0.5665 0.3959 -0.0322
## s.e. 0.0307 0.0351 0.0320 0.0898
##
## sigma^2 estimated as 1.068: log likelihood = -1452.41, aic = 2914.83Tal y como se esperaba, la serie estimada tiene valores estimados en parámetros casi idénticos a los valores simulados, todos significativos. Además, como se esperaba el intercepto no es singificativamente diferentes de su valor de parámetros subyacente de cero.
El Proceso Arma
Recapitulando los conceptos anteriores, un proceso autoregresivo de orden \(p\), es un AR(p) si: \[ x_t = \alpha_1x_{t-1}+\alpha_2x_{t-2}+...+\alpha_px_{t-p}+w_t \] Donde \(w_t\) es un ruido blanco y los \(\alpha_i\) son parámetros del modelo. Una clase muy útil de modelo son obtenidos cuando los términos AR y MA son puestos juntos en la misma expresión. Una serie de tiempo \(x_t\) sigue un proceso autoregresivo de medias móviles (ARMA) de orden (p,q), denotado como ARMA(p,q), cuando:
\[ x_t = \alpha_1 x_{t-1} + \alpha_2 x_{t-p} + w_t + \beta_1 w_{t-1} + \beta_2 w_{t-2} + ... + \beta_q w_{t-q} \] Donde \(w_t\) es el ruido blanco. Esta ecuación puede ser representada en términos del operador de rezago y reordenada en la forma polinomial más concisa.
\[ \phi(B)x_t = \phi_q(B)w_t \] Los siguientes puntos acerca del proceso ARMA(p,q) deben ser considerados:
a.- El proceso estacionario si las raíces de \(\theta\) exceden a la unidad en un valor absoluto.
b.- El proceso es invertible cuando las raíces de \(\phi\) exceden todas a la unidad en un valor absoluto.
c.- El modelo AR(p) es un caso especial de un proceso ARMA(p,0).
d.- El modelo MA(q) es un caso especial de un modelo ARMA(0,q).
e.- Parsimonia paramétrica. Cuando se ajustan los datos, un modelo ARMA será a menudo más eficiente de los parámetros que un modelo AR o MA.
Simulación y ajuste en R
El proceso ARMA, puede ser simulado usando a función arima.sin, que toma una lista de coeficientes representando a los parámetros AR y MA. Un modelo AMRA(p,q) puede ser ajustado haciendo uso de la función arima con la función order fijada en c(p,0,q). El algoritmo de ajuste procede de manera similar al proceso MA. El código debajo simula un proceso ARMA(1,1) con parámetros \(\alpha = 0.6\) y \(\beta = 0.5\) y estima los parámetros del modelo para la serie simulada.
La estimación vendrá dada por:
## ar1 ma1 intercept
## -0.596966371 0.502703368 -0.006571345Como era de esperarse, las estimaciones de los parámetros \(\alpha\) y \(\beta\) tienen valores muy cercanos a los parámetros subyacentes del modelo.
Pronósticos
La función predict puede ser usada para predecir valores de una regresión, esta función estima también los errores estándar asociados con el modelo de regresión. Para los modelos ARMA usando la función arma. La función predict solamente requiere el número de períodos adelante para el pronóstico deseado.
Una vez que se importan los datos, se crea la serie de tiempo y se estima el modelo usando a función arima. Luego de esto se hace uso de la función predict para hacer las predicciones a partir del modelo estimado.
gnp.ts <- ts(gnp96$gnp96, start = 1967, freq = 4)
gnp.arma <- arima(x = gnp.ts, order = c(1,0,1), method = "CSS")
predict.arma <- predict(gnp.arma, n.ahead=12)Luego se crea una serie de tiempo que contiene a la predicción, especificando el valor e inicio y a frecuencia de la misma, que este caso es trimestral.
Los valores de la serie pueden ser visualizados de la siguiente manera:
## Qtr1 Qtr2 Qtr3 Qtr4
## 2002 9737.662 9811.786
## 2003 9886.538 9961.924 10037.948 10114.615
## 2004 10191.933 10269.905 10348.538 10427.836
## 2005 10507.807 10588.454Los errores estándar pueden también ser observados, y se notará en estos que a medida que avanza el tiempo pronósticado el error es mayor.
## Qtr1 Qtr2 Qtr3 Qtr4
## 2002 44.37924 70.82840
## 2003 90.08688 106.14778 120.30243 133.16353
## 2004 145.07873 156.26756 166.87906 177.01986
## 2005 186.76929 196.18829Finalmente, para observar gráficamente el pronóstico y los valores observados hacemos uso de la siguiente sintáxis:
Modelos No Estacionarios
Ya se ha visto que existen diversas series de tiempo que no son estacionarias a causa de las tendencias o los efectos estacionales que presentan. En particular, los caminos aleatorios, que caracterizan a muchos tipos de series, son no estacionarios, pero pueden ser transformadas a series estacionarias haciendo diferenciación de primer orden para las mismas. En esta sesión se extenderán los modelos ARMA. Como las series diferenciadas necesitan agregarse (o integrarse) para recuperar a la serie original, el proceso estocástico subyacente es llamado autorregresivo integrado de medias móviles, que abreviadamente es ARIMA.
Modelos Arima No Estacionales
Diferenciación de Series
Diferenciando una serie, se pueden remover tendencias, ya sean estas tendencias estocásticas, como un camino aleatorio, o determinísticas como en el caso de una tendencia lineal. El caso de un proceso random walk, \(x_t = x_{t-1}+w_t\), la serie diferenciada de primer orden es ruido blanco {wt}, es decir, \(\triangle x_t = x_t - x_{t-1} = w_t\), por lo tanto, es estacionario. La función arima en R no permite que los modelos diferenciados que se ajustan incluyan una constante. Si se quiere ajustar un modelo diferenciado a una tendencia determinística en R, se necesita hacer una diferenciación, la media ajusta la serie diferenciada para que la media sea cero, y entonces se estima un modelo ARMA a la serie diferenciada y ajustada usando la función arima con la opción include.mean fijada en FALSE y con d=0. El código siguiente muestra una serie con tendencia en donde se obtendrá la diferenciación y se realizará un gráfico de comparación entre la serie diferenciada y la observada.
Y en la serie diferenciada, se puede observar que para el gráfico diferenciado, que ya no puede ubicar una tendencia ni ningún patrón estacional en dicha serie.
El Modelo Integrado
Una serie {\(x_t\)} es integrada de orden \(d\), denorada como \(I(d)\), si la diferencia de orden \(d\) de la serie es un ruido blanco, es decir \(\triangle^d x_t = w_t\). Donde \(\triangle^d = (1-B)^d\), donde B es el operador de rezago, una serie será integrada de orden \(d\) si:
\[ (1-B)^d x_t = w_t \] El proceso de camino aleatorio es un caso especial de \(I(1)\). El comando diff puede ser usado para obtener diferenciaciones de mayor orden, ya sea por la aplicación repetida de este o el ajuste de \(d\) en los valores requeridos. Por ejemplo, diff(diff(x)) y diff(x, d=2) producirán ambos series diferenciadas de segundo orden para x. La diferenciación de segundo orden puede ser útil a veces para suavizar una serie con una curva de tendencia subyacente en el ruido blanco.
El parámetro lag puede ser usado para fijar el rezago de la diferenciación. Por defecto, lag está fijado en la unidad, pero otros valores pueden ser útiles para remover los efectos aditivos estacionales. Por ejemplo, dif(x, lag=2) removerá la línea de tendencia y los efectos aditivos estacionales en una serie con frecuencia mensual.
Definición y Ejemplos
Una serie de tiempo {\(x_t\)} sigue un proceso ARIMA(p,d,q) si las diferencias de orden d para la serie \(d\) siguen un proceso ARMA(p,q). Si se introduce \(y_t = (1-B)^d x_t\), entonces \(\phi_p (B)y_t = \phi_q(B)w_t\). Podemos sustituir para \(y_t\) y se obtiene la forma más breve del proceso ARIMA(p,d,q):
\[ \phi_p(B)(1-B)^d x_t = \phi_q(B) w_t \] Donde \(\phi_p\) y \(\phi_q\) son polinomios de orden \(p\) y \(q\), respectivamente. Algunos ejemplos de modelos ARIMA son:
- \(x_t = x_{t-1} + w_t + \beta w_{t-1}\), donde \(\beta\) es un parámetro del modelo. Este modelo es un ARIMA(0,1,1), que es llamado a veces modelo integrado de medias móviles, denotado como IMA(1,1). En general, ARIMA(0,d,q) = IMA(d,q).
- \(x_t = \alpha x_{t-1} + x_{t-1} + \alpha x_{t-2} + w_t\), donde \(\alpha\) es un parámetro del modelo. Este modelo es un ARIMA(1,1,0), que es llamado a veces modelo autoregresivo integrado, denotado como ARI(1,1). En general, ARIMA(p,d,0) = ARI(p,d).
Simulación y ajuste
Un proceso ARIMA(p,d,q) puede ser ajustado haciendo uso de la función arima con el parámetro order fijado en c(p,d,q). Estos procesos también pueden ser simulados en R escribiendo el código apropiado. Por ejemplo, en el código debajo, data para el modelo \(ARIMA(1,1,1) x_t = 0.5x_{t-1} + x_{t-1} - 0.5x_{t-2} + w_t + 0.3w_{t-1}\) es simulada y posteriormente ajustada.
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))##
## Call:
## arima(x = x, order = c(1, 1, 1))
##
## Coefficients:
## ar1 ma1
## 0.4235 0.3308
## s.e. 0.0433 0.0450
##
## sigma^2 estimated as 1.067: log likelihood = -1450.13, aic = 2906.26La serie se puede graficar de la siguiente forma:
El escribir un código propio tiene la ventaja de que se puede comprender mejor el modelo estimado. Sin embargo, R ofrece una función que permite simular procesos ARIMA, esta función es arima.sim, donde los parámetros modelo y n especifican el modelo y la longitud de la simulación, respectivamente:
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))##
## Call:
## arima(x = x, order = c(1, 1, 1))
##
## Coefficients:
## ar1 ma1
## 0.5567 0.2502
## s.e. 0.0372 0.0437
##
## sigma^2 estimated as 1.079: log likelihood = -1457.45, aic = 2920.91Modelos Arima Estacionales
Un modelo ARIMA estacional usa la diferenciación en el rezago igual al número de períodos \(s\) para remover efectos aditivos estacionales. El modelo ARIMA estacional incluye términos autoregresivos y de medias móviles en el rezago \(s\). El modelo \(ARIMA(p,d,q)(P,D,Q)_S\) puede resumirse en la siguiente expresión usando el operador rezago:
\[ \phi_p(B)^s\phi_p(1-B^s)^D(1-B)^d x_t = \phi_Q(B^s)\phi_q(B)w_t \] Donde \(\phi_P, \phi_p, \phi_Q\) y \(\phi_q\) son polinomios de orden \(P, p, Q\) y \(q\), respectivamente. En general, el modelo es no estacionario, aunque si D = d = 0 y las raíces de los polinomios de la ecuación exceden todas a la unidad en términos absolutos, el modelo resultante será estacionario. Como ejemplo, consideremos una serie de tiempo trimestral que es la suma de una tendencia lineal, 4 términos aditivos y un término de ruido blanco:
\[ x_t = a + b_t + s_{[t]} + w_t \] Donde [t] es el resto después de dividir \(t\) entre 4, entonces \(s_{[t]} = S_{[t-4]}\). Primero, consideremos la diferenciación de primer orden en el rezago 4 solamente:
\[ (1-B^4) x_t = x_t - x_{t-4} \\ = a + bt - (a+b(t-4)) + s_{[t]} - s_{[t-4]} + w_t - w_{t-4} \\ = 4b + w_t - w_{t-4} \] Formalmente, el modelo puede ser expresado como un \(ARIMA(0,0,0)(0,1,1)_4\) con un término constante de \(4b\). Ahora suponiendo que se aplica una diferenciación de primer orden en el rezago 1 antes de diferenciar el rezago 4, entonces:
\[ (1-B^4)(1-B)x_t = (1-B^4)(b + s_{[t]}- s_{[t-1]} + w_t + w_{t-1}) \\ = w_t - w_{t-1} + w_{t-4} + w_{t-5} \] Que es un \(ARIMA(0,1,1)(0,1,1,)_4\) sin término constante.
Estimación en R
Los modelos 𝐴𝑅𝐼𝑀𝐴 estacionales pueden tener un número largo de parámetros y combinaciones de términos. Por lo tanto, es apropiado que se prueben cierto número de modelos cuando se ajustan los datos y se deberá elegir el modelo de mejor ajuste haciendo uso del criterio AIC. En R, para ajusta modelos 𝐴𝑅𝐼𝑀𝐴 se hace uso de la función arima, especificando la opción seasonal para indicar el componente estacional. Para hacer la comparación de los modelos podemos estimar el AIC y hacer la comparación de los modelos.
## [1] -946.6604## [1] -930.0554Esta Sintaxis permite comparar dos modelos en base a su AIC, el primer modelo que tiene la forma \(ARIMA(1,1,0)(1,0,0)_4\) es el que tiene una mejor capacidad explicativa.
Modelos ARCH
Estos modelos sirven para hacer modelamiento de la volatilidad, ya que se requiere un modelo que permita cambios condicionales en la varianza. Una aproximación a esto es hacer uso de un modelo autoregresivo para la varianza. Esto lleva a la siguiente definición: una serie {\(e_t\), es autoregresiva de heterocedasticidad de primer orden, denotada como \(ARCH(1)\), si:
\[ e_t = w_t\sqrt{\alpha_0 + \alpha_1e_{t-1}^2} \]
Para ver como se introduce la volatilidad, se usa la ecuación anterior para calcular la varianza:
\[ Var(e_t) = E(e_t^2) \\ = E(w_t^2)E(\alpha_0+\alpha_1 e_{t-1}^2) \\ = E(\alpha_0 + \alpha_1 e_{t-1}^2) \\ = \alpha_0 + \alpha_1Var(e_{t-1}) \] Si se compara esta ecuación con el proceso \(AR(1)x_t=\alpha_0+\alpha_1x_{t-1}+w_t\) se podrá observar que la varianza de un proceso ARCH(1) se comporta de la misma manera que un proceso AR(1).
Extensiones y modelos GARCH
Un modelo ARCH de primer orden puede ser extendido hasta el p-ésimo orden incluyendo rezagos más altos. Un proceso ARCH(p) está dado por:
\[ e_t = w_t\sqrt{\alpha_0 + \sum_{t=i}^p \alpha_p e_{t-i}^2} \] Donde \(w_t\) es de nuevo un ruido blanco con media 0 y varianza unitaria. Una extensión adicional, muy usado en aplicaciones financieras es el modelo ARCH generalizado denorado como GARCH(q,p), que tiene al modelo ARCH(p) como un caso especial GARCH(0,p). Una serie sigue un proceso GARCH(q,p)si:
\[ e_t = w_t \sqrt{h_t} \] Donde:
\[ h_t = \alpha_0 + \sum_{i=1}^p \alpha_i e_{t-i}^2 + \sum_{j=1}^q \beta_j h_{t-j} \] Y \(\alpha_i\) y \(\beta_j(i=0,1,...,p; j=1,...,q)\) son parámetros del modelo. En R un modelo GARCH puede ser ajustado haciendo uso de la función GARCH del paquete tseries.
Simulación y ajuste del modelo
El siguiente código simula n GARCH(1,1) de forma \(a_t=w_t\sqrt{h_t}\), donde \(h_t=\alpha_0 +\alpha_1 a_{t-1} + \beta_1 h_{t-1}\) con \(\alpha_1 + \beta_1 < 1\) para asegurar estabilidad. La serie simulada tomará lugar en el vector ay se dibujarán los correlogramas asociados a este proceso.
set.seed(1)
alpha0 <- 0.1
alpha1 <- 0.4
beta1 <- 0.2
w <- rnorm(10000)
a <- rep(0, 10000)
h <- rep(0, 10000)
for (i in 2:10000){
h[i] <- alpha0 + alpha1*(a[i-1]^2)+ beta1*h[i-1]
a[i] <- w[i]*sqrt(h[i])
}El proceso tendrá el siguiente gráfico:
Los correlogramas serán los siguientes
La serie a exhibe características de una serie 𝐺𝐴𝑅𝐶𝐻 no correlacionada,pero el cuadrado de los valores de la misma presenta correlación. La siguiente sintaxis presenta el ajuste de la serie simulada haciendo uso de la función garch. La función por defecto es 𝐺𝐴𝑅𝐶𝐻(1,1), que suele proveer un modelo adecuado, pero modelos de orden superior pueden ser especificados con el parámetro order=c(p, q) para elegir a 𝑝 y a 𝑞.
## Warning: package 'tseries' was built under R version 4.0.3## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo##
## Call:
## garch(x = a, grad = "numerical", trace = FALSE)
##
## Coefficient(s):
## a0 a1 b1
## 0.09876 0.36659 0.24415En el ejemplo, se ha usado el parámetro Trace=FALSEpara suprimir la salida y el valor numérico de la gradiente grad=numerical, que es ligeramente más robusto que lo que está por defecto.
Bibliografía
Cowpertwait, P. S., & Metcalfe, A. V. (2009). Introductory time series with R. Springer Science & Business Media.
Shumway, R. H., & Stoffer, D. S. (2010). Time series analysis and its applications: with R examples. Springer Science & Business Media.