library(tseries)
library(forecast)
library(fUnitRoots)
library(nortest)
library(dplyr)
library(lubridate)

Componentes de una Serie de Tiempo

Tendencia (Trend):
- La línea de tendencia muestra la dirección general del tiempo. En este gráfico, la tendencia es claramente ascendente, lo que indica que los valores observados aumentan con el tiempo.
Estacionalidad (Seasonal):
- La componente estacional captura patrones recurrentes a intervalos regulares, como estaciones del año, meses o días de la semana. En el gráfico, se observa una patrón repetitivo anual, indicando una fuerte estacionalidad con picos y valles que se repiten cada año.
Aleatorio (Random/Residual):
- La componente aleatoria muestra las fluctuaciones impredecibles o ruido en los datos después de remover la tendencia y la estacionalidad. Este componente captura variaciones que no siguen un patrón específico y puede deberse a factores externos o errores de medición.

pasajeros = AirPassengers
plot(decompose(pasajeros))

Ruido Blanco

El ruido blanco en series temporales es una secuencia de datos aleatorios que no contiene patrones discernibles ni estructura predecible. Se utiliza como una referencia ideal en el análisis de series temporales para determinar si los datos residuales de un modelo han capturado adecuadamente todas las tendencias y componentes estacionales.

Propiedades del Ruido Blanco

Media Cero:
- Los valores promedio de las observaciones son cero. En un conjunto suficientemente grande de datos, las desviaciones positivas y negativas se equilibran entre sí.
Varianza Constante:
- La varianza de las observaciones es constante a lo largo del tiempo, lo que indica que la dispersión de los datos alrededor de la media no cambia.
No Autocorrelación:
- Las observaciones no están correlacionadas entre sí. Esto significa que no hay patrones predecibles en los datos, y el valor de una observación en un punto en el tiempo no proporciona información sobre los valores en otros puntos en el tiempo.

# Generamos una serie de ruido blanco de 1000 observaciones.
n <- 1000
whitenoise <- rnorm(n, mean = 0, sd = 1)

plot(whitenoise, type = "l", 
     main = "White Noise", 
     xlab = "Tiempo", ylab = "Valor")

Random Walk

El “random walk” o caminata aleatoria es un concepto utilizado en estadística y teoría de probabilidades para describir una secuencia de pasos aleatorios. En el contexto de series temporales, un random walk es un proceso donde cada valor en la serie es igual al valor anterior más un cambio aleatorio. Este concepto se usa frecuentemente para modelar fenómenos como los precios de las acciones, donde se supone que los cambios diarios son impredecibles y no están influenciados por movimientos anteriores.

Propiedades del Random Walk

Dependencia en el Valor Anterior:
- Cada observación en una caminata aleatoria depende del valor inmediatamente anterior más un término de error aleatorio. Matemáticamente, puede representarse como \(X_t = X_{t-1} + \epsilon_t\), donde \(\epsilon_t\) es un término de error aleatorio con media cero.
No Estacionariedad:
- La caminata aleatoria no es estacionaria, lo que significa que sus propiedades estadísticas, como la media y la varianza, cambian con el tiempo. La varianza de los valores tiende a aumentar con el tiempo, lo que refleja una mayor dispersión a medida que avanza la serie.
Martingala:
- Un random walk es un tipo de martingala, lo que implica que el valor esperado del próximo paso es igual al valor actual. En otras palabras, no hay “memoria” en el proceso, y los pasos futuros no se ven afectados por los pasos pasados.
Impredecibilidad:
- Debido a su naturaleza aleatoria, no es posible predecir con precisión los valores futuros de una caminata aleatoria basándose en los valores pasados. Este es un principio clave en finanzas y otros campos donde los modelos de random walk se utilizan para representar comportamientos impredecibles.

# Generamos una serie de random walk de 1000 observaciones.
n <- 1000
rwalk <- NULL
rwalk[1] <- rnorm(1, mean = 0, sd = 1)

for (i in 2:n) {
  rwalk[i] <- rwalk[i-1] + rnorm(1, mean = 0, sd = 1) 
}

# Generamos una serie de random walk con drift de 1000 observaciones.
rwalk_drift <- NULL
rwalk_drift[1] <- rwalk[1]

for (i in 2:n) {
  rwalk_drift[i] <- rwalk_drift[i-1] + rnorm(1, mean = 0, sd = 1) + 0.1
}

plot(rwalk, type = "l", main = "Random Walk", 
     xlab = "Tiempo", ylab = "Valor", ylim = c(-20,50))
