Trabajo Práctico Integrador - Series de Tiempo

Descripción de los datos

El ICCASal es una medida de frecuencia mensual, de disponibilidad periódica y frecuente, que permite realizar un monitoreo constante de la evolución de la actividad económica de la provincia de Salta.

El índice compuesto resume en una sola medida periódica el comportamiento de un conjunto de variables que se mueven de forma simultánea, y se relacionan de modo homogéneo y coincidente al cicl económico general.

Está confeccionado por 11 series económicas, representativas de los sectores de mayor importancia en la economía de Salta, que atravesaron un proceso de selección, cumpliendo los criterios de conformidad, sincronía temporal, significancia económica, adecuación estadística y actualización oportuna.

La información que proporciona el índice es de libre disponibilidad y está dirigida a un amplio conjunto de usuarios, abriendo las puertas para su uso por parte de diversos grupos de interés, incluyendo investigadores, economistas, así como personas que toman decisiones gubernamentales y empresariales. Además, permite su comparabilidad con otros indicadores de acitividad económica del mismo tenor que se desarrollan en otras provincias y regiones del país.

Objetivo de estudio	Actividad económica, ciclo económico provincial
Alcance geográfico	Provincia de Salta
Metodología	The Conference Board adaptada (NBER)
Frecuencia	Mensual
Serie de referencia	Estimador Mensual de Actividad Económica de Argentina
Periodo inicial	Enero 2004
Unidad de medida	Base 2004 = promedio 100
Disponibilidad	50-70 días

Carga de los datos

