Este proyecto busca predecir la evolución del ELO del joven ajedrecista Faustino Oro utilizando un modelo de Random Walk con Filtro de Kalman, estimación por máxima verosimilitud (EM/MLE) y suavizado Rauch–Tung–Striebel (RTS). Finalmente, se extiende el modelo para proyectar 10 períodos futuros con intervalos de confianza.
El objetivo es ilustrar cómo los modelos de estado oculto pueden capturar la dinámica estocástica de un sistema en evolución, incluso con una serie corta de datos.
# --- Importar CSV ---
datos <- read.csv("./RatingFaustino_2024.csv")
# Seleccionando la columna de ELO
y <- datos$Rating
# Visualización inicial del ELO
g1 <- ggplot(datos, aes(x = 1:length(y), y = y)) +
geom_line(color = "steelblue", size = 1) +
geom_point(color = "darkblue", size = 2, alpha = 0.7) +
labs(title = "Evolución del ELO de Faustino Oro",
x = "Tiempo",
y = "ELO observado") +
theme_minimal()
g1
El modelo de espacio de estados propuesto es un Random Walk con ruido gaussiano:
\[ \begin{aligned} x_t &= x_{t-1} + w_t, \quad w_t \sim N(0, Q) \\ y_t &= x_t + v_t, \quad v_t \sim N(0, R) \end{aligned} \]
Donde: - \(x_t\): ELO verdadero (estado latente) - \(y_t\): ELO observado - \(Q\): varianza del proceso - \(R\): varianza de observación
# --- Filtro de Kalman ---
kalman_filter <- function(y, Q, R) {
T <- length(y)
x_pred <- rep(NA, T)
P_pred <- rep(NA, T)
x_filt <- rep(NA, T)
P_filt <- rep(NA, T)
x_filt[1] <- y[1]
P_filt[1] <- 1
for (t in 2:T) {
# Predicción
x_pred[t] <- x_filt[t-1]
P_pred[t] <- P_filt[t-1] + Q
# Actualización
K <- P_pred[t] / (P_pred[t] + R)
x_filt[t] <- x_pred[t] + K * (y[t] - x_pred[t])
P_filt[t] <- (1 - K) * P_pred[t]
}
list(x_filt = x_filt, P_filt = P_filt, x_pred = x_pred, P_pred = P_pred)
}
Los parámetros \(Q\) y \(R\) se estiman maximizando la verosimilitud basada en innovaciones. La optimización se realiza en escala logarítmica para garantizar positividad.
# --- Log-Verosimilitud ---
log_likelihood <- function(params, y) {
Q <- exp(params[1]) # parámetros en escala log para mantener positividad
R <- exp(params[2])
T <- length(y)
x_prev <- y[1]
P_prev <- 1
ll <- 0
for (t in 2:T) {
# Predicción
x_pred <- x_prev
P_pred <- P_prev + Q
# Innovación
v <- y[t] - x_pred
S <- P_pred + R
ll <- ll - 0.5 * (log(2*pi*S) + (v^2)/S)
# Actualización
K <- P_pred / S
x_prev <- x_pred + K * v
P_prev <- (1 - K) * P_pred
}
return(-ll) # Negativo para minimización
}
# --- Estimación por MLE ---
fit <- optim(par = log(c(1, 1)), fn = log_likelihood, y = y, method = "L-BFGS-B")
Q_hat <- exp(fit$par[1])
R_hat <- exp(fit$par[2])
cat("\nEstimaciones MLE:\nQ =", Q_hat, "\nR =", R_hat, "\n")
##
## Estimaciones MLE:
## Q = 1550.211
## R = 0.001492023
# --- Filtro con parámetros estimados ---
res <- kalman_filter(y, Q_hat, R_hat)
Este paso refina las estimaciones del estado en cada período considerando información futura, proporcionando trayectorias más estables.
# --- Smoothing de Rauch–Tung–Striebel ---
kalman_smoother <- function(y, Q, R, kf) {
T <- length(y)
x_smooth <- rep(NA, T)
P_smooth <- rep(NA, T)
x_smooth[T] <- kf$x_filt[T]
P_smooth[T] <- kf$P_filt[T]
for (t in (T-1):1) {
A <- kf$P_filt[t] / (kf$P_pred[t+1])
x_smooth[t] <- kf$x_filt[t] + A * (x_smooth[t+1] - kf$x_pred[t+1])
P_smooth[t] <- kf$P_filt[t] + A^2 * (P_smooth[t+1] - kf$P_pred[t+1])
}
list(x_smooth = x_smooth, P_smooth = P_smooth)
}
sm <- kalman_smoother(y, Q_hat, R_hat, res)
Con el último estado suavizado, se proyectan 10 pasos hacia adelante y se construye una banda de confianza del 95%.
# --- Predicción futura (10 pasos) ---
predict_future <- function(last_state, last_P, Q, steps = 10) {
x_future <- rep(NA, steps)
P_future <- rep(NA, steps)
x_future[1] <- last_state
P_future[1] <- last_P + Q
for (t in 2:steps) {
x_future[t] <- x_future[t-1]
P_future[t] <- P_future[t-1] + Q
}
data.frame(
tiempo = (length(y)+1):(length(y)+steps),
predicho = x_future,
lower = x_future - 1.96 * sqrt(P_future),
upper = x_future + 1.96 * sqrt(P_future)
)
}
future_pred <- predict_future(sm$x_smooth[length(y)], sm$P_smooth[length(y)], Q_hat, 10)
El gráfico combina los datos observados, estimaciones filtradas, suavizadas y predicciones futuras.
# --- Resultados ---
resultados <- data.frame(
tiempo = 1:length(y),
observado = y,
filtrado = res$x_filt,
suavizado = sm$x_smooth
)
# --- Visualización ---
ggplot() +
geom_line(data = resultados, aes(x = tiempo, y = observado), color = "gray", size = 0.8, alpha = 0.6) +
geom_line(data = resultados, aes(x = tiempo, y = filtrado), color = "blue", size = 1) +
geom_line(data = resultados, aes(x = tiempo, y = suavizado), color = "red", size = 1) +
geom_line(data = future_pred, aes(x = tiempo, y = predicho), color = "darkgreen", size = 1, linetype = "dashed") +
geom_ribbon(data = future_pred, aes(x = tiempo, ymin = lower, ymax = upper), alpha = 0.2, fill = "darkgreen") +
labs(title = "Predicción del ELO de Faustino Oro (Filtro de Kalman + Suavizado + Predicción)",
y = "ELO", x = "Tiempo",
caption = "Azul: filtrado | Rojo: suavizado | Verde: predicción 10 pasos") +
theme_minimal()
A continuación se agregan visuales complementarios para enriquecer el análisis:
# 1. Dispersión del error de innovación
errores <- y - res$x_pred
g_error <- ggplot(data.frame(t = 1:length(errores), e = errores), aes(x = t, y = e)) +
geom_line(color = "tomato", size = 0.8) +
geom_hline(yintercept = 0, color = "gray50", linetype = "dashed") +
labs(title = "Errores de innovación (observado - predicho)", y = "Error", x = "Tiempo") +
theme_minimal()
# 2. Distribución de los errores
g_hist <- ggplot(data.frame(e = errores), aes(x = e)) +
geom_histogram(fill = "darkorange", color = "white", bins = 15, alpha = 0.7) +
labs(title = "Distribución de los errores de innovación", x = "Error", y = "Frecuencia") +
theme_minimal()
# 3. Evolución de la incertidumbre
g_var <- ggplot(data.frame(t = 1:length(res$P_filt), var = res$P_filt), aes(x = t, y = var)) +
geom_line(color = "purple", size = 0.8) +
labs(title = "Evolución de la incertidumbre (P_filt)", y = "Varianza del estado", x = "Tiempo") +
theme_minimal()
gridExtra::grid.arrange(g_error, g_hist, g_var, ncol = 1)
Este enfoque es extensible: se pueden agregar covariables, shocks estructurales o modelos de tendencia más complejos (p. ej., Local Linear Trend o ARIMA-Kalman).