lines(rwalk_drift, type = "l", col = "red")
legend("topleft",legend = c("Sin Drift", "Con Drift"),
       col = c("black","red"),
       pch = 15)

Prueba de Dickey-Fuller Aumentada

La prueba de Dickey-Fuller Aumentada (ADF, por sus siglas en inglés) es un test estadístico utilizado para determinar si una serie temporal es estacionaria o no. La estacionariedad es una propiedad importante de las series temporales porque las técnicas de modelado de series temporales comunes, como ARIMA, suponen que los datos son estacionarios.

Hipótesis de la Prueba

Hipótesis Nula (\(H_0\)): La serie tiene una raíz unitaria (la serie no es estacionaria).
Hipótesis Alternativa (\(H_a\)): La serie no tiene una raíz unitaria (la serie es estacionaria).

Proceso de la Prueba

Preparación de los Datos:
- Se toma la serie temporal que se quiere analizar.
Ejecución del Test:
- Se aplica la prueba ADF a los datos de la serie temporal.
- La prueba calcula un valor estadístico de Dickey-Fuller y un valor p (p-value).
Interpretación de los Resultados:
- Si el valor p es menor que el nivel de significancia (por ejemplo, 0.05):
  - Rechazamos la hipótesis nula (\(H_0\)).
  - Concluimos que la serie es estacionaria.
- Si el valor p es mayor que el nivel de significancia:
  - No podemos rechazar la hipótesis nula (\(H_0\)).
  - Concluimos que la serie no es estacionaria.

Sin embargo, aunque los resultados indiquen estacionariedad en algunos casos, es importante entender que la prueba ADF evalúa la existencia de raíces unitarias en el modelo, lo que implica que los cambios en la serie de tiempo tienen efectos duraderos y la serie no retorna a su nivel promedio. Es posible que los datos tengan una tendencia subyacente, lo cual también debe ser considerado al analizar los resultados.

# Aplicamos la Prueba de Dickey-Fuller Aumentada a la serie white noise
adf.test(whitenoise)

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

## 
##  Augmented Dickey-Fuller Test
## 
## data:  whitenoise
## Dickey-Fuller = -10.351, Lag order = 9, p-value = 0.01
## alternative hypothesis: stationary

# Aplicamos la Prueba de Dickey-Fuller Aumentada a las series random walk con y sin drift 
adf.test(rwalk)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  rwalk
## Dickey-Fuller = -2.6209, Lag order = 9, p-value = 0.3155
## alternative hypothesis: stationary

adf.test(rwalk_drift)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  rwalk_drift
## Dickey-Fuller = -2.2442, Lag order = 9, p-value = 0.4749
## alternative hypothesis: stationary

adf.test(AirPassengers)

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

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

adfTest(AirPassengers, type = "nc")

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -0.3524
##   P VALUE:
##     0.5017 
## 
## Description:
##  Thu Aug  8 18:42:41 2024 by user: franc

adfTest(AirPassengers, type = "c")

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -2.345
##   P VALUE:
##     0.1861 
## 
## Description:
##  Thu Aug  8 18:42:41 2024 by user: franc

adfTest(AirPassengers, type = "ct")

## Warning in adfTest(AirPassengers, type = "ct"): p-value smaller than printed
## p-value

## 
## Title:
##  Augmented Dickey-Fuller Test
## 
## Test Results:
##   PARAMETER:
##     Lag Order: 1
##   STATISTIC:
##     Dickey-Fuller: -7.6523
##   P VALUE:
##     0.01 
## 
## Description:
##  Thu Aug  8 18:42:41 2024 by user: franc

Diferenciación

La diferenciación es una técnica utilizada en el análisis de series temporales para hacer que una serie no estacionaria se vuelva estacionaria. Esto se logra calculando las diferencias entre valores consecutivos de la serie temporal.

Proceso de Diferenciación

Serie Original:
- Comienza con la serie temporal original que deseas analizar. Observa si la serie presenta tendencias o patrones que indiquen que no es estacionaria.
Aplicación de la Diferenciación:
- Calcula las diferencias entre valores consecutivos de la serie temporal. La fórmula para la diferenciación de primer orden es \(\Delta y_t = y_t - y_{t-1}\).
Serie Diferenciada:
- La serie resultante después de la diferenciación se considera la serie diferenciada. La diferenciación debería eliminar tendencias y hacer que la serie sea más estacionaria.
Análisis de la Función de Autocorrelación (ACF):
- Calcula la función de autocorrelación (ACF) de la serie original y de la serie diferenciada.
- En la serie original, la ACF generalmente disminuye lentamente, lo que indica una dependencia a largo plazo y no estacionariedad.
- En la serie diferenciada, la ACF debería disminuir rápidamente, lo que es indicativo de una serie más estacionaria.