En primera instancia, se procede a cargar los datos de la serie de tiempo descargados desde la página web de la Dirección General de Estadísticas de Salta (http://estadisticas.salta.gov.ar/web/partes-de-prensa/6). Para llevar a cabo esta tarea se utiliza el paquete readr, asegurando que los datos se cargan correctamente.

# Cargamos la librería
library(readr)

## Warning: package 'readr' was built under R version 4.2.3

# Escribimos el directorio en donde está el archivo:
Dir <- paste0(dirname(rstudioapi::getActiveDocumentContext()$path),"/dataset/ICCASAL_2021.csv")
Data <- read.csv(Dir)

A continuación es necesario convertir los datos, que se almacenan como un simple vector, a una serie temporal, mediante la función base de R para conversión de series de tiempo, y estableciendo los parámetros comunicados por la DGE.

# Convertimos los datos en una serie temporal
Data <- ts(Data, start=c(2004, 1), frequency=12)

str(Data)

##  Time-Series [1:216, 1] from 2004 to 2022: 97.1 98.6 99.4 99.8 99.8 ...
##  - attr(*, "dimnames")=List of 2
##   ..$ : NULL
##   ..$ : chr "ICCASal"

Los datos se tratan de una serie de tiempo con 216 observaciones, que abarcan desde el mes de Enero del año 2004 hasta el año 2022. Es una serie de tiempo representada en una única columna, lo que sugiere que es univariante, es decir, un único valor registrado en cada punto temporal.

Análisis Exploratorio

En esta primera sección, nos enfocamos en el análisis exploratorio de datos aplicado a la serie de tiempo bajo estudio. El propósito de esta etapa es examinar de manera rigurosa la estructura y características temporales de los datos. A lo largo de esta sección, se abordan aspectos esenciales como la visualización de la serie de tiempo, la detección de valores atípicos, la evaluación de la estacionariedad y la exploración de la autocorrelación. Estos pasos sientan las bases para la preparación de los datos para futuros análisis de modelado y pronóstico.

Identificación de valores faltantes

# Conteo de valores faltantes
sum(is.na(Data))

## [1] 0

Como el resultado es 0, esto indica que no se han identificado valores faltantes en la serie de tiempo.

Ciclos

# Verificar los ciclos de la serie
cycle(Data)

##      Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
## 2004   1   2   3   4   5   6   7   8   9  10  11  12
## 2005   1   2   3   4   5   6   7   8   9  10  11  12
## 2006   1   2   3   4   5   6   7   8   9  10  11  12
## 2007   1   2   3   4   5   6   7   8   9  10  11  12
## 2008   1   2   3   4   5   6   7   8   9  10  11  12
## 2009   1   2   3   4   5   6   7   8   9  10  11  12
## 2010   1   2   3   4   5   6   7   8   9  10  11  12
## 2011   1   2   3   4   5   6   7   8   9  10  11  12
## 2012   1   2   3   4   5   6   7   8   9  10  11  12
## 2013   1   2   3   4   5   6   7   8   9  10  11  12
## 2014   1   2   3   4   5   6   7   8   9  10  11  12
## 2015   1   2   3   4   5   6   7   8   9  10  11  12
## 2016   1   2   3   4   5   6   7   8   9  10  11  12
## 2017   1   2   3   4   5   6   7   8   9  10  11  12
## 2018   1   2   3   4   5   6   7   8   9  10  11  12
## 2019   1   2   3   4   5   6   7   8   9  10  11  12
## 2020   1   2   3   4   5   6   7   8   9  10  11  12
## 2021   1   2   3   4   5   6   7   8   9  10  11  12

Esta salida indica que la serie de tiempo estudiada exhibe un patrón cíclico mensual.

Resúmen estadístico

# Resúmen descriptivo
summary(Data)

##     ICCASal      
##  Min.   : 97.12  
##  1st Qu.:138.88  
##  Median :145.53  
##  Mean   :144.63  
##  3rd Qu.:158.74  
##  Max.   :166.35

Este resultado revela que, a pesar de la amplia extensión entre el valor máximo y mínimo, la dispersión de los datos es relativamente baja, y esta observación se ve respaldada por la proximidad de la media y la mediana, ambas alrededor de 145.

Gráfico

# Librerías
library(ggplot2)

## Warning: package 'ggplot2' was built under R version 4.2.3

library(ggfortify)

## Warning: package 'ggfortify' was built under R version 4.2.3

# Gráfico
autoplot(Data) +
  labs(x="Fecha", y="ICCASal", title="ICCASal (2004-2022)")

La serie exhibe un crecimiento gradual a lo largo de la segunda mitad de la década de 2000, seguido de un patrón relativamente estable en la década de 2010, con picos notables en el año 2015, y una disminución considerable al finalizar esta década, alcanzando un mínimo relativo a comienzos del año 2020. Posteriormente retoma una tendencia positiva, llegando al mes de diciembre del año 2021 en donde alcanza valores similares a los del año 2019, implicando una leve recuperación superando los efectos del contexto de una fuerte incidencia de la pandemia de COVID-19 en la economía.

Estacionalidad de la serie

Para visualizar la estacionalidad de la serie temporal, es factible realizar una descomposición por ciclos, lo que nos permitirá observar cómo se comporta la serie en cada mes a lo largo de su evolución.

# Librerías
library(ggplot2)

# Armar un dataframe con la serie y sus ciclos
Data_df <- data.frame(Data,
                      Cycle=as.factor(cycle(Data)))

ggplot(data=Data_df, aes(x=Cycle, y=ICCASal)) +
  stat_boxplot(geom="errorbar", width=0.5) +
  geom_violin(alpha=0.2) + geom_boxplot(width=0.4) + 
  labs(x="Mes", title="Dispersión por ciclos mensuales del ICCAsal") + 
  geom_hline(aes(yintercept=mean(Data), color='Media'), linetype='dashed') +
  geom_hline(aes(yintercept=median(Data), color='Mediana'), linetype='dashed') +
  geom_hline(aes(yintercept=min(Data), color='Extremos'), linetype='dashed') +
  geom_hline(aes(yintercept=max(Data), color='Extremos'), linetype='dashed') +
  scale_color_manual(name = "Estadísticos", values = c(Mediana="blue", Media="red", Extremos="orange"))

Observamos en este gráfico varios aspectos importantes. En primer lugar, se identifican outliers inferiores, que son claramente los primeros valores de la serie y que destacan como valores inusuales. Las mayores variaciones se concentran en los cuartiles inferiores, con notables picos entre los meses de febrero y junio, donde la dispersión de los datos es más amplia, alcanzando valores cercanos a 110. Aunque las medianas se mantienen relativamente constantes a lo largo del año, se aprecia un ligero aumento en la mediana desde el mes de mayo. Asimismo, los cuartiles superiores exhiben un incremento notable en sus límites durante el mes de abril, que luego disminuyen en julio.

Se grafican la media (línea roja punteada), la mediana (línea azul punteada) y los valores extremos superiores e inferiores (líneas naranjas punteadas) con propósito ilustrativo.

Sin embargo, es importante señalar que este gráfico no muestra una estacionalidad clara, y de hecho, más bien sugiere una falta de patrón estacional, especialmente teniendo en cuenta los valores atípicos inferiores. La serie parece caracterizarse por su variabilidad y la ausencia de ciclos estacionales distintivos a lo largo del año.

Correlogramas

La representación gráfica de la función de autocorrelación simple (ACF) y parcial (PACF) resume la correlación de una serie temporal en varios retardos. A continuación, se realiza la representación de la ACF y PACF para la serie bajo estudio.

Correlograma de ACF

# Librerías
library(ggplot2)
library(forecast)

## Warning: package 'forecast' was built under R version 4.2.3

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

## Registered S3 methods overwritten by 'forecast':
##   method                 from     
##   autoplot.Arima         ggfortify
##   autoplot.acf           ggfortify
##   autoplot.ar            ggfortify
##   autoplot.bats          ggfortify
##   autoplot.decomposed.ts ggfortify
##   autoplot.ets           ggfortify
##   autoplot.forecast      ggfortify
##   autoplot.stl           ggfortify
##   autoplot.ts            ggfortify
##   fitted.ar              ggfortify
##   fortify.ts             ggfortify
##   residuals.ar           ggfortify

# Gráficos
ggAcf(Data) +
  labs(title="Serie: ICCASal")

En este caso, el correlograma de la función de autocorrelación simple indica una fuerte autocorrelación en la serie, valores significativamente altos que disminuyen lentamiente a medida que se analizan los retardos. Esto indica que la serie original dista mucho de ser completamente estocástica. Incluso es observable que no existe una relación clara cada 12 ciclos en la serie, lo cual podría haber indicado una estacionalidad. Por otra parte, al haber un decaimiento muy lento, es lógico suponer una falta de estacionariedad en la serie.

Correlograma de PACF

# Gráficos
ggPacf(Data) +
  labs(title="Serie: ICCASal")

En este contexto, dado que en la función de autocorrelación parcial se observa solamente un valor significativamente alto y el resto de los valores son notoriamente bajos y parecen seguir un patrón aleatorio.

Descomposición de la serie

A continuación podemos descomponer la serie para analizarla por sus componentes, utilizando la función base decompose().

# Librerías
library(cowplot)

## Warning: package 'cowplot' was built under R version 4.2.3

# Descomposiciones
Data_add <- decompose(Data, type="additive")
Data_mult <- decompose(Data, type="multiplicative")

# Gráficos
corr_data_add <- autoplot(Data_add) + labs(title="Descomposición de la serie aditiva")
corr_data_mult <- autoplot(Data_mult) + labs(title="Descomposición de la serie multiplicativa")

plot_grid(corr_data_add, corr_data_mult, nrow=1)

En estos gráficos de la serie descompuesta, podemos observar que ambas descomposiciones son similares entre sí, por lo que elegir uno por sobre otro es irrelevante. La parte estacional de la serie no es relevante, en ninguno de los dos tipos de descomposición, como se había observado a partir del gráfico de boxplot estacional, por lo cual no será necesaria una diferencia para corregir la estacionalidad.

Podemos observar los residuos con más detalle.

# Librería
library(cowplot)

# Descomponer la serie
Data_desc <- decompose(Data)
Data_random <- Data_desc$random

# Gráficos
corr_Data_random <- autoplot(Data_random) + labs(title="Componente aleatoria de la serie")
Data_random_acf <- ggAcf(Data_random) + labs(title="ACF de la componente aleatoria de la serie")
Data_random_pacf <- ggPacf(Data_random) + labs(title="PACF de la componente aleatoria de la serie")

plot_grid(corr_Data_random, Data_random_acf, Data_random_pacf, ncol=1)

En primer lugar, el gráfico de la componente aleatoria de la serie indica que estos residuos son estacionarios en media y en varianza, por lo tanto se puede realizar un análisis de los correlogramas. Aunque el gráfico de la componente aleatoria de la serie nos muestra una notable aleatoriedad a lo largo de la serie, el correlograma de la función de autocorrelación simple indica un valor significativamente alto en los primeros tres retardos, seguido de un comportamiento ondulatorio entre valores positivos y negativos, propio de un término autorregresivo en la serie. Por otra parte, el correlograma de la función de autocorrelación parcial muestra solo dos valores significativamente altos continuos, indicios de un modelo autorregresivo, como se había indicado anteriormente.

De todas las observaciones realizadas, la más importante es notar que la serie es no estacionaria, la cual es una característica a considerar para la modelización.

Identificación

El propósito de esta fase de modelización es seleccionar el modelo ARIMA(p,d,q) más adecuado para la serie, de manera que pueda reproducir con precisión sus características. La identificación de este modelo se realiza en dos etapas:

1. Análisis de la estacionariedad, que implica determinar las transformaciones necesarias para lograr que la serie sea estacionaria. Esta etapa consta de dos partes:

• Estacionariedad en la varianza: se buscan transformaciones que estabilicen la varianza de la serie.

• Estacionariedad en la media: se determina el número de diferencias \(d\) necesarias para hacer que la serie sea estacionaria en términos de su media.

2. Selección de los órdenes \(p\) y \(q\). Una vez que se ha logrado estacionarizar la serie, el objetivo es identificar, de ser posible, un proceso estacionario ARMA(p,q) que la generó.

Estacionariedad en varianza

Una serie se considerará estacionaria en términos de varianza cuando se cumpla la suposición de que existe una única varianza constante para toda la serie temporal. Esto implica que la variabilidad de la serie alrededor de su media se mantenga constante a lo largo del tiempo. Si la serie no cumple con este requisito de estacionariedad en varianza, se recurre a las transformaciones destinadas a estabilizar la varianza, como las transformaciones Box-Cox,

\[ Y_t^{(\lambda)} = \begin{cases} \frac{Y_t^{\lambda} - 1}{\lambda} & \lambda \ne 0\\ \ln Y_t & \lambda = 0 \end{cases} \]

Las transformaciones Box-Cox incluye una familia de funciones: raíz cuadrada, inversa, etc. En el contexto de las series económicas, es común que estas sean siempre positivas y no contengan valores iguales a cero. Por lo tanto, en primer lugar, podemos aplicar la transformación de Box-Cox para valores distintos de lambda. El paquete forecast dispone de una función llamada BoxCox() que facilita esta transformación, permitiendo especificar el valor de lambda deseado.

# Librerías
library(cowplot)

# Valores de lambda
lambda_1 <- -2
lambda_2 <- -1
lambda_3 <- 0
lambda_4 <- 1
lambda_5 <- 2

# Transformaciones
Data_boxcox_1 <- BoxCox(Data, lambda=lambda_1)
Data_boxcox_2 <- BoxCox(Data, lambda=lambda_2)
Data_boxcox_3 <- BoxCox(Data, lambda=lambda_3)
Data_boxcox_4 <- BoxCox(Data, lambda=lambda_4)
Data_boxcox_5 <- BoxCox(Data, lambda=lambda_5)

# Gráficos de la serie transformada
gg_Data_boxcox_1 <- autoplot(Data_boxcox_1)
gg_Data_boxcox_2 <- autoplot(Data_boxcox_2)
gg_Data_boxcox_3 <- autoplot(Data_boxcox_3)
gg_Data_boxcox_4 <- autoplot(Data_boxcox_4)
gg_Data_boxcox_5 <- autoplot(Data_boxcox_5)

plot_grid(gg_Data_boxcox_1, gg_Data_boxcox_2, gg_Data_boxcox_3, gg_Data_boxcox_4, gg_Data_boxcox_5, ncol=1)

Pero, por otro lado, el mismo paquete forecast proporciona la función BoxCox.lambda() que permite seleccionar automáticamente un parámetro ajustado según un método establecido. Si el método es "loglik", el valor de lambda se elige para maximizar la verosimilitud logarítmica de un modelo lineal ajustado. Para datos no estacionales se ajusta una tendencia temporal lineal, mientras que para datos estacionales, se utiliza una tendencia temporal lineal con variables estacionales.

BoxCox.lambda(Data, method="loglik", lower=-2)

## [1] 2

Lo cual nos asegura que el mejor parámetro que podríamos escoger es 2, correspondiente a la última serie graficada anteriormente.

Estacionariedad en media

En esta sección, se busca determinar si la serie es estacionaria en términos de media, es decir, si oscila alrededor de un nivel constante o no. Para tomar esta decisión, se consideran las características distintivas entre series estacionarias y no estacionarias.

Características de una serie estacionaria:

• Una serie se considera estacionaria en términos de media cuando se puede asumir que existe una única media constante para toda la serie temporal, es decir, cuando fluctúa alrededor de una media constante.

• La función de autocorrelación teórica de un proceso estacionario en media decae rápidamente, siguiendo un patrón exponencial.

Características de una serie no estacionaria:

• La serie no es estacionaria cuando presenta tendencia o varios tramos con medias diferentes.

• En general, un proceso con alguna raíz unitaria presenta una función de autocorrelación muestral con un decaimiento muy lento, no siendo necesario que se mantenga próximo a la unidad.

Si la serie no es estacionaria en términos de media, se puede lograr la estacionariedad mediante la aplicación de diferenciación. En consecuencia, si la serie no es estacionaria en media, se aplicarán \(d\) diferenciaciones sucesivas de orden 1 a la serie hasta obtener una serie estacionaria:

\[ Z_t = (1 - B)^{d} Y_t \]

Es crucial la precaución al identificar el número exacto de diferencias necesario para lograr que la serie sea estacionaria, ya que pueden surgir diversos problemas, incluida la sobrediferenciación. La serie se considera sobrediferenciada cuando se aplican más diferenciaciones de las requeridas. Por ejemplo, esto ocurre si se selecciona un orden de integración \(d\) cuando la serie \(\nabla^{d-1} Y_t\) ya es estacionaria.

En este caso, es importante tener en cuenta que si se diferencia un proceso que ya es estacionario, este continuará siendo estacionario. El objetivo en esta etapa de la modelización ARIMA es determinar el número mínimo de diferencias \(d\) necesario para transformar la serie en una serie estacionaria. Por otra parte, a partir de la observación empírica, se ha observado que las series económicas generalmente requieren órdenes de diferenciación 0, 1 o 2.

La función base de R, diff(), permite realizar múltiples diferenciaciones estableciendo el orden de estas.

# Librerías
library(cowplot)

# Transformación BoxCox
lambda <- BoxCox.lambda(Data, method="loglik")
Data_tr <- BoxCox(Data, lambda=lambda)

# Diferenciaciones
Data_diff_1 <- diff(Data_tr, 1) # Primera diferenciación
Data_diff_2 <- diff(Data_tr, 2) # Segunda diferenciación

# Gráficos
gg_Data_diff_0 <- autoplot(Data) + labs(x="Tiempo", title="Serie original")
gg_Data_diff_1 <- autoplot(Data_diff_1) + labs(x="Tiempo", title="Primera diferenciación de la serie, d=1")
gg_Data_diff_2 <- autoplot(Data_diff_2) + labs(x="Tiempo", title="Segunda diferenciación de la serie, d=2")
plot_grid(gg_Data_diff_0, gg_Data_diff_1, gg_Data_diff_2, ncol=1)

A partir de este resultado, se observa que la primera diferenciación es suficiente para transformar la serie en estacionaria en términos de media, ya que sus valores fluctúan alrededor de una constante. Podemos apoyar esta observación generando los correlogramas correspondientes:

# Librerías
library(cowplot)

# Transformación BoxCox
lambda <- BoxCox.lambda(Data, method="loglik")
Data_tr <- BoxCox(Data, lambda=lambda)

# Diferenciación
Data_diff <- diff(Data_tr, 1)

# Gráficos
gg_Data_diff_0 <- ggAcf(Data) + labs(title="ACF de la serie original")
gg_Data_diff_1 <- ggAcf(Data_diff_1) + labs(title="ACF de la primera diferenciación de la serie, d=1")
gg_Data_diff_2 <- ggAcf(Data_diff_2) + labs(title="ACF de la primera diferenciación de la serie, d=2")
plot_grid(gg_Data_diff_0, gg_Data_diff_1, gg_Data_diff_2, ncol=1)

Los correlogramas no respaldan por completo nuestra observación, si bien en ambos correlogramas de las diferencias se observa valores significativamente altos a partir del cuarto o quinto retardo, en el correlograma se observa una disminución exponencial en la segunda diferenciación, mientras que no es del todo clara la disminución de los valores en el correlograma de la primera diferenciación. Por lo tanto decidimos utilizar una segunda diferenciación.

Por otra parte, los constrastes de raíces unitarias proporcionan contrastes estadísticos que permiten a partir del conjunto de información hacer inferencia sobre la existencia o no de una raíz unitaria en una serie, es decir, sobre la no estacionariedad de la serie. Si se rechaza la hipótesis nula de existencia de raíz unitaria en \(\nabla^{d-1} Y_t\), no se diferenciará más la serie. En caso contrario, si no se rechaza la hipótesis nula se tomará una diferencia más de orden 1.

El contraste de Dickey-Fuller es una generalización de un contraste para un modelo AR(1), sabiendo que la serie puede seguir procesos más generales ARMA(p,q). Como es bien conocido, todo proceso ARMA(p,q) se puede aproximar hasta el grado de bondad necesario mediante un AR(p): \[ Y_t = \phi_1 Y_{t-1} + \phi_1 Y_{t-2} + ... + \phi_p Y_{t-p} + a_t \qquad a_t \sim NIID(0, \sigma^2) \]

Este modelo se puede reparametrizar como sigue: \[ \nabla Y_t = \beta Y_{t-1} + \alpha_1 \nabla Y_{t-1} + ... + \alpha_{p-1} \nabla Y_{t-p+1} + a_t \]

donde \[ \beta = \sum_{i = 1}^{p} \phi_i - 1 \qquad \text{y} \qquad \alpha_i = \sum_{j=1}^{i} \phi_{p-i+j} \]

Dado que un proceso AR(p) tiene una raíz unitaria cuando \(\sum_{i=1}^{p} \phi_i = 1\), contrastar la hipótesis nula de existencia de raíz unitaria es equivalente a contrastar la \(H_0: \beta = 0\) en la regresión anterior. Este contraste se basa en la estimación por mínimos cuadrados ordinarios del parámetro \(\beta\) en el modelo y en el correspondiente estadístico \(t\).

La librería tseries contiene una función para realizar este contraste de hipótesis para una serie establecida. Por ejemplo, para la serie original:

# Librerías
library(tseries)

## Warning: package 'tseries' was built under R version 4.2.3

# Test de Dickey-Fuller aumentado
adf.test(Data)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  Data
## Dickey-Fuller = -2.1963, Lag order = 5, p-value = 0.4937
## alternative hypothesis: stationary

Una “rule of thumb” nos permite concluir que, como el p-value es mucho mayor que 0.05, no existe evidencia suficiente para rechazar la hipótesis nula, y concluimos que la serie original es no estacionaria.

Y el mismo contraste para la serie transformada, con una diferenciación de \(d = 2\).

# Librerías
library(tseries)

# Transformación BoxCox
lambda <- BoxCox.lambda(Data, method="loglik")
Data_tr <- BoxCox(Data, lambda=lambda)

# Diferenciación
Data_tr <- diff(Data_tr, 2)

# Test de Dickey-Fuller aumentado
adf.test(Data_tr)

## Warning in adf.test(Data_tr): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  Data_tr
## Dickey-Fuller = -5.3762, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary

Utilizando la misma regla propuesta anteriormente, al observar que el p-value de \(0.01 < 0.05\), podemos concluir que existe evidencia suficiente al nivel de significancia establecido para rechazar la hipótesis nula, aceptamos la hipótesis alternativa, y concluimos que la serie transformada es estacionaria.

El análisis anterior nos permite concluir que, si bien la serie original no es estacionaria ni en varianza ni en media, sí será estacionaria la siguiente serie transformada: \[ Z_t = \nabla_{12} \nabla^2 \biggl( \frac{Y_t^2 - 1}{2} \biggl) \]

Identificación del modelo estacionario

Una vez que hemos determinado el orden de diferenciación \(d\), procedemos a obtener la transformación estacionaria de la serie \(Z_t\) como \((1 - B)^{d} Y_t\). Esta transformación nos permite representar la serie como un proceso ARMA(p,q) estacionario. Durante esta etapa, nuestro objetivo es identificar los valores de \(p\) y \(q\) que mejor capturen las características de la serie estacionaria, además de analizar la inclusión del parámetro asociado a la media, \(\delta\).

La función de autocorrelación, ACF, juega un papel fundamental en la descripción de las propiedades dinámicas del proceso estacionario. Se utiliza el ACF como la herramienta principal para determinar los valores óptimos de \(p\) y \(q\) en el modelo ARMA que mejor se ajuste a las características de la serie estacionaria \(Z_t\).

Para identificar los ordenes \(p\) y \(q\), se comparan las funciones de autocorrelación muestrales simples y parciales con las ACF y PACF teóricas de los modelos ARMA, cuyas características son conocidas. De todos modos, es importante destacar que este método no se trata de identificar el modelo correcto, sino de acotar un subconjunto de modelos ARIMA que han podido generar la serie. Posteriormente, se puede volver a plantear la identificación del proceso de forma iterativa hasta obtener resultados aceptables.

# Librerías
library(cowplot)

# Transformación BoxCox
lambda <- BoxCox.lambda(Data, method="loglik")
Data_tr <- BoxCox(Data, lambda=lambda)

# Diferenciación
Data_tr <- diff(Data_tr, 2)

# Gráficos
Data_tr_acf <- ggAcf(Data_tr) + labs(title="ACF de la serie transformada")
Data_tr_pacf <- ggPacf(Data_tr) + labs(title="PACF de la serie transformada")
plot_grid(Data_tr_acf, Data_tr_pacf, ncol=1)

El correlograma ACF muestral de la serie \(Z_t\) presenta una estructura infinita decreciendo en forma de onda amortiguada, lo cual sugiere un modelo AR(p) o ARMA(p,q) con parte autorregresiva de orden mayor o igual que 2. El análisis del correlograma PACF muestra un valor significativamente grande en los primeros retardos, siendo estos valores positivos y negativos, seguido también de una estructura finita que decrece en forma amortiguada alternándose entre valores positivos y negativos. Esta observación indica un término MA(q) de valor alto en módulo, cuyos coeficientes probablemente sean negativos.

Según esta interpretación, y teniendo en cuenta que tanto el ACF como el PACF tienen valores significativos en la misma cantidad de retardos, lo correcto sería pensar que existe al menos un modelo ARMA(2,2) que se puede ajustar a la serie, pero también es posible un modelo ARMA(2,3).

Estimación

Una vez identificado el proceso estocástico que ha podido generar la serie temporal \(Y_t\), la siguiente etapa consiste en estimar los parámetros desconocidos de dichos modelos: \[ \beta = (\delta, \phi_1, ..., \phi_p, \theta_1, ..., \theta_q)' \qquad \text{y} \qquad \sigma^2_a \]

Estos parámetros pueden ser estimados de diversas formas, pero las funciones base de R, arima() y arma() permiten ajustar un modelo a la serie transformada, estableciendo los parámetros \(p\), \(d\), y \(q\) que se identificaron previamente.

En este paso, también tenemos la oportunidad de comparar varios modelos con diferentes conjuntos de parámetros y tomar la decisión sobre cuál es el más adecuado. Para esta tarea, podemos hacer uso del criterio de información de Akaike (AIC), el cual evalúa la bondad de ajuste del modelo ARIMA(p,d,q). El AIC es una medida relativa de la calidad de un modelo estadístico en relación a un conjunto específico de datos, y nos proporciona una herramienta para la selección del modelo más apropiado.

En el caso general, el AIC es \[ AIC = 2 - 2 \ln (L) \] donde \(k\) es el número de parámetros en el modelo estadístico, y \(L\) es el máximo valor de la función de verosimilitud para el modelo estimado.

Entre un conjunto de modelos candidatos para los datos, el modelo preferido es aquel que presenta el valor más bajo en el criterio de información de Akaike (AIC). El AIC no solo refleja la calidad del ajuste, sino que también incorpora una penalización que aumenta con el número de parámetros estimados. Esta penalización tiene la finalidad de desalentar el sobreajuste, es decir, la tendencia de aumentar el número de parámetros libres en el modelo para mejorar el ajuste, sin considerar el número real de parámetros involucrados en el proceso de generación de datos.

# Estimaciones ARIMA(p,d,q)
Data_arima_1 <- Arima(Data, order=c(1,2,1))
Data_arima_2 <- Arima(Data, order=c(1,2,2))
Data_arima_3 <- Arima(Data, order=c(1,2,3))
Data_arima_4 <- Arima(Data, order=c(2,2,1))
Data_arima_5 <- Arima(Data, order=c(2,2,2))
Data_arima_6 <- Arima(Data, order=c(2,2,3))

# Valores del criterio de Akaike
paste0("Para un modelo ARIMA(1,2,1), el valor del AIC es: ", Data_arima_1$aic)

## [1] "Para un modelo ARIMA(1,2,1), el valor del AIC es: 612.327240067026"

paste0("Para un modelo ARIMA(1,2,2), el valor del AIC es: ", Data_arima_2$aic)

## [1] "Para un modelo ARIMA(1,2,2), el valor del AIC es: 611.990002841039"

paste0("Para un modelo ARIMA(1,2,3), el valor del AIC es: ", Data_arima_3$aic)

## [1] "Para un modelo ARIMA(1,2,3), el valor del AIC es: 609.633084396685"

paste0("Para un modelo ARIMA(2,2,1), el valor del AIC es: ", Data_arima_4$aic)

## [1] "Para un modelo ARIMA(2,2,1), el valor del AIC es: 613.713160616812"

paste0("Para un modelo ARIMA(2,2,2), el valor del AIC es: ", Data_arima_5$aic)

## [1] "Para un modelo ARIMA(2,2,2), el valor del AIC es: 610.219607477346"

paste0("Para un modelo ARIMA(2,2,3), el valor del AIC es: ", Data_arima_6$aic)

## [1] "Para un modelo ARIMA(2,2,3), el valor del AIC es: 609.700661273426"

El análisis de los resultados del criterio de información de Akaike (AIC) revela que los modelos ARIMA(1,2,3) y ARIMA(2,2,3) presentan los valores más bajos de AIC. La decisión de un modelo ARMA(2,2,3) por sobre el otro se basa exclusivamente en el correlograma mostrado anteriormente.

A continuación, podemos analizar las estimaciones generadas por el ajuste del modelo establecido.

# Estimación ARIMA(2,2,3)
Arima(Data, order=c(2,2,3))

## Series: Data 
## ARIMA(2,2,3) 
## 
## Coefficients:
##          ar1     ar2      ma1      ma2     ma3
##       0.3693  0.2466  -0.7631  -0.4954  0.2807
## s.e.  0.3158  0.1662   0.3094   0.2553  0.1208
## 
## sigma^2 = 0.9708:  log likelihood = -298.85
## AIC=609.7   AICc=610.11   BIC=629.9

El modelo estimado para la serie es el siguiente: \[ z_t = \delta + 0.37 z_{t-1} + 0.25 z_{t-2} + a_t - 0.76 a_{t-1} - 0.50 a_{t-2} + 0.28 a_{t-3} \]

Para estimar la constante \(\delta\) es preferible utilizar la estimación dada por Gretl, que es la media estimada de la serie estacionaria. Si el modelo estimado es AR(p) o ARMA(p,q): \[ \hat{\mu} = \frac{\hat{\delta}}{1 - \hat{\phi_1} - \hat{\phi_2}} = \frac{\hat{\delta}}{1-0.37 - 0.25} \] En donde \(\hat{\mu}\) es:

round(mean(Data),2)

## [1] 144.63

Y despejando, \[ \hat{\delta} = 144.63 \cdot 0.38 = 54.96 \]

Finalmente, el modelo completo será \[ z_t = 54.96 + 0.37 z_{t-1} + 0.25 z_{t-2} + a_t - 0.76 a_{t-1} - 0.50 a_{t-2} + 0.28 a_{t-3} \]

Validación del modelo

En la etapa de validación o diagnosis se procede a evaluar la adecuación de los modelos estimados a los datos. Se tiene en cuenta:

1. Se verifica si las estimaciones de los coeficientes del modelo son estadísticamente significativas y si cumplen con las condiciones de estacionariedad e invertibilidad que son requisitos esenciales para los parámetros del modelo.

2. Se analiza si los residuos del modelo exhiben un comportamiento similar al de las innovaciones, lo que significa que deben comportarse como ruido blanco.

Análisis de coeficientes estimados

Sabiendo que los coeficientes estimados son significativamente grandes, podemos verificar si cumplen las condiciones de estacionariedad e invertibilidad. Para eso tenemos que verificar si los módulos de las raíces del polinomio autorregresivo y de medias móviles están por fuera del círculo unitario.

Para el polinomio autorregresivo, podemos reacomodar el modelo:

\[ (1 - 0.37 B - 0.25 B^2) z_t = ... \]

La función polyroot() verificar las raíces de este polinomio

polyroot(c(1,-0.37,-0.25))

## [1]  1.39251-0i -2.87251+0i

Y la función uc.check() del paquete UnitCircle nos permite graficar el círculo unitario para tener una mejor visualización de las raíces.

library(UnitCircle)

uc.check(pol=c(1,-0.37,-0.25), print_output=FALSE)

Como ambas raíces se encuentran por fuera del círculo unitario, concluimos que la parte autorregresiva cumple con la condición de estacionariedad.

Para el polinomio de medias móviles, podemos reacomodar el modelo: \[ ... = (1 - 0.76 B - 0.50 B^2 + 0.28 B^3) a_t \]

La función polyroot() verificar las raíces de este polinomio

polyroot(c(1, -0.76, -0.50, 0.28))

## [1]  1.021920-0i -1.526156+0i  2.289951+0i

Y la función uc.check() del paquete UnitCircle nos permite graficar el círculo unitario para tener una mejor visualización de las raíces.

library(UnitCircle)

uc.check(pol=c(1, -0.76, -0.50, 0.28), print_output=FALSE)

Por lo tanto, según los resultados anteriores, podemos concluir que el modelo cumple con las condiciones de estacionariedad e invertibilidad.

Análisis de residuos

Si el modelo ARMA(p,q) elegido para la serie estacionaria \(Z_t\), por ejemplo \[ \Phi_p (B) Z_t = \Theta_q (B) a_t \] es correcto, entonces \[ a_t = \frac{\Phi_p(B)}{\Theta_q (B)} Z_t \] es un proceso ruido blanco, y los residuos de la estimación de este modelo son estimaciones de \(a_t\). El análisis de residuos consiste en una serie de contrastes de diagnóstico con el objetivo de determinar si los residuos replican el comportamiento de un ruido blanco, es decir, si su media es cero, su varianza constante y las autocorrelaciones nulas.

La función ggtsdiag() del paquete ggfortify nos permite realizar un diagnóstico del modelo.

# Librerías
library(ggfortify)

# Estimación ARIMA(1,1,2)
Data_arima <- Arima(Data, order=c(2,2,3))

# Diagnóstico de residuos
ggtsdiag(Data_arima)

checkresiduals(Data_arima)

## 
##  Ljung-Box test
## 
## data:  Residuals from ARIMA(2,2,3)
## Q* = 19.962, df = 19, p-value = 0.3969
## 
## Model df: 5.   Total lags used: 24

Estos diagnóstico arroja resultados interesantes. En primer lugar observamos que los residuos (Standarized Residuals) presentan una media en 0, y una varianza casi constante, por lo menos en casi toda la serie residual, y su histograma muestra que se distribuyen bastante normal. Por otra parte, el ACF de los residuos muestra un único valor significativamente grande, correspondiente a la autocorrelación de orden 1, y muchos valores significativamente pequeños que se los puede interpretar como valores aleatorios, respondiendo a un ruido blanco.

Finalmente, el test de Ljung-Box nos permite comprobar analíticamente mediante un constraste si es que existe autocorrelación en la serie, presentando como hipótesis nula que los residuos se distribuyen de manera independiente, y como hipótesis alternativa que los residuos no se distribuyen independientemente, y por tanto exhiben correlación. Valores p demasiado grandes a 0.05 (línea azul punteada) nos indican que no existe evidencia suficiente para rechazar la hipótesis nula, por lo tanto debemos concluir que los residuos se distribuyen de manera independiente.

Como conclusión, el modelo elegido es un buen ajuste a la serie temporal.

Predicción

El objetivo es obtener predicciones óptimas de \(Y_t\) en algún momento futuro basadas en un conjunto de información dado, que en el caso del análisis de series temporales univariante está formado por el pasado disponible de la serie temporal \[ I_T = \{ Y_t, Y_{T-1}, Y_{T-2}, ... \} \]

Supongamos que se observa una serie temporal denotada por \(Y_t\), para \(t=0\) hasta \(t=T\), donde \(T\) es la última observación de que disponemos. La tarea de predecir series temporales implica anticipar el valor que tomará la serie en momentos futuros \(T + l\), donde \(l\) representa el número de períodos en el futuro que estamos considerando. Dado que estamos asumiendo que la serie temporal es una realización de un proceso estocástico estacionario, el valor que deseamos predecir, \(Y_{T+l}\), también es una variable aleatoria.

La predicción de una variable aleatoria implica la predicción de su función de distribución, aunque en la mayoría de los casos, determinar completamente la forma de la función de densidad resulta extremadamente difícil sin realizar suposiciones poco realistas. Un objetivo más realista y práctico es diseñar intervalos de confianza alrededor del valor \(Y_{T+l}\). Es importante diferenciar entre la predicción por intervalos, que se enfoca en la construcción de estos intervalos de predicción, y la predicción puntual, que implica simplemente asignar un único valor a \(Y_{T+l}\) que de alguna manera represente toda la distribución de posibles valores.

De manera práctica, el paquete forecast permite mediante una función homónima, realizar una predicción.

# Librerías
library(forecast)

# Estimación ARIMA(2,2,3)
Data_arima <- Arima(Data, order=c(2,2,3))

# Predicción
Data_for <- forecast(Data_arima, level=c(80,95), h=12)

# Gráfico
autoplot(Data_for) +
  labs(title="Modelo ARIMA(2,2,3) ajustado a la serie", x="Tiempo", y="ICCASal")

Esta predicción se caracteriza por su enfoque conservador en múltiples aspectos, lo que resulta en una proyección que se asemeja a una línea prácticamente recta con una leve tendencia negativa. Las áreas azules indican las estimaciones por intervalos que se mencionaron anteriormente, en este caso el área más clara representa una estimación del 95% de confianza, mientras que la más oscura indica una estimación del 80% de confianza.

Ajuste automático del modelo

La librería forecast incluye la función auto.arima(), diseñada con la finalidad de efectuar una calibración automática de un modelo ARIMA en una serie temporal. Para llevar a cabo esta tarea, se requiere la especificación de los valores máximos para el ajuste del modelo, la determinación de la presencia o ausencia de patrones estacionales en la serie temporal, la elección de un criterio para el ajuste del modelo y, opcionalmente, la asignación de un valor para el parámetro \(\lambda\) que permita realizar una transformación Box-Cox si se considera necesario.

# Librería forecast
library(forecast)

auto.arima(Data, max.p=4, max.q=4, max.P=4, max.Q=4, max.d=3, max.order=5, seasonal=FALSE, ic='aic', lambda=2)

## Series: Data 
## ARIMA(2,2,2) 
## Box Cox transformation: lambda= 2 
## 
## Coefficients:
##           ar1     ar2      ma1      ma2
##       -0.2252  0.3358  -0.1591  -0.7918
## s.e.   0.1984  0.1416   0.1737   0.1709
## 
## sigma^2 = 19619:  log likelihood = -1360.23
## AIC=2730.45   AICc=2730.74   BIC=2747.28

La discrepancia en la elección de modelos se debe a la naturaleza conservadora de la función auto.arima(), la cual tiende a evitar la selección de modelos cuyas raíces estén próximas al círculo unitario en el espacio de parámetros. Esto se hace para garantizar la estacionariedad e invertibilidad del modelo, que son requisitos fundamentales en el análisis de series temporales. Aunque esta enfoque es prudente para evitar problemas de convergencia y estabilidad en la estimación de parámetros, también puede llevar a la selección de modelos más simples o menos flexibles de lo que se podría obtener a través de un enfoque más manual.

Predicción con el modelo ajustado automático

Podemos intentar generar una predicción utilizando el modelo automático ajustado por la función auto.arima().

# Librerías
library(forecast)

# Estimación ARIMA(2,2,3)
Data_arima <- auto.arima(Data, max.p=4, max.q=4, max.P=4, max.Q=4, max.d=3, max.order=5, seasonal=FALSE, ic='aic', lambda=2)

# Predicción
Data_for <- forecast(Data_arima, level=c(80,95), h=12)

# Gráfico
autoplot(Data_for) +
  labs(title="Modelo ARIMA(2,2,2) ajustado a la serie", x="Tiempo", y="ICCASal")

En este caso vemos que la predicción es muy similar a la realizada con el modelo encontrado manualmente, aunque en este caso la predicción resulta en una tendencia un poco más negativa que la anterior.

Es también interesante experimentar qué sucede cuando establecemos que la serie contiene un comportamiento estacional, a pesar que no fue observado en todo el análisis previo.

auto.arima(Data, max.p=4, max.q=4, max.P=4, max.Q=4, max.d=3, max.order=5, seasonal=TRUE, ic='aic', lambda=2)

## Series: Data 
## ARIMA(0,2,2)(0,0,2)[12] 
## Box Cox transformation: lambda= 2 
## 
## Coefficients:
##           ma1      ma2     sma1     sma2
##       -0.4076  -0.4548  -0.1645  -0.1050
## s.e.   0.0664   0.0768   0.0698   0.0781
## 
## sigma^2 = 19607:  log likelihood = -1360.23
## AIC=2730.46   AICc=2730.75   BIC=2747.29

Es curioso que se ajusta un modelo un poco más complejo, sin parte autorregresiva, y con modelo subyacente ARIMA(0,0,2) que ajusta la estacionalidad que se comporta, en este caso, de manera anual.

# Librerías
library(forecast)

# Estimación ARIMA(2,2,3)
Data_arima <- auto.arima(Data, max.p=4, max.q=4, max.P=4, max.Q=4, max.d=3, max.order=5, seasonal=TRUE, ic='aic', lambda=2)

# Predicción
Data_for <- forecast(Data_arima, level=c(80,95), h=12)

# Gráfico
autoplot(Data_for) +
  labs(title="Modelo ARIMA(0,2,2)(0,0,2)[12] ajustado a la serie", x="Tiempo", y="ICCASal")

Ahora la predicción es completamente diferente a la generada anteriormente, en este caso tiene una tendencia ligeramente positiva, y se observa un crecimiento en cierto periodo.

Conclusión

# Librerías
library(forecast)

# Estimación ARIMA(2,2,3)
Data_arima <- Arima(Data, order=c(2,2,3))

# Predicción
Data_for <- forecast(Data_arima, level=c(80,95), h=12)

# Gráfico
autoplot(Data_for) +
  labs(title="Modelo ARIMA(2,2,3) ajustado a la serie", x="Tiempo", y="ICCASal") +
  xlim(2019,2023) +
  ylim(120,160)

## Scale for x is already present.
## Adding another scale for x, which will replace the existing scale.

## Warning: Removed 180 rows containing missing values (`geom_line()`).

Durante los primeros meses del año 2021, la economía provincial experimentó una serie de contracciones, manifestándose con caídas moderadas en los niveles de actividad. No obstante, a partir de junio de 2021, se inició una secuencia continua de mejoras, y en diciembre de ese año, se alcanzó el nivel más alto de actividad desde mediados de 2019. Esto sugiere que la actividad económica gradualmente se encamina hacia una recuperación, superando los efectos de la pandemia. Las proyecciones indican que el índice podría mostrar una ligera mejoría o mantenerse en un patrón casi constante con una pequeña fluctuación negativa, pero es poco probable su recuperación inmediata y alcanzar valores semejantes a los que alcanzó antes de la pandemia.