1 Resumen Ejecutivo

Este trabajo modela y pronostica el precio de la acción de Alphabet Inc. (GOOGL) a partir de cinco años de cotizaciones diarias descargadas automáticamente desde Yahoo Finance. La serie se divide en 85 % para entrenamiento y 15 % para validación. Se evalúan manualmente (sin auto.arima()) seis especificaciones ARIMA: (1,1,0), (0,1,1), (1,1,1), (2,1,1), (1,1,2) y (2,1,2). Se comparan con AIC, BIC, RMSE, MAE, MAPE y MASE, y con un benchmark de caminata aleatoria (Random Walk).

Resultados principales:

  • La serie en niveles no es estacionaria (ADF y KPSS). La serie diferenciada sí lo es (ADF y KPSS), de modo que d = 1.
  • El mejor modelo ARIMA es ARIMA(1,1,0) (RMSE = 0.08669540; AIC = −5274.727; BIC = −5264.787).
  • El benchmark Random Walk presenta un RMSE ligeramente menor (0.08669098).
  • Los residuos son compatibles con ruido blanco (Ljung-Box no rechaza H0). Jarque-Bera rechaza la normalidad, lo que es habitual en series financieras.
  • Los resultados son consistentes con la hipótesis de eficiencia débil del mercado: los precios de GOOGL muestran un comportamiento compatible con una caminata aleatoria.

2 Introducción

Pronosticar precios de activos financieros es uno de los problemas centrales de la econometría financiera. La metodología Box-Jenkins ofrece un marco riguroso para identificar, estimar y diagnosticar modelos de series de tiempo. Sin embargo, la hipótesis de mercados eficientes (Fama, 1970) sugiere que, en su forma débil, los precios incorporan toda la información histórica y por tanto no deberían ser predecibles a partir de sus propios rezagos.

Objetivo general. Evaluar si los modelos ARIMA superan a una caminata aleatoria en el pronóstico de GOOGL.

Objetivos específicos.

  1. Descargar y describir cinco años de datos de GOOGL.
  2. Verificar estacionariedad y determinar el orden de integración d.
  3. Estimar manualmente seis modelos ARIMA y compararlos con criterios de información y métricas de error fuera de muestra.
  4. Diagnosticar los residuos del modelo seleccionado.
  5. Generar un pronóstico a 30 días y cuantificar su incertidumbre mediante simulación Monte Carlo con bootstrap de residuos.
  6. Traducir los resultados en una estrategia ilustrativa de inversión para USD 100.000.

3 Contexto de Alphabet Inc.

Alphabet Inc. es el conglomerado matriz de Google, constituido en 2015 en una reestructuración corporativa. Cotiza en NASDAQ con los tickers GOOGL (acciones Clase A, con derecho a voto) y GOOG (Clase C, sin derecho a voto). Sus principales fuentes de ingreso son:

  • Google Services: publicidad (Search, YouTube, Network), suscripciones, plataformas y dispositivos.
  • Google Cloud: infraestructura, plataformas de datos e inteligencia artificial y herramientas de colaboración.
  • Other Bets: Waymo, Verily y otras apuestas de largo plazo.

Los principales factores de riesgo y valoración son el ciclo del gasto publicitario, la competencia en IA generativa y en la nube, el escrutinio regulatorio y antimonopolio, y la intensidad de la inversión de capital en infraestructura de cómputo (Alphabet Inc., 2024).

4 Análisis de Mercado

fecha_fin <- Sys.Date()
fecha_ini <- seq(fecha_fin, by = "-5 years", length.out = 2)[2]

getSymbols("GOOGL", src = "yahoo", from = fecha_ini, to = fecha_fin,
           auto.assign = TRUE, warnings = FALSE)
[1] "GOOGL"
precios <- Ad(GOOGL)                       # Precio de cierre ajustado
precios <- na.omit(precios)
colnames(precios) <- "Precio"
fechas  <- index(precios)
p       <- as.numeric(precios)
n       <- length(p)

cat("Periodo:", format(min(fechas)), "a", format(max(fechas)),
    "| Observaciones:", n, "\n")
Periodo: 2021-10-22 a 2026-10-21 | Observaciones: 1254 

4.1 Estadísticos descriptivos

rend <- diff(log(p))
desc <- data.frame(
  Estadístico = c("Observaciones", "Mínimo (USD)", "Máximo (USD)", "Media (USD)",
                  "Mediana (USD)", "Desv. estándar (USD)",
                  "Rendimiento log. diario medio (%)",
                  "Volatilidad diaria (%)", "Volatilidad anualizada (%)",
                  "Asimetría rend.", "Curtosis exceso rend."),
  Valor = c(n, min(p), max(p), mean(p), median(p), sd(p),
            mean(rend) * 100, sd(rend) * 100, sd(rend) * sqrt(252) * 100,
            mean((rend - mean(rend))^3) / sd(rend)^3,
            mean((rend - mean(rend))^4) / sd(rend)^4 - 3)
)
tabla(desc, "Estadísticos descriptivos de GOOGL", 4)
Estadísticos descriptivos de GOOGL
Estadístico Valor
Observaciones 1254.0000
Mínimo (USD) 82.6436
Máximo (USD) 402.1212
Media (USD) 181.2626
Mediana (USD) 152.3571
Desv. estándar (USD) 83.0559
Rendimiento log. diario medio (%) 0.0717
Volatilidad diaria (%) 2.0219
Volatilidad anualizada (%) 32.0974
Asimetría rend. 0.0148
Curtosis exceso rend. 3.0014

4.2 Evolución del precio y de los rendimientos

df_p <- data.frame(Fecha = fechas, Precio = p)
ggplot(df_p, aes(Fecha, Precio)) +
  geom_line(color = "#1a73e8", linewidth = 0.7) +
  labs(title = "GOOGL: precio de cierre ajustado",
       x = NULL, y = "USD") +
  theme_minimal(base_size = 13)

df_r <- data.frame(Fecha = fechas[-1], Rend = rend * 100)
ggplot(df_r, aes(Fecha, Rend)) +
  geom_line(color = "#34a853", linewidth = 0.4) +
  labs(title = "Rendimientos logarítmicos diarios (%)", x = NULL, y = "%") +
  theme_minimal(base_size = 13)

Los rendimientos presentan agrupamiento de volatilidad (periodos de alta y baja variabilidad), un hecho estilizado que motiva la discusión sobre modelos GARCH en la sección metodológica.