Ventajas de la Diferenciación

Eliminación de Tendencias: La diferenciación ayuda a eliminar las tendencias en los datos, haciendo que la serie sea estacionaria.
Facilita el Modelado: Al hacer que la serie sea estacionaria, se facilita la aplicación de modelos de series temporales como ARIMA.

Consideraciones

Orden de Diferenciación: A veces puede ser necesario aplicar diferenciación de segundo orden (diferenciar la serie ya diferenciada) si la serie sigue sin ser estacionaria después de la primera diferenciación.
Sobre-Diferenciación: Diferenciar excesivamente puede introducir ruido y hacer que la serie pierda su estructura significativa.

# Calculamos la primera diferenciación de la random walk
delta_rwalk = diff(rwalk)

# Graficamos la serie random walk original 
plot(rwalk,type = "l", main = "Random Walk Original")

# Graficamos la serie random walk diferenciada
plot(delta_rwalk,type = "l",main = "Primera diferenciación")

# Graficamos el correlograma de la serie original
acf(rwalk, main = "Correlograma de Autocorrelación")

# Graficamos el correlograma de la serie diferenciada
acf(delta_rwalk, main = "Correlograma de Autocorrelación")

Modelo ARIMA

El modelo ARIMA (AutoRegressive Integrated Moving Average) es una técnica de modelado estadístico usada principalmente para el análisis y pronóstico de series temporales. ARIMA combina tres componentes clave:

AR (AutoRegresivo): Utiliza la relación entre una observación y una serie de valores anteriores (retrasos) de la misma serie temporal. Se expresa con el término \(p\), que indica el número de términos autoregresivos.
I (Integrado): Se refiere a la diferenciación de las observaciones para hacer que la serie temporal sea estacionaria, es decir, que sus propiedades estadísticas, como la media y la varianza, no cambien con el tiempo. El término \(d\) indica el número de veces que los datos han sido diferenciados.
MA (Media Móvil): Utiliza la dependencia entre una observación y un error residual de un modelo de media móvil aplicado a retrasos anteriores. El término \(q\) indica el número de términos de media móvil.

El modelo ARIMA se denota generalmente como ARIMA(\(p, d, q\)).

La metodología para construir un modelo ARIMA típicamente incluye los siguientes pasos:

1. Identificación: Determinar si los datos son estacionarios y diferenciar la serie si es necesario. También identificar los posibles valores de \(p\) y \(q\) usando la función de autocorrelación (ACF) y la función de autocorrelación parcial (PACF).

2. Estimación: Ajustar el modelo a los datos mediante la estimación de los parámetros de \(p\), \(d\) y \(q\).

3. Diagnóstico: Evaluar el modelo ajustado revisando los residuos para asegurar que se comportan como ruido blanco (sin estructura aparente).

4. Pronóstico: Utilizar el modelo ajustado para predecir valores futuros de la serie temporal.

Modelo SARIMA

El modelo SARIMA (Seasonal Autoregressive Integrated Moving Average) es una extensión del modelo ARIMA que incorpora componentes estacionales para capturar patrones repetitivos a lo largo del tiempo. Este modelo es particularmente útil cuando los datos muestran patrones que se repiten a intervalos regulares, como días, semanas, meses o años.

Componentes del Modelo SARIMA

Un modelo SARIMA se denota como SARIMA(\(p,d,q,P,Q,D\))\([s]\), donde:

\(p\): Orden del componente autoregresivo (AR).
\(d\): Grado de diferenciación no estacional necesario para hacer estacionaria la serie.
\(q\): Orden del componente de media móvil (MA).
\(P\): Orden del componente autoregresivo estacional (SAR).
\(D\): Grado de diferenciación estacional necesario para hacer estacionaria la serie.
\(Q\): Orden del componente de media móvil estacional (SMA).
\(s\): Periodo de estacionalidad, es decir, el número de observaciones en un ciclo completo de la estacionalidad (por ejemplo, 12 para datos mensuales con un ciclo anual).

Uso de `auto.arima` en R

La función auto.arima en R es una herramienta del paquete forecast que automatiza el proceso de identificación, estimación y diagnóstico de un modelo ARIMA para una serie temporal. Selecciona el mejor modelo ARIMA basado en un criterio de información (como AIC, AICc o BIC).

Funcionamiento de `auto.arima`

