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:
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.
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:
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).
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
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í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 |
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.
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)| 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.
El flujo de trabajo sigue la metodología Box-Jenkins (Box et al., 2015):
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.
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}\)).
\[ 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.
\[ 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.
\[ \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 . \]
Incorpora regresores exógenos \(x_{t}\):
\[ \phi(B)\,(1-B)^d\,y_t = c + \boldsymbol{\beta}^{\prime}\mathbf{x}_t + \theta(B)\,\varepsilon_t . \]
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\).
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).
\[ 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.
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().
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)| 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:
###############################################
# 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.
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")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.
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)| 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 |
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)| 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.
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")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)| 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 |
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)
AIC = -3873.91 | BIC = -3863.97
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.
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)| 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)Conclusiones del diagnóstico.
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)| 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}\).
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)| 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.
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)| % 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.
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