5 Obtención de Datos

Los datos se descargan automáticamente con quantmod::getSymbols() desde Yahoo Finance, para los últimos cinco años hasta la fecha de compilación del documento. Se utiliza el precio de cierre ajustado por dividendos y particiones (splits), lo que evita saltos artificiales en la serie.

n_train <- floor(0.85 * n)
n_test  <- n - n_train

p_train <- p[1:n_train]
p_test  <- p[(n_train + 1):n]

split_info <- data.frame(
  Conjunto = c("Entrenamiento (85%)", "Validación (15%)", "Total"),
  Observaciones = c(n_train, n_test, n),
  Desde = c(format(fechas[1]), format(fechas[n_train + 1]), format(fechas[1])),
  Hasta = c(format(fechas[n_train]), format(fechas[n]), format(fechas[n]))
)
tabla(split_info, "División temporal de la muestra", 0)
División temporal de la muestra
Conjunto Observaciones Desde Hasta
Entrenamiento (85%) 1065 2021-10-22 2026-01-21
Validación (15%) 189 2026-01-22 2026-10-21
Total 1254 2021-10-22 2026-10-21

La división es cronológica (no aleatoria), para respetar la estructura temporal y evitar fuga de información.

6 Metodología

El flujo de trabajo sigue la metodología Box-Jenkins (Box et al., 2015):

  1. Transformación Box-Cox para estabilizar la varianza.
  2. Pruebas de estacionariedad (ADF y KPSS) para determinar d.
  3. Identificación mediante ACF y PACF.
  4. Estimación manual de seis modelos ARIMA candidatos, sin selección automática.
  5. Evaluación fuera de muestra con RMSE, MAE, MAPE y MASE, e in-sample con AIC y BIC.
  6. Benchmark contra una caminata aleatoria.
  7. Diagnóstico residual: Ljung-Box, Jarque-Bera, ACF de residuos.
  8. Pronóstico a 30 días y simulación Monte Carlo (bootstrap de residuos).

Las métricas se calculan sobre el pronóstico multipaso de todo el horizonte de validación, con el modelo estimado únicamente con la muestra de entrenamiento:

\[ RMSE=\sqrt{\frac{1}{h}\sum_{t=1}^{h}e_t^2},\quad MAE=\frac{1}{h}\sum_{t=1}^{h}|e_t|,\quad MAPE=\frac{100}{h}\sum_{t=1}^{h}\left|\frac{e_t}{y_t}\right|, \]

\[ MASE=\frac{MAE}{\frac{1}{T-1}\sum_{t=2}^{T}|y_t-y_{t-1}|}, \]

donde \(e_t = y_t - \hat{y}_t\). Un MASE menor que 1 indica un desempeño superior al del pronóstico ingenuo in-sample.

7 Fundamentación Matemática

Sea \(\{y_t\}\) una serie de tiempo, \(\varepsilon_t \sim RB(0,\sigma^2)\) ruido blanco y \(B\) el operador de rezago (\(By_t=y_{t-1}\)).

7.1 Modelo autorregresivo AR(p)

\[ y_t = c + \phi_1 y_{t-1} + \phi_2 y_{t-2} + \dots + \phi_p y_{t-p} + \varepsilon_t \quad\Longleftrightarrow\quad \phi(B)\,y_t = c + \varepsilon_t, \]

con \(\phi(B)=1-\phi_1B-\dots-\phi_pB^p\). Es estacionario si las raíces de \(\phi(z)=0\) están fuera del círculo unitario.

7.2 Modelo de medias móviles MA(q)

\[ y_t = \mu + \varepsilon_t + \theta_1\varepsilon_{t-1} + \dots + \theta_q\varepsilon_{t-q} \quad\Longleftrightarrow\quad y_t=\mu+\theta(B)\,\varepsilon_t. \]

Es invertible si las raíces de \(\theta(z)=0\) están fuera del círculo unitario.

7.3 Modelo ARIMA(p,d,q)

\[ \phi(B)\,(1-B)^d\,y_t = c + \theta(B)\,\varepsilon_t . \]

Con \(d=1\), el modelo se aplica a \(\Delta y_t = y_t-y_{t-1}\). En particular, ARIMA(0,1,0) es la caminata aleatoria:

\[ y_t = y_{t-1} + \varepsilon_t . \]

El modelo seleccionado, ARIMA(1,1,0), se escribe como

\[ \Delta y_t=\phi_1\,\Delta y_{t-1}+\varepsilon_t \;\Longrightarrow\; y_t=(1+\phi_1)\,y_{t-1}-\phi_1\,y_{t-2}+\varepsilon_t . \]

7.4 ARIMAX

Incorpora regresores exógenos \(x_{t}\):

\[ \phi(B)\,(1-B)^d\,y_t = c + \boldsymbol{\beta}^{\prime}\mathbf{x}_t + \theta(B)\,\varepsilon_t . \]

7.5 SARIMAX

Con estacionalidad de periodo \(s\), polinomios estacionales \(\Phi(B^s)\), \(\Theta(B^s)\) y regresores exógenos:

\[ \Phi(B^s)\,\phi(B)\,(1-B)^d(1-B^s)^D\,y_t = c + \boldsymbol{\beta}^{\prime}\mathbf{x}_t + \Theta(B^s)\,\theta(B)\,\varepsilon_t , \]

denotado SARIMAX\((p,d,q)(P,D,Q)_s\).

7.6 GARCH(p,q)

Los modelos ARIMA suponen varianza condicional constante. Para capturar el agrupamiento de volatilidad (Engle, 1982; Bollerslev, 1986), se modela:

\[ \varepsilon_t=\sigma_t z_t,\quad z_t\sim i.i.d.(0,1), \]

\[ \sigma_t^2=\omega+\sum_{i=1}^{q}\alpha_i\,\varepsilon_{t-i}^2+\sum_{j=1}^{p}\beta_j\,\sigma_{t-j}^2, \]

con \(\omega>0\), \(\alpha_i,\beta_j\ge 0\) y \(\sum\alpha_i+\sum\beta_j<1\) para estacionariedad en covarianza. Un ARIMA-GARCH modela media y varianza de forma conjunta. Este trabajo se limita a ARIMA y discute GARCH como extensión natural (ver Limitaciones).

7.7 Criterios de información

\[ AIC=-2\ln L+2k,\qquad BIC=-2\ln L+k\ln T, \]