Determinación de Estacionariedad: Determina si la serie temporal es estacionaria y aplica las diferencias necesarias para hacerla estacionaria, encontrando el valor adecuado de \(d\).
Selección de Modelos ARIMA: Explora diferentes combinaciones de \(p\) y \(q\) para modelos AR y MA.
Evaluación de Modelos: Para cada combinación de parámetros, estima el modelo ARIMA y calcula los criterios de información.
Selección del Mejor Modelo: Selecciona el modelo que minimiza el criterio de información.
Validación del Modelo: Revisa los residuos del modelo seleccionado para asegurar que se comportan como ruido blanco.

Parámetros Adicionales

auto.arima también permite ajustar varios parámetros para personalizar el proceso de selección del modelo:

seasonal: Especifica si se debe considerar un componente estacional.
stepwise: Determina si se utiliza un método de búsqueda escalonado (por defecto es TRUE, lo que acelera el proceso).
approximation: Indica si se deben utilizar aproximaciones para acelerar el cálculo (útil para series muy grandes).

Ejemplo:

Pronóstico de demanda mensual de azúcar

Fuente de datos: https://www.kaggle.com/code/sharmeens/supply-chain-analysis-and-forecasting-time-series/input

data = read.csv("all_agricultural_products_data.csv")
head(data)

##   ticker commodity       date open high low close volume
## 1   CC=F     Cocoa 2000-01-03  840  846 820   830   2426
## 2   CC=F     Cocoa 2000-01-04  830  841 823   836   1957
## 3   CC=F     Cocoa 2000-01-05  840  850 828   831   3975
## 4   CC=F     Cocoa 2000-01-06  830  847 824   841   3454
## 5   CC=F     Cocoa 2000-01-07  848  855 836   853   5008
## 6   CC=F     Cocoa 2000-01-10  848  855 837   839   3900

# Convertir la columna date a tipo Date
data$date <- as.Date(data$date, format="%Y-%m-%d")

# Crear una nueva columna year_month con el año y el mes
data <- data %>%
  mutate(year_month = format(date, "%Y-%m"))

# Verificar las primeras filas del dataframe
head(data)

##   ticker commodity       date open high low close volume year_month
## 1   CC=F     Cocoa 2000-01-03  840  846 820   830   2426    2000-01
## 2   CC=F     Cocoa 2000-01-04  830  841 823   836   1957    2000-01
## 3   CC=F     Cocoa 2000-01-05  840  850 828   831   3975    2000-01
## 4   CC=F     Cocoa 2000-01-06  830  847 824   841   3454    2000-01
## 5   CC=F     Cocoa 2000-01-07  848  855 836   853   5008    2000-01
## 6   CC=F     Cocoa 2000-01-10  848  855 837   839   3900    2000-01

# Agrupar por year_month y calcular las métricas deseadas
monthly_data <- data %>%
  group_by(year_month,commodity) %>%
  summarize(
    open = mean(open, na.rm=TRUE),
    high = mean(high, na.rm=TRUE),
    low = mean(low, na.rm=TRUE),
    close = mean(close, na.rm=TRUE),
    volume = sum(volume, na.rm=TRUE)
  )

## `summarise()` has grouped output by 'year_month'. You can override using the
## `.groups` argument.

# Verificar las primeras filas del dataframe agrupado
head(monthly_data)

## # A tibble: 6 × 7
## # Groups:   year_month [2]
##   year_month commodity  open  high   low close volume
##   <chr>      <chr>     <dbl> <dbl> <dbl> <dbl>  <int>
## 1 2000-01    Cocoa     842.  851.  830.  840.   88465
## 2 2000-01    Coffee    116.  118.  114.  115.  120093
## 3 2000-01    Cotton     54.7  55.5  54.2  55.1 180359
## 4 2000-02    Cocoa     766.  777.  755.  764.   73201
## 5 2000-02    Coffee    107.  108.  106.  106.   98735
## 6 2000-02    Cotton     56.7  57.5  56.1  56.8 108820

sugar_volume = monthly_data$volume[monthly_data$commodity == "Sugar"]

sugar_ts = ts(sugar_volume, frequency = 12, start = c(2000,1))

plot(sugar_ts,
     type ="l",
     ylab = "Volume",
     main = "Sugar")

# Recortar la serie temporal para que empiece en 2010
start_year <- 2010
sugar_ts_2010 <- window(sugar_ts, start=c(start_year, 1))

# Graficar la nueva serie temporal
plot(sugar_ts_2010,
     type="l",
     ylab="Volume",
     main="Sugar desde 2010")

plot(decompose(sugar_ts_2010))

model = auto.arima(sugar_ts_2010)
model

## Series: sugar_ts_2010 
## ARIMA(1,1,1)(2,0,0)[12] 
## 
## Coefficients:
##           ar1      ma1    sar1    sar2
##       -0.0837  -0.6460  0.2477  0.3240
## s.e.   0.0824   0.0749  0.0404  0.0544
## 
## sigma^2 = 3.666e+10:  log likelihood = -2241.42
## AIC=4492.84   AICc=4493.22   BIC=4508.37

residuals <- residuals(model)
Box.test(residuals, lag = 20, type = "Ljung-Box")

## 
##  Box-Ljung test
## 
## data:  residuals
## X-squared = 19.354, df = 20, p-value = 0.4989

pronostico <- forecast(model, h = 30)
plot(pronostico)
lines(sugar_ts_2010, type = "l")

pronostico

##          Point Forecast     Lo 80   Hi 80    Lo 95   Hi 95
## Nov 2023        1277900 1032537.7 1523263 902650.7 1653150
## Dec 2023        1288385 1034218.7 1542552 899671.0 1677100
## Jan 2024        1375438 1108585.7 1642290 967322.8 1783553
## Feb 2024        1285614 1007015.6 1564212 859534.5 1711694
## Mar 2024        1300858 1010959.9 1590756 857497.1 1744219
## Apr 2024        1401369 1100597.9 1702141 941379.1 1861360
## May 2024        1263282  952016.8 1574547 787242.9 1739321
## Jun 2024        1302633  981216.0 1624049 811068.3 1794197
## Jul 2024        1387193 1055936.4 1718451 880579.5 1893807
## Aug 2024        1302834  962020.5 1643648 781604.7 1824063
## Sep 2024        1388460 1038350.5 1738569 853013.9 1923905
## Oct 2024        1282111  922947.1 1641276 732817.0 1831406
## Nov 2024        1355037  969213.8 1740860 764971.4 1945102
## Dec 2024        1328786  931061.1 1726511 720518.2 1937054
## Jan 2025        1409959  999786.5 1820132 782654.2 2037264
## Feb 2025        1310665  888484.6 1732844 664996.1 1956333
## Mar 2025        1352729  918868.4 1786590 689196.3 2016262
## Apr 2025        1425818  980582.7 1871053 744889.4 2106746
## May 2025        1237093  780767.6 1693419 539203.2 1934984
## Jun 2025        1355825  888671.6 1822979 641375.5 2070275
## Jul 2025        1395171  917435.7 1872907 664537.7 2125805
## Aug 2025        1307191  819102.1 1795279 560723.6 2053658
## Sep 2025        1329262  831035.8 1827488 567290.7 2091233
## Oct 2025        1332990  824828.3 1841152 555823.6 2110157
## Nov 2025        1353927  809793.8 1898060 521747.0 2186107
## Dec 2025        1350821  792390.7 1909252 496775.3 2204867
## Jan 2026        1399138  825473.0 1972803 521793.1 2276483
## Feb 2026        1345434  757036.3 1933832 445557.4 2245310
## Mar 2026        1360794  758015.4 1963573 438923.5 2282665
## Apr 2026        1411469  794644.9 2028294 468117.7 2354821

Temas de Investigación de Operaciones

Introducción al Análisis de Series de Tiempo

Johnny F. Linares Martínez

2024-08-08

Componentes de una Serie de Tiempo

Ruido Blanco

Propiedades del Ruido Blanco

Random Walk

Propiedades del Random Walk

Prueba de Dickey-Fuller Aumentada

Hipótesis de la Prueba

Proceso de la Prueba

Diferenciación

Proceso de Diferenciación

Ventajas de la Diferenciación

Consideraciones

Modelo ARIMA

Modelo SARIMA

Componentes del Modelo SARIMA

Uso de `auto.arima` en R

Funcionamiento de `auto.arima`

Parámetros Adicionales

Ejemplo:

Temas de Investigación de Operaciones

Introducción al Análisis de Series de Tiempo

Johnny F. Linares Martínez

2024-08-08

Componentes de una Serie de Tiempo

Ruido Blanco

Propiedades del Ruido Blanco

Random Walk

Propiedades del Random Walk

Prueba de Dickey-Fuller Aumentada

Hipótesis de la Prueba

Proceso de la Prueba

Diferenciación

Proceso de Diferenciación

Ventajas de la Diferenciación

Consideraciones

Modelo ARIMA

Modelo SARIMA

Componentes del Modelo SARIMA

Uso de auto.arima en R

Funcionamiento de auto.arima

Parámetros Adicionales

Ejemplo:

Uso de `auto.arima` en R

Funcionamiento de `auto.arima`