K_GRID <- 1:6
LAMBDAS <- c(0.3, 0.1, 0.03, 0.01, 0.003)
ALPHAS <- c(0, 0.5, 1)
SVR_GRID <- tidyr::expand_grid(cost = c(0.1, 1, 10), epsilon = c(0.05, 0.2))
# Cada función devuelve una matriz 2 x p: filas = h3, h12; columnas = instrumentos
como_2xp <- function(M) { M <- matrix(M, nrow = 2); colnames(M) <- codigos; M }
f_naive <- function(t) {
u <- apply(Yobs[1:t, , drop = FALSE], 2, function(x) { v <- x[!is.na(x)]; v[length(v)] })
como_2xp(rbind(u, u))
}
f_arima <- function(t) {
como_2xp(sapply(seq_len(p), function(j) {
fit <- forecast::auto.arima(ts(Yf[1:t, j], frequency = 12), seasonal = FALSE,
allowdrift = FALSE, stepwise = TRUE, approximation = FALSE)
as.numeric(forecast::forecast(fit, h = 12)$mean)[c(3, 12)]
}))
}
f_var_dif <- function(t) {
dY <- diff(Yf[1:t, ])
sel <- vars::VARselect(dY, lag.max = 4, type = "none")$selection[["SC(n)"]]
fit <- vars::VAR(dY, p = max(1, sel), type = "none")
dfc <- sapply(predict(fit, n.ahead = 12)$fcst, function(m) m[, "fcst"])
niv <- sweep(apply(dfc, 2, cumsum), 2, Yf[t, ], "+")
como_2xp(niv[c(3, 12), ])
}
f_var_niv <- function(t) {
sel <- vars::VARselect(Yf[1:t, ], lag.max = 2, type = "const")$selection[["SC(n)"]]
fit <- vars::VAR(Yf[1:t, ], p = max(1, sel), type = "const")
como_2xp(sapply(predict(fit, n.ahead = 12)$fcst, function(m) m[c(3, 12), "fcst"]))
}
f_vecm <- function(t) {
cj <- urca::ca.jo(Yf[1:t, ], type = "trace", ecdet = "const", K = 2, spec = "transitory")
vv <- vars::vec2var(cj, r = r_vecm)
como_2xp(sapply(predict(vv, n.ahead = 12)$fcst, function(m) m[c(3, 12), "fcst"]))
}
# Factores: pronóstico + reconstrucción + residuo actual persistente
reconstruir_con_residuo <- function(Ff, aj, t) {
res_t <- Yf[t, ] - reconstruir(aj$F[t, , drop = FALSE], aj)
como_2xp(sweep(reconstruir(Ff, aj), 2, as.numeric(res_t), "+"))
}
f_pca_ar <- function(t, k) {
aj <- pca_ajuste(Yf[1:t, ], k)
Ff <- sapply(seq_len(k), function(i) {
f <- aj$F[, i]
tryCatch(as.numeric(predict(arima(f, order = c(1, 0, 0)), n.ahead = 12)$pred)[c(3, 12)],
error = function(e) rep(f[length(f)], 2))
})
reconstruir_con_residuo(matrix(Ff, nrow = 2), aj, t)
}
f_pca_var <- function(t, k) {
aj <- pca_ajuste(Yf[1:t, ], k)
if (k == 1) { # VAR univariado = AR(1)
return(f_pca_ar(t, 1))
}
fit <- vars::VAR(aj$F, p = 1, type = "const")
Ff <- sapply(predict(fit, n.ahead = 12)$fcst, function(m) m[c(3, 12), "fcst"])
reconstruir_con_residuo(matrix(Ff, nrow = 2), aj, t)
}
# ML ---------------------------------------------------------------------------
crear_features <- function(X) {
X <- as.matrix(X)
rez3 <- rbind(matrix(NA, 3, ncol(X)), X[1:(nrow(X) - 3), , drop = FALSE])
out <- cbind(X, X - rez3)
colnames(out) <- c(paste0(colnames(X), "_niv"), paste0(colnames(X), "_d3"))
out
}
ml_origen <- function(t) {
bloques <- list(raw = crear_features(Yf[1:t, , drop = FALSE]))
for (k in K_GRID) bloques[[paste0("k", k)]] <- crear_features(pca_ajuste(Yf[1:t, ], k)$F)
out <- list()
for (nb in names(bloques)) {
Fe <- bloques[[nb]]
xt <- Fe[t, , drop = FALSE]
etiqueta <- if (nb == "raw") "tasas" else "PCA"
for (h in H) {
if (t + h > N) next
s_idx <- 4:(t - h)
for (j in seq_len(p)) {
ok <- s_idx[!is.na(Yobs[s_idx + h, j])]
x <- Fe[ok, , drop = FALSE]
y <- Yobs[ok + h, j] - Yf[ok, j]
# Elastic Net
for (a in ALPHAS) {
pr <- tryCatch({
fit <- glmnet::glmnet(x, y, alpha = a, lambda = LAMBDAS, standardize = TRUE)
as.numeric(predict(fit, newx = xt, s = LAMBDAS))
}, error = function(e) rep(0, length(LAMBDAS)))
out[[length(out) + 1]] <- data.frame(
modelo = paste0("EN-", etiqueta), config = paste0(nb, "_a", a, "_l", LAMBDAS),
origen = t, h = h, codigo = codigos[j], pronostico = Yf[t, j] + pr)
}
# SVR
pr_svr <- vapply(seq_len(nrow(SVR_GRID)), function(g) {
tryCatch({
fit <- e1071::svm(x = x, y = y, type = "eps-regression", kernel = "radial",
cost = SVR_GRID$cost[g], epsilon = SVR_GRID$epsilon[g], scale = TRUE)
as.numeric(predict(fit, xt))
}, error = function(e) 0)
}, numeric(1))
out[[length(out) + 1]] <- data.frame(
modelo = paste0("SVR-", etiqueta),
config = paste0(nb, "_C", SVR_GRID$cost, "_e", SVR_GRID$epsilon),
origen = t, h = h, codigo = codigos[j], pronostico = Yf[t, j] + pr_svr)
}
}
}
bind_rows(out)
}