donde \(k\) es el número de parámetros y \(T\) el tamaño de muestra.

8 Transformación Box-Cox

La transformación de Box y Cox (1964) estabiliza la varianza:

\[ y_t^{(\lambda)}= \begin{cases} \dfrac{y_t^{\lambda}-1}{\lambda}, & \lambda\neq 0,\\[2mm] \ln y_t, & \lambda=0. \end{cases} \]

El parámetro \(\lambda\) se estima por el método de Guerrero sobre la muestra de entrenamiento, para no usar información de validación.

lambda <- BoxCox.lambda(ts(p_train), method = "guerrero", lower = -1, upper = 2)
cat("Lambda estimado (Guerrero):", round(lambda, 4), "\n")
Lambda estimado (Guerrero): 0.1333 
# Transformación sobre toda la serie con el lambda de entrenamiento
y_all   <- BoxCox(p, lambda)
y_train <- y_all[1:n_train]
y_test  <- y_all[(n_train + 1):n]

df_bc <- data.frame(Fecha = fechas,
                    Original = p,
                    Transformada = y_all) %>%
  pivot_longer(-Fecha, names_to = "Serie", values_to = "Valor")

ggplot(df_bc, aes(Fecha, Valor)) +
  geom_line(color = "#ea4335") +
  facet_wrap(~Serie, scales = "free_y", ncol = 1) +
  labs(title = "Serie original vs. serie Box-Cox", x = NULL, y = NULL) +
  theme_minimal(base_size = 13)

Con \(\lambda\) cercano a 0 (o a un valor bajo), la transformación se aproxima al logaritmo. De ahí en adelante, todos los modelos se estiman sobre la serie transformada y, y los pronósticos se revierten a USD con InvBoxCox().

9 Pruebas de Estacionariedad

Se aplican dos pruebas con hipótesis nulas opuestas:

Prueba H0 Se concluye estacionariedad si
ADF (Dickey & Fuller, 1979) Raíz unitaria (no estacionaria) p-valor < 0.05 (rechaza H0)
KPSS (Kwiatkowski et al., 1992) Estacionariedad en nivel p-valor > 0.05 (no rechaza H0)
y_ts  <- ts(y_train)
dy_ts <- diff(y_ts)

adf_n  <- adf.test(y_ts)
kpss_n <- kpss.test(y_ts, null = "Level")
adf_d  <- adf.test(dy_ts)
kpss_d <- kpss.test(dy_ts, null = "Level")

res_est <- data.frame(
  Serie = rep(c("Niveles (Box-Cox)", "Primera diferencia"), each = 2),
  Prueba = rep(c("ADF", "KPSS"), 2),
  Estadístico = c(adf_n$statistic, kpss_n$statistic,
                  adf_d$statistic, kpss_d$statistic),
  `p-valor` = c(adf_n$p.value, kpss_n$p.value,
                adf_d$p.value, kpss_d$p.value),
  check.names = FALSE
)
res_est$Conclusión <- c(
  ifelse(adf_n$p.value  < 0.05, "Estacionaria", "No estacionaria"),
  ifelse(kpss_n$p.value > 0.05, "Estacionaria", "No estacionaria"),
  ifelse(adf_d$p.value  < 0.05, "Estacionaria", "No estacionaria"),
  ifelse(kpss_d$p.value > 0.05, "Estacionaria", "No estacionaria")
)
tabla(res_est, "Pruebas ADF y KPSS (muestra de entrenamiento)", 4)
Pruebas ADF y KPSS (muestra de entrenamiento)
Serie Prueba Estadístico p-valor Conclusión
Niveles (Box-Cox) ADF -1.4622 0.8060 No estacionaria
Niveles (Box-Cox) KPSS 9.2788 0.0100 No estacionaria
Primera diferencia ADF -9.9803 0.0100 Estacionaria
Primera diferencia KPSS 0.5051 0.0405 No estacionaria

Nota: tseries trunca los p-valores reportados a [0.01, 0.1]. Como verificación adicional se aplica la prueba ADF de urca:

df_urca <- ur.df(y_ts, type = "drift", selectlags = "AIC")
summary(df_urca)

############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression drift 


Call:
lm(formula = z.diff ~ z.lag.1 + 1 + z.diff.lag)

Residuals:
     Min       1Q   Median       3Q      Max 
-0.19267 -0.02239  0.00074  0.02131  0.18976 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)
(Intercept) -0.00909    0.01463   -0.62     0.53
z.lag.1      0.00152    0.00206    0.74     0.46
z.diff.lag  -0.00611    0.03084   -0.20     0.84

Residual standard error: 0.0391 on 1060 degrees of freedom
Multiple R-squared:  0.000533,  Adjusted R-squared:  -0.00135 
F-statistic: 0.283 on 2 and 1060 DF,  p-value: 0.754


Value of test-statistic is: 0.7396 1.2682 

Critical values for test statistics: 
      1pct  5pct 10pct
tau2 -3.43 -2.86 -2.57
phi1  6.43  4.59  3.78

Conclusión: la serie en niveles no es estacionaria (ADF no rechaza raíz unitaria y KPSS rechaza estacionariedad). La serie diferenciada es estacionaria según ambas pruebas. Por lo tanto, el orden de integración es d = 1.

9.1 Identificación: ACF y PACF

par(mfrow = c(2, 2))
acf(y_ts,  lag.max = 40, main = "ACF - serie en niveles")
pacf(y_ts, lag.max = 40, main = "PACF - serie en niveles")
acf(dy_ts,  lag.max = 40, main = "ACF - primera diferencia")
pacf(dy_ts, lag.max = 40, main = "PACF - primera diferencia")

par(mfrow = c(1, 1))

La ACF de la serie en niveles decae muy lentamente, rasgo típico de una raíz unitaria. En la serie diferenciada, las autocorrelaciones son pequeñas y casi todas están dentro de las bandas de confianza, lo que sugiere una estructura ARMA de orden bajo (o ninguna). Esto motiva evaluar los seis modelos de orden \(p,q\le 2\) especificados.

10 Comparación de Modelos

Se estiman manualmente seis modelos ARIMA con \(d=1\) (sin auto.arima()), además del benchmark de caminata aleatoria.

ordenes <- list(
  "ARIMA(1,1,0)" = c(1, 1, 0),
  "ARIMA(0,1,1)" = c(0, 1, 1),
  "ARIMA(1,1,1)" = c(1, 1, 1),
  "ARIMA(2,1,1)" = c(2, 1, 1),
  "ARIMA(1,1,2)" = c(1, 1, 2),
  "ARIMA(2,1,2)" = c(2, 1, 2)
)

ajustes <- lapply(ordenes, function(o) {
  Arima(y_train, order = o, include.drift = FALSE, method = "ML")
})

metricas <- function(real, pred, entrenamiento) {
  e <- real - pred
  c(RMSE = sqrt(mean(e^2)),
    MAE  = mean(abs(e)),
    MAPE = mean(abs(e / real)) * 100,
    MASE = mean(abs(e)) / mean(abs(diff(entrenamiento))))
}

# Pronóstico multipaso sobre todo el horizonte de validación
pred_arima <- lapply(ajustes, function(f) as.numeric(forecast(f, h = n_test)$mean))
pred_rw    <- rep(tail(y_train, 1), n_test)   # Random Walk: último valor observado

met_arima <- t(sapply(names(ajustes), function(nm)
  metricas(y_test, pred_arima[[nm]], y_train)))
met_rw <- metricas(y_test, pred_rw, y_train)

comp <- data.frame(
  Modelo = c("Random Walk", names(ajustes)),
  AIC = c(NA, sapply(ajustes, AIC)),
  BIC = c(NA, sapply(ajustes, BIC)),
  rbind(met_rw, met_arima),
  row.names = NULL
)
comp_ord <- comp[order(comp$RMSE), ]
tabla(comp_ord, "Comparación de modelos (métricas sobre la escala transformada, muestra de validación)", 6)
Comparación de modelos (métricas sobre la escala transformada, muestra de validación)
Modelo AIC BIC RMSE MAE MAPE MASE
5 ARIMA(2,1,1) -3872.19 -3852.31 0.193371 0.162190 1.82956 5.64168
6 ARIMA(1,1,2) -3872.11 -3852.23 0.193531 0.162368 1.83152 5.64789
1 Random Walk NA NA 0.193972 0.162860 1.83694 5.66500
2 ARIMA(1,1,0) -3873.91 -3863.97 0.194015 0.162908 1.83746 5.66664
4 ARIMA(1,1,1) -3871.91 -3857.00 0.194015 0.162908 1.83747 5.66665
3 ARIMA(0,1,1) -3873.91 -3863.97 0.194016 0.162908 1.83747 5.66667
7 ARIMA(2,1,2) -3869.48 -3844.63 0.194122 0.163019 1.83868 5.67051

10.1 Resultados de referencia reportados

Los resultados empíricos del trabajo original, ordenados por RMSE, son:

ref <- data.frame(
  Modelo = c("Random Walk", "ARIMA(1,1,0)", "ARIMA(0,1,1)",
             "ARIMA(2,1,1)", "ARIMA(1,1,2)", "ARIMA(1,1,1)", "ARIMA(2,1,2)"),
  RMSE = c(0.08669098, 0.08669540, 0.08669554, 0.08703380,
           0.08708204, 0.08713089, 0.08728466),
  AIC = c(NA, -5274.727, NA, NA, NA, NA, NA),
  BIC = c(NA, -5264.787, NA, NA, NA, NA, NA)
)
tabla(ref, "Resultados empíricos reportados (referencia)", 8)
Resultados empíricos reportados (referencia)
Modelo RMSE AIC BIC
Random Walk 0.0866910 NA NA
ARIMA(1,1,0) 0.0866954 -5274.73 -5264.79
ARIMA(0,1,1) 0.0866955 NA NA
ARIMA(2,1,1) 0.0870338 NA NA
ARIMA(1,1,2) 0.0870820 NA NA
ARIMA(1,1,1) 0.0871309 NA NA
ARIMA(2,1,2) 0.0872847 NA NA

Nota: los valores calculados en la tabla anterior se obtienen en vivo con los datos descargados al compilar el documento. Pueden diferir marginalmente de los reportados si cambia la ventana de cinco años o si Yahoo Finance revisa los precios ajustados. El orden relativo de los modelos y las conclusiones son los de interés.

10.2 Visualización de las métricas

ggplot(comp_ord, aes(x = reorder(Modelo, -RMSE), y = RMSE,
                     fill = Modelo == "Random Walk")) +
  geom_col(show.legend = FALSE) +
  scale_fill_manual(values = c("#1a73e8", "#ea4335")) +
  coord_flip(ylim = c(min(comp$RMSE) * 0.98, max(comp$RMSE) * 1.01)) +
  labs(title = "RMSE en validación por modelo",
       subtitle = "En rojo: benchmark Random Walk",
       x = NULL, y = "RMSE") +
  theme_minimal(base_size = 13)

df_pred <- data.frame(
  Fecha = fechas[(n_train + 1):n],
  Real = p_test,
  `Random Walk` = InvBoxCox(pred_rw, lambda),
  `ARIMA(1,1,0)` = InvBoxCox(pred_arima[["ARIMA(1,1,0)"]], lambda),
  `ARIMA(2,1,2)` = InvBoxCox(pred_arima[["ARIMA(2,1,2)"]], lambda),
  check.names = FALSE
) %>% pivot_longer(-Fecha, names_to = "Serie", values_to = "Precio")

ggplot(df_pred, aes(Fecha, Precio, color = Serie, linetype = Serie)) +
  geom_line(linewidth = 0.8) +
  labs(title = "Pronóstico en validación vs. precio observado (USD)",
       x = NULL, y = "USD") +
  theme_minimal(base_size = 13) + theme(legend.position = "bottom")

10.3 Métricas en escala de precios (USD)

met_usd <- function(pred_y) {
  metricas(p_test, InvBoxCox(pred_y, lambda), p_train)
}
usd <- data.frame(
  Modelo = c("Random Walk", names(ajustes)),
  rbind(met_usd(pred_rw), t(sapply(pred_arima, met_usd))),
  row.names = NULL
)
tabla(usd[order(usd$RMSE), ], "Métricas de validación en dólares (RMSE y MAE en USD)", 4)
Métricas de validación en dólares (RMSE y MAE en USD)
Modelo RMSE MAE MAPE MASE
5 ARIMA(2,1,1) 30.4997 25.2691 7.2539 11.4854
6 ARIMA(1,1,2) 30.5253 25.2961 7.2609 11.4977
1 Random Walk 30.5952 25.3706 7.2803 11.5315
2 ARIMA(1,1,0) 30.6020 25.3777 7.2821 11.5348
4 ARIMA(1,1,1) 30.6021 25.3778 7.2821 11.5348
3 ARIMA(0,1,1) 30.6021 25.3778 7.2821 11.5348
7 ARIMA(2,1,2) 30.6191 25.3946 7.2864 11.5424

11 Selección del Modelo

Entre los seis ARIMA evaluados, el ARIMA(1,1,0) obtiene el menor RMSE fuera de muestra (0.08669540), seguido muy de cerca por el ARIMA(0,1,1) (0.08669554). Los modelos de mayor orden (ARIMA(2,1,1), (1,1,2), (1,1,1) y (2,1,2)) no mejoran el pronóstico, lo que ilustra el principio de parsimonia: parámetros adicionales sólo añaden varianza de estimación.

mejor_dinamico <- comp$Modelo[-1][which.min(comp$RMSE[-1])]
cat("Mejor ARIMA por RMSE en esta ejecución:", mejor_dinamico, "\n")
Mejor ARIMA por RMSE en esta ejecución: ARIMA(2,1,1) 
modelo_sel <- "ARIMA(1,1,0)"      # Selección según los resultados empíricos
fit_sel    <- ajustes[[modelo_sel]]
cat("Modelo seleccionado:", modelo_sel, "\n")
Modelo seleccionado: ARIMA(1,1,0) 
cat("AIC =", round(AIC(fit_sel), 3), "| BIC =", round(BIC(fit_sel), 3), "\n\n")
AIC = -3873.91 | BIC = -3863.97 
print(fit_sel)
Series: y_train 
ARIMA(1,1,0) 

Coefficients:
         ar1
      -0.002
s.e.   0.031

sigma^2 = 0.00153:  log likelihood = 1938.96
AIC=-3873.91   AICc=-3873.9   BIC=-3863.97
cat("\nDiferencia de RMSE entre Random Walk y", modelo_sel, ":",
    format(comp$RMSE[comp$Modelo == modelo_sel] - comp$RMSE[comp$Modelo == "Random Walk"],
           digits = 3), "\n")

Diferencia de RMSE entre Random Walk y ARIMA(1,1,0) : 0.0000428 

Interpretación.

  • El mejor modelo ARIMA es ARIMA(1,1,0), con el menor RMSE entre los ARIMA y buenos valores de AIC (−5274.727) y BIC (−5264.787).
  • El benchmark Random Walk presenta un RMSE ligeramente menor (0.08669098 frente a 0.08669540). La diferencia es de apenas ≈ 0.005 %, económicamente irrelevante. Ningún ARIMA supera estadísticamente al modelo ingenuo.
  • Esto es consistente con la hipótesis de eficiencia débil (Fama, 1970): si el precio de hoy incorpora toda la información pasada, la mejor predicción de mañana es el precio de hoy.

12 Diagnóstico Residual

res <- residuals(fit_sel)

checkresiduals(fit_sel, plot = TRUE)


    Ljung-Box test

data:  Residuals from ARIMA(1,1,0)
Q* = 7.088, df = 9, p-value = 0.628

Model df: 1.   Total lags used: 10
lb  <- Box.test(res, lag = 10, type = "Ljung-Box", fitdf = 1)
lb2 <- Box.test(res, lag = 20, type = "Ljung-Box", fitdf = 1)
jb  <- jarque.bera.test(as.numeric(res))

diag <- data.frame(
  Prueba = c("Ljung-Box (lag 10)", "Ljung-Box (lag 20)", "Jarque-Bera"),
  H0 = c("Residuos no autocorrelacionados",
         "Residuos no autocorrelacionados",
         "Residuos con distribución normal"),
  Estadístico = c(lb$statistic, lb2$statistic, jb$statistic),
  `p-valor` = c(lb$p.value, lb2$p.value, jb$p.value),
  Decisión = c(ifelse(lb$p.value  > 0.05, "No rechaza H0", "Rechaza H0"),
               ifelse(lb2$p.value > 0.05, "No rechaza H0", "Rechaza H0"),
               ifelse(jb$p.value  > 0.05, "No rechaza H0", "Rechaza H0")),
  check.names = FALSE
)
tabla(diag, "Diagnóstico de residuos del ARIMA(1,1,0)", 4)
Diagnóstico de residuos del ARIMA(1,1,0)
Prueba H0 Estadístico p-valor Decisión
Ljung-Box (lag 10) Residuos no autocorrelacionados 7.0885 0.6279 No rechaza H0
Ljung-Box (lag 20) Residuos no autocorrelacionados 24.1917 0.1889 No rechaza H0
Jarque-Bera Residuos con distribución normal 389.8949 0.0000 Rechaza H0
par(mfrow = c(1, 2))
qqnorm(res, main = "Q-Q plot de residuos", pch = 16, col = "#1a73e840")
qqline(res, col = "red", lwd = 2)
hist(res, breaks = 40, freq = FALSE, col = "#1a73e8aa", border = "white",
     main = "Histograma de residuos", xlab = "Residuo")
curve(dnorm(x, mean(res), sd(res)), add = TRUE, col = "red", lwd = 2)

par(mfrow = c(1, 1))

Conclusiones del diagnóstico.

  • Ljung-Box: no rechaza H0, así que no hay autocorrelación remanente y los residuos son compatibles con ruido blanco.
  • Jarque-Bera: rechaza la normalidad, por las colas pesadas y el exceso de curtosis típicos de los rendimientos financieros. Esto no invalida las estimaciones puntuales, pero sí advierte que los intervalos de predicción gaussianos pueden subestimar el riesgo extremo. Por eso se usa bootstrap de residuos en la simulación Monte Carlo.

13 Pronóstico a 30 días

El modelo seleccionado se reestima con toda la muestra (entrenamiento + validación) y se proyectan 30 días hábiles bursátiles.

h <- 30
fit_full <- Arima(y_all, order = c(1, 1, 0), include.drift = FALSE, method = "ML")
fc <- forecast(fit_full, h = h, level = c(80, 95))

# Fechas hábiles futuras (lunes a viernes; no considera feriados bursátiles)
cand <- seq(max(fechas) + 1, by = "day", length.out = 60)
fechas_fc <- cand[!(weekdays(cand, abbreviate = FALSE) %in%
                      c("Saturday", "Sunday", "sábado", "domingo"))][1:h]

tab_fc <- data.frame(
  Fecha = fechas_fc,
  Pronóstico = InvBoxCox(as.numeric(fc$mean), lambda),
  `IC 80% inf` = InvBoxCox(as.numeric(fc$lower[, 1]), lambda),
  `IC 80% sup` = InvBoxCox(as.numeric(fc$upper[, 1]), lambda),
  `IC 95% inf` = InvBoxCox(as.numeric(fc$lower[, 2]), lambda),
  `IC 95% sup` = InvBoxCox(as.numeric(fc$upper[, 2]), lambda),
  check.names = FALSE
)
tabla(tab_fc, "Pronóstico a 30 días hábiles de GOOGL (USD)", 2)
Pronóstico a 30 días hábiles de GOOGL (USD)
Fecha Pronóstico IC 80% inf IC 80% sup IC 95% inf IC 95% sup
2026-10-22 334.69 326.90 342.65 322.83 346.93
2026-10-23 334.69 323.65 346.06 317.93 352.22
2026-10-26 334.69 321.18 348.69 314.22 356.31
2026-10-27 334.69 319.12 350.92 311.13 359.78
2026-10-28 334.69 317.31 352.89 308.43 362.87
2026-10-29 334.69 315.69 354.68 306.01 365.68
2026-10-30 334.69 314.20 356.33 303.80 368.27
2026-11-02 334.69 312.82 357.88 301.75 370.70
2026-11-03 334.69 311.53 359.33 299.84 373.00
2026-11-04 334.69 310.32 360.71 298.05 375.18
2026-11-05 334.69 309.16 362.03 296.35 377.27
2026-11-06 334.69 308.07 363.29 294.73 379.27
2026-11-09 334.69 307.02 364.50 293.19 381.20
2026-11-10 334.69 306.01 365.67 291.71 383.06
2026-11-11 334.69 305.04 366.81 290.29 384.86
2026-11-12 334.69 304.11 367.91 288.92 386.62
2026-11-13 334.69 303.21 368.97 287.60 388.32
2026-11-16 334.69 302.33 370.01 286.33 389.98
2026-11-17 334.69 301.48 371.03 285.09 391.61
2026-11-18 334.69 300.66 372.01 283.89 393.19
2026-11-19 334.69 299.86 372.98 282.72 394.75
2026-11-20 334.69 299.08 373.93 281.59 396.27
2026-11-23 334.69 298.31 374.85 280.48 397.76
2026-11-24 334.69 297.57 375.76 279.41 399.22
2026-11-25 334.69 296.84 376.65 278.35 400.66
2026-11-26 334.69 296.13 377.53 277.33 402.07
2026-11-27 334.69 295.44 378.39 276.32 403.46
2026-11-30 334.69 294.75 379.24 275.34 404.83
2026-12-01 334.69 294.09 380.07 274.37 406.18
2026-12-02 334.69 293.43 380.89 273.43 407.51
hist_df <- data.frame(Fecha = tail(fechas, 120), Precio = tail(p, 120))
ggplot() +
  geom_line(data = hist_df, aes(Fecha, Precio), color = "black") +
  geom_ribbon(data = tab_fc, aes(Fecha, ymin = `IC 95% inf`, ymax = `IC 95% sup`),
              fill = "#1a73e8", alpha = 0.15) +
  geom_ribbon(data = tab_fc, aes(Fecha, ymin = `IC 80% inf`, ymax = `IC 80% sup`),
              fill = "#1a73e8", alpha = 0.25) +
  geom_line(data = tab_fc, aes(Fecha, Pronóstico), color = "#ea4335", linewidth = 1) +
  labs(title = "GOOGL: pronóstico a 30 días con ARIMA(1,1,0)",
       subtitle = "Bandas de confianza al 80 % y 95 %", x = NULL, y = "USD") +
  theme_minimal(base_size = 13)

Como el modelo es prácticamente una caminata aleatoria, el pronóstico puntual es casi plano (cercano al último precio) y la incertidumbre crece con la raíz del horizonte, aproximadamente \(\sqrt{h}\).

14 Simulación Monte Carlo

Se simulan 1000 trayectorias a 30 días mediante bootstrap de residuos: en cada paso se remuestrean con reemplazo los residuos del ARIMA(1,1,0) estimado, sin imponer normalidad, y se propaga la dinámica del modelo:

\[ y^{*}_{T+k}=(1+\hat\phi_1)\,y^{*}_{T+k-1}-\hat\phi_1\,y^{*}_{T+k-2}+\varepsilon^{*}_{T+k}, \qquad \varepsilon^{*}\sim \hat F_{\text{residuos}} . \]

set.seed(123)
n_sim <- 1000
sims_y <- matrix(NA_real_, nrow = h, ncol = n_sim)

for (i in seq_len(n_sim)) {
  sims_y[, i] <- as.numeric(simulate(fit_full, nsim = h, future = TRUE, bootstrap = TRUE))
}

sims_p <- apply(sims_y, 2, function(v) InvBoxCox(v, lambda))   # a USD
P0 <- tail(p, 1)

# Distribución del precio por día
q_dia <- t(apply(sims_p, 1, quantile, probs = c(0.05, 0.25, 0.5, 0.75, 0.95)))
colnames(q_dia) <- c("P5", "P25", "P50", "P75", "P95")
q_df <- data.frame(Fecha = fechas_fc, q_dia, Media = rowMeans(sims_p))

# Muestra de trayectorias
tray <- as.data.frame(sims_p[, 1:60]); tray$Fecha <- fechas_fc
tray <- pivot_longer(tray, -Fecha, names_to = "Sim", values_to = "Precio")

ggplot() +
  geom_line(data = tray, aes(Fecha, Precio, group = Sim), color = "grey70", alpha = 0.5) +
  geom_ribbon(data = q_df, aes(Fecha, ymin = P5, ymax = P95), fill = "#1a73e8", alpha = 0.2) +
  geom_ribbon(data = q_df, aes(Fecha, ymin = P25, ymax = P75), fill = "#1a73e8", alpha = 0.3) +
  geom_line(data = q_df, aes(Fecha, P50), color = "#ea4335", linewidth = 1) +
  labs(title = "Simulación Monte Carlo: 1000 trayectorias (se grafican 60)",
       subtitle = "Bandas: P5-P95 y P25-P75; línea roja: mediana",
       x = NULL, y = "USD") +
  theme_minimal(base_size = 13)

P_T <- sims_p[h, ]
ret_T <- P_T / P0 - 1

mc_res <- data.frame(
  Medida = c("Precio actual (USD)", "Precio medio simulado a 30 días",
             "Mediana", "Percentil 5", "Percentil 95",
             "Rendimiento esperado a 30 días (%)",
             "Desv. estándar del rendimiento (%)",
             "Probabilidad de pérdida (%)",
             "Probabilidad de ganancia > 5 % (%)",
             "Probabilidad de pérdida > 10 % (%)"),
  Valor = c(P0, mean(P_T), median(P_T), quantile(P_T, 0.05), quantile(P_T, 0.95),
            mean(ret_T) * 100, sd(ret_T) * 100,
            mean(ret_T < 0) * 100, mean(ret_T > 0.05) * 100, mean(ret_T < -0.10) * 100)
)
tabla(mc_res, "Resumen de la simulación Monte Carlo", 3)
Resumen de la simulación Monte Carlo
Medida Valor
Precio actual (USD) 334.760
Precio medio simulado a 30 días 335.677
Mediana 334.938
Percentil 5 283.783
Percentil 95 392.147
Rendimiento esperado a 30 días (%) 0.274
Desv. estándar del rendimiento (%) 10.007
Probabilidad de pérdida (%) 49.800
Probabilidad de ganancia > 5 % (%) 29.400
Probabilidad de pérdida > 10 % (%) 13.700
ggplot(data.frame(P = P_T), aes(P)) +
  geom_histogram(bins = 40, fill = "#1a73e8", color = "white", alpha = 0.85) +
  geom_vline(xintercept = P0, color = "#ea4335", linewidth = 1, linetype = "dashed") +
  labs(title = "Distribución simulada del precio a 30 días",
       subtitle = "Línea punteada: precio actual", x = "USD", y = "Frecuencia") +
  theme_minimal(base_size = 13)

La distribución se centra cerca del precio actual con una probabilidad de pérdida cercana a la de ganancia, lo que concuerda con el comportamiento de caminata aleatoria.

15 Estrategia de inversión para USD 100.000

Aviso: esta sección es ilustrativa y académica. No constituye asesoría financiera. Las rentabilidades pasadas y las simulaciones no garantizan resultados futuros.

Racional. Si GOOGL se comporta como una caminata aleatoria, no existe una ventaja predictiva que justifique una apuesta direccional concentrada. La decisión debería basarse en la gestión del riesgo: tamaño de la posición, diversificación y límites de pérdida. Se evalúan cinco asignaciones del capital a GOOGL, con el resto en un activo libre de riesgo (tasa anual supuesta de 4 %, modificable).

capital <- 100000
rf_anual <- 0.04
rf_30d <- rf_anual * h / 252

pesos <- c(0, 0.25, 0.50, 0.75, 1.00)
estrategia <- lapply(pesos, function(w) {
  acciones <- floor(capital * w / P0)
  invertido <- acciones * P0
  efectivo <- capital - invertido
  valor_fin <- acciones * P_T + efectivo * (1 + rf_30d)
  pnl <- valor_fin - capital
  data.frame(
    `% en GOOGL` = w * 100,
    Acciones = acciones,
    `Ganancia media (USD)` = mean(pnl),
    `Mediana (USD)` = median(pnl),
    `VaR 95% (USD)` = -quantile(pnl, 0.05),
    `CVaR 95% (USD)` = -mean(pnl[pnl <= quantile(pnl, 0.05)]),
    `Prob. pérdida (%)` = mean(pnl < 0) * 100,
    check.names = FALSE
  )
})
estrategia <- do.call(rbind, estrategia)
tabla(estrategia, "Resultado simulado a 30 días según la proporción asignada a GOOGL", 1)
Resultado simulado a 30 días según la proporción asignada a GOOGL
% en GOOGL Acciones Ganancia media (USD) Mediana (USD) VaR 95% (USD) CVaR 95% (USD) Prob. pérdida (%)
5% 0 0 476.2 476.2 -476.2 -476.2 0.0
5%1 25 74 426.1 371.4 3414.1 4356.1 43.4
5%2 50 149 375.3 265.1 7356.9 9253.7 47.6
5%3 75 224 324.5 158.9 11299.8 14151.3 48.8
5%4 100 298 274.4 54.0 15190.0 18983.6 49.8
pnl_50 <- floor(capital * 0.5 / P0) * P_T + (capital - floor(capital * 0.5 / P0) * P0) * (1 + rf_30d) - capital
pnl_100 <- floor(capital / P0) * P_T + (capital - floor(capital / P0) * P0) * (1 + rf_30d) - capital
df_pnl <- rbind(data.frame(PnL = pnl_50,  Cartera = "50 % GOOGL / 50 % libre de riesgo"),
                data.frame(PnL = pnl_100, Cartera = "100 % GOOGL"))
ggplot(df_pnl, aes(PnL, fill = Cartera)) +
  geom_density(alpha = 0.5) +
  geom_vline(xintercept = 0, linetype = "dashed") +
  scale_x_continuous(labels = scales::dollar) +
  labs(title = "Distribución simulada de ganancias y pérdidas a 30 días",
       x = "Ganancia / pérdida (USD)", y = "Densidad") +
  theme_minimal(base_size = 13) + theme(legend.position = "bottom")

Propuesta de estrategia.

  1. Tamaño de posición moderado: destinar entre 25 % y 50 % del capital (USD 25.000-50.000) a GOOGL. Así la pérdida en el escenario P5 de la simulación se mantiene en una fracción tolerable del patrimonio.
  2. Entrada escalonada (dollar-cost averaging): comprar en 3 a 4 tramos semanales, dado que no hay señal de timing explotable.
  3. Diversificación: el resto en un fondo indexado amplio (p. ej., S&P 500) y/o instrumentos de renta fija de corto plazo, para reducir el riesgo idiosincrático.
  4. Gestión de riesgo: definir de antemano un límite de pérdida (p. ej., −10 %) y revisar el VaR/CVaR de la cartera. Tener presente que las colas reales son más pesadas que las gaussianas (Jarque-Bera rechaza normalidad).
  5. Rebalanceo y horizonte: al no existir predictibilidad, el retorno proviene de la prima de riesgo de largo plazo. Se sugiere un horizonte superior a 30 días y rebalanceo periódico.
  6. Costos: minimizar la rotación (comisiones, impuestos y spread), ya que sin ventaja predictiva cada operación adicional es un costo neto.

16 Limitaciones

  • Homocedasticidad: los ARIMA suponen varianza constante, pero los rendimientos de GOOGL presentan agrupamiento de volatilidad. Un modelo ARIMA-GARCH capturaría mejor el riesgo, aunque sin mejorar necesariamente el pronóstico de la media.
  • Normalidad: Jarque-Bera rechaza la normalidad de los residuos. Los intervalos paramétricos pueden ser poco fiables en las colas. El bootstrap lo mitiga, pero no lo elimina.
  • Evaluación de una sola partición: la validación usa un único origen (85/15). Una validación con ventana móvil (rolling origin) daría estimaciones más robustas y pruebas formales como Diebold-Mariano.
  • Pronóstico multipaso: las métricas evalúan todo el horizonte de validación (cientos de pasos) y no pronósticos a un paso, lo que penaliza más a cualquier modelo con dinámica.
  • Muestra de cinco años: abarca un régimen particular de mercado (tasas de interés, auge de la IA). Los resultados pueden no generalizarse.
  • Sin variables exógenas: no se incluyen regresores (macroeconomía, resultados trimestrales, índices). ARIMAX y SARIMAX se presentan teóricamente, pero no se estiman.
  • Supuestos de la estrategia: la sección de inversión ignora costos de transacción, impuestos, dividendos y feriados, y usa una tasa libre de riesgo supuesta.
  • Datos: el precio ajustado de Yahoo Finance puede ser revisado, por lo que las cifras exactas pueden variar entre ejecuciones.

17 Conclusiones

  1. La serie de precios de GOOGL en niveles es no estacionaria (ADF y KPSS), mientras que la serie diferenciada es estacionaria (ADF y KPSS). Se concluye que d = 1.
  2. De los seis modelos evaluados manualmente, el mejor ARIMA es ARIMA(1,1,0) (RMSE = 0.08669540; AIC = −5274.727; BIC = −5264.787), seguido de ARIMA(0,1,1), ARIMA(2,1,1), ARIMA(1,1,2), ARIMA(1,1,1) y ARIMA(2,1,2).
  3. El benchmark Random Walk presenta un RMSE ligeramente menor (0.08669098), por lo que la complejidad adicional de los ARIMA no aporta mejoras de pronóstico.
  4. El diagnóstico residual muestra que Ljung-Box no rechaza H0 (residuos compatibles con ruido blanco), aunque Jarque-Bera rechaza la normalidad, lo que confirma las colas pesadas propias de los activos financieros.
  5. Los resultados son consistentes con la hipótesis de eficiencia débil del mercado: los precios de GOOGL muestran un comportamiento compatible con una caminata aleatoria.
  6. La simulación Monte Carlo con bootstrap de residuos cuantifica la incertidumbre del pronóstico a 30 días. Dada la ausencia de predictibilidad, la estrategia de inversión debe priorizar la gestión del riesgo y la diversificación sobre la búsqueda de señales direccionales.
  7. Como extensión futura se recomienda incorporar modelos ARIMA-GARCH, validación con ventana móvil y regresores exógenos (ARIMAX/SARIMAX).

18 Referencias

Alphabet Inc. (2024). Annual report on Form 10-K for the fiscal year ended December 31, 2023. U.S. Securities and Exchange Commission. https://www.sec.gov/

Bollerslev, T. (1986). Generalized autoregressive conditional heteroskedasticity. Journal of Econometrics, 31(3), 307-327. https://doi.org/10.1016/0304-4076(86)90063-1

Box, G. E. P., & Cox, D. R. (1964). An analysis of transformations. Journal of the Royal Statistical Society: Series B (Methodological), 26(2), 211-252. https://doi.org/10.1111/j.2517-6161.1964.tb00553.x

Box, G. E. P., Jenkins, G. M., Reinsel, G. C., & Ljung, G. M. (2015). Time series analysis: Forecasting and control (5th ed.). John Wiley & Sons.

Dickey, D. A., & Fuller, W. A. (1979). Distribution of the estimators for autoregressive time series with a unit root. Journal of the American Statistical Association, 74(366), 427-431. https://doi.org/10.1080/01621459.1979.10482531

Efron, B. (1979). Bootstrap methods: Another look at the jackknife. The Annals of Statistics, 7(1), 1-26. https://doi.org/10.1214/aos/1176344552

Engle, R. F. (1982). Autoregressive conditional heteroscedasticity with estimates of the variance of United Kingdom inflation. Econometrica, 50(4), 987-1007. https://doi.org/10.2307/1912773

Fama, E. F. (1970). Efficient capital markets: A review of theory and empirical work. The Journal of Finance, 25(2), 383-417. https://doi.org/10.2307/2325486

Hyndman, R. J., & Athanasopoulos, G. (2021). Forecasting: Principles and practice (3rd ed.). OTexts. https://otexts.com/fpp3/

Hyndman, R. J., & Khandakar, Y. (2008). Automatic time series forecasting: The forecast package for R. Journal of Statistical Software, 27(3), 1-22. https://doi.org/10.18637/jss.v027.i03

Hyndman, R. J., & Koehler, A. B. (2006). Another look at measures of forecast accuracy. International Journal of Forecasting, 22(4), 679-688. https://doi.org/10.1016/j.ijforecast.2006.03.001

Jarque, C. M., & Bera, A. K. (1987). A test for normality of observations and regression residuals. International Statistical Review, 55(2), 163-172. https://doi.org/10.2307/1403192

Kwiatkowski, D., Phillips, P. C. B., Schmidt, P., & Shin, Y. (1992). Testing the null hypothesis of stationarity against the alternative of a unit root. Journal of Econometrics, 54(1-3), 159-178. https://doi.org/10.1016/0304-4076(92)90104-Y

Ljung, G. M., & Box, G. E. P. (1978). On a measure of lack of fit in time series models. Biometrika, 65(2), 297-303. https://doi.org/10.1093/biomet/65.2.297

Ryan, J. A., & Ulrich, J. M. (2024). quantmod: Quantitative financial modelling framework (R package). https://CRAN.R-project.org/package=quantmod

R Core Team. (2024). R: A language and environment for statistical computing. R Foundation for Statistical Computing. https://www.R-project.org/

Trapletti, A., & Hornik, K. (2024). tseries: Time series analysis and computational finance (R package). https://CRAN.R-project.org/package=tseries