Módulo 2·Regresión polinomial, logarítmica y exponencial

Modelos de Tendencia Determinística

Modelos polinomiales, log-polinomiales y exponencial-polinomiales para capturar la tendencia de largo plazo. Estimación por MCO, inferencia y diagnóstico de residuos.

Modelos polinomiales grado pModelo log-polinomialModelo exponencial-polinomialEstimación por MCO (lm)Diagnóstico de residuos

Anterior: Fundamentos y Descomposición Siguiente: Modelos de Estacionalidad Determinística

¿Por qué modelar la tendencia?

Cuando una serie de tiempo tiene una tendencia clara —sube año a año, o baja, o tiene una curva— los estadísticos y las pruebas de hipótesis clásicas dejan de ser válidas. La media no es constante, la varianza puede cambiar, y las autocorrelaciones están contaminadas por la tendencia. El primer paso siempre es remover la tendencia para poder estudiar lo que queda.

En este módulo aprenderás a modelar la tendencia —no solo a eliminarla— mediante polinomios y transformaciones logarítmicas. La diferencia es importante: un modelo de tendencia da una ecuación interpretable, permite extrapolar y da residuos con una distribución conocida. Un simple diferenciado quita la tendencia pero no te dice qué forma tenía.

2.1 Modelos según la descomposición

La forma funcional de la tendencia depende directamente del tipo de descomposición que elegiste en el Módulo 1. Si la serie es aditiva (varianza estable), la tendencia entra de forma lineal o polinomial en la ecuación. Si es multiplicativa (varianza crece con el nivel), trabajas con $\log(Y_t)$ y la tendencia también queda lineal o polinomial en esa escala.

Modelo polinomial (aditivo)

Para $Y_t = T_t + E_t$ , la tendencia se aproxima con un polinomio de grado $p$ :

$T_t = \beta_0 + \sum_{j=1}^{p} \beta_j \, t^j$

El modelo completo queda:

$Y_t = \beta_0 + \beta_1 t + \beta_2 t^2 + \cdots + \beta_p t^p + E_t, \quad E_t \overset{\text{iid}}{\sim} N(0, \sigma^2)$

Los parámetros $\beta_j$ se estiman por mínimos cuadrados ordinarios (MCO) con lm(). Su interpretación es directa: $\beta_1$ es la pendiente lineal de la tendencia por unidad de tiempo; $\beta_2$ controla la curvatura.

Modelo log-polinomial (multiplicativo)

Para la descomposición multiplicativa $Y_t = T_t \cdot \exp(\varepsilon_t)$ , donde el error actúa proporcionalmente al nivel, aplicamos logaritmo natural y obtenemos un modelo lineal en $W_t = \log(Y_t)$ :

$W_t = \log(Y_t) = \underbrace{\beta_0 + \sum_{j=1}^{p} \beta_j \, t^j}_{\log T_t} + \varepsilon_t, \quad \varepsilon_t \overset{\text{iid}}{\sim} N(0, \sigma^2)$

El modelo transformado es lineal y se ajusta con lm() sobre $W_t$ . Para volver a la escala original se aplica la corrección por sesgo de la distribución log-normal:

$\hat{Y}_t = \exp\!\left(\hat{W}_t + \frac{MSE}{2}\right) = \exp\!\left(\widehat{\log Y_t}\right) \cdot \exp\!\left(\frac{MSE}{2}\right)$

¿Por qué exp(MSE/2)? — La intuición

Cuando $W_t \sim N(\mu_t, \sigma^2)$ , la variable $Y_t = \exp(W_t)$ sigue una distribución log-normal cuya media es:

$E[Y_t] = \exp\!\left(\mu_t + \frac{\sigma^2}{2}\right)$

Sustituimos $\mu_t \leftarrow \hat{W}_t$ y $\sigma^2 \leftarrow MSE$ y obtenemos el estimador correcto de $E[Y_t]$ . Omitir el factor $\exp(MSE/2)$ es estimar la mediana de la distribución, no la media. Si $MSE = 0.04$ la diferencia es ~2% (tolerable); si $MSE = 0.25$ es ~13% (ya importa). La corrección siempre sube el pronóstico respecto al naive $\exp(\hat{W}_t)$ .

Modelo exponencial-polinomial (parcialmente multiplicativo)

Cuando se desea que la tendencia sea intrínsecamente exponencial pero con error aditivo en la escala original:

$Y_t = \exp\!\left(\beta_0 + \sum_{j=1}^{p} \beta_j \, t^j\right) + E_t, \quad E_t \overset{\text{iid}}{\sim} N(0, \sigma^2)$

La diferencia con el modelo log-polinomial es sutil pero importante: aquí la varianza de $E_t$ es constante en escala original (no en escala log). La estimación requiere mínimos cuadrados no lineales —función nls() en R o la función regexponencialv02() del repositorio de la profesora Nelfi González.

2.2 Estimación por mínimos cuadrados

Para el modelo polinomial, la ecuación en forma matricial es:

$\mathbf{Y} = \mathbf{X}\boldsymbol{\beta} + \mathbf{E}$

donde $\mathbf{X}$ es la matriz de Vandermonde:

$\mathbf{X} = \begin{pmatrix} 1 & t_1 & t_1^2 & \cdots & t_1^p \\ 1 & t_2 & t_2^2 & \cdots & t_2^p \\ \vdots & & & & \vdots \\ 1 & t_n & t_n^2 & \cdots & t_n^p \end{pmatrix}$

El estimador de mínimos cuadrados es:

$\hat{\boldsymbol{\beta}} = (\mathbf{X}^T \mathbf{X})^{-1} \mathbf{X}^T \mathbf{Y}$

Cómo leer el summary(lm()) en R

Cuando ejecutas summary(mod), los números clave son:

Estimate (β̂): el valor ajustado de cada coeficiente.
Std. Error: la incertidumbre del estimador.
Pr(>|t|): el p-valor de la prueba $H_0: \beta_j = 0$ . Si es menor a 0.05, el término aporta información real.
Residual standard error: es $\hat{\sigma}$ , la raíz del MSE. $MSE = \hat{\sigma}^2$ , que necesitas para la corrección log-normal.
Adjusted R²: proporción de varianza explicada, penalizada por el número de parámetros. Más útil que el R² simple para comparar modelos de distinto grado.

Parsimonia: el principio que guía la selección de modelos

En series de tiempo —como en estadística en general— existe una tensión entre ajuste y complejidad. Un polinomio de grado 10 siempre se ajusta mejor al conjunto de entrenamiento que uno de grado 2, pero sus predicciones fuera de la muestra suelen ser peores: está aprendiendo el ruido, no la señal.

El principio de parsimonia (navaja de Occam aplicada a modelos) dice: el mejor modelo es el más simple que explica adecuadamente los datos. En la práctica, esto se operacionaliza con AIC y BIC.

2.3 AIC, BIC y selección de modelos

AIC (Criterio de Información de Akaike) y BIC (Criterio Bayesiano de Schwarz) son las herramientas estándar para comparar modelos con distinto número de parámetros. Penalizan la log-verosimilitud por la complejidad del modelo:

$\text{AIC} = -2\ell(\hat{\boldsymbol{\beta}}) + 2k \qquad \text{BIC} = -2\ell(\hat{\boldsymbol{\beta}}) + k\log(n)$

donde $\ell$ es la log-verosimilitud maximizada, $k$ es el número de parámetros estimados y $n$ es el tamaño de la muestra. Un AIC/BIC más bajo indica un mejor modelo. BIC penaliza más fuerte la complejidad (especialmente con $n$ grande), por lo que tiende a elegir modelos más parsimoniosos.

¿Cuándo usar AIC vs BIC?

Usa AIC cuando el objetivo principal es la predicción —quieres capturar toda la señal posible aunque el modelo sea un poco más complejo.

Usa BIC cuando el objetivo es identificar el verdadero proceso generador de datos —quieres el modelo más simple que sea consistente con los datos.

En la práctica, si AIC y BIC apuntan al mismo modelo, la decisión está clara. Si discrepan, analiza los residuos de ambos modelos y decide con criterio estadístico.

2.4 Ejemplo A — Lago Hurón: tendencia polinomial

El dataset LakeHuron (base R) registra el nivel anual del lago Hurón en pies sobre el nivel del mar, de 1875 a 1972 ( $n = 98$ observaciones). La serie muestra una tendencia decreciente con cierta curvatura — un caso ideal para comparar polinomios de distintos grados con AIC/BIC.

¿Qué vas a ver al ejecutar este código?

Una tabla comparando grados 1 a 4: AIC, BIC, R² ajustado y RMSE. Busca el grado donde AIC/BIC deja de bajar —ahí está el punto de equilibrio entre ajuste y parsimonia.
El resumen del modelo de grado 2: coeficientes, p-valores y R² ajustado.
Un gráfico de la serie original con la tendencia cuadrática superpuesta.

▶ LakeHuron — ajuste polinomial con Mipoly() y selección por AIC/BIC(ggplot2)

# ── Función auxiliar Mipoly (autocontenida) ───────────────
# Construye la matriz [t, t^2, ..., t^p] para usar con lm()
Mipoly <- function(tiempo, grado) {
df <- data.frame(matrix(nrow = length(tiempo), ncol = grado))
for (j in 1:grado) df[, j] <- tiempo^j
names(df) <- paste0("t", 1:grado)
df
}

# ── Dataset ───────────────────────────────────────────────
yt <- LakeHuron              # base R — serie ts anual, frequency = 1
n  <- length(yt)             # 98 observaciones
t  <- 1:n

# ── Comparar grados 1 a 4 con AIC, BIC, R²adj, RMSE ──────
criterios <- lapply(1:4, function(g) {
mod <- lm(as.numeric(yt) ~ ., data = Mipoly(t, g))
data.frame(
  Grado = g,
  AIC   = round(AIC(mod), 2),
  BIC   = round(BIC(mod), 2),
  R2adj = round(summary(mod)$adj.r.squared, 4),
  RMSE  = round(sqrt(mean(residuals(mod)^2)), 3)
)
})
cat("── Comparación de modelos (menor AIC/BIC = mejor) ──\n")
print(do.call(rbind, criterios))

# ── Modelo seleccionado: grado 2 ──────────────────────────
cat("\n── Modelo grado 2 (seleccionado por parsimonia) ───\n")
mod_p2 <- lm(as.numeric(yt) ~ ., data = Mipoly(t, grado = 2))
summary(mod_p2)

# ── Gráfico con ggplot2 ───────────────────────────────────
library(ggplot2)

df_plot <- data.frame(
Año    = as.numeric(time(yt)),
Obs    = as.numeric(yt),
Ajuste = fitted(mod_p2)
)

ggplot(df_plot, aes(x = Año)) +
geom_line(aes(y = Obs,    colour = "Observado"),
          linewidth = 0.9, alpha = 0.85) +
geom_line(aes(y = Ajuste, colour = "Tendencia cuadrática"),
          linewidth = 1.4) +
scale_colour_manual(
  values = c("Observado" = "#78716c", "Tendencia cuadrática" = "#1d4ed8")
) +
labs(
  title    = "Lago Hurón — Tendencia cuadrática (1875–1972)",
  subtitle = expression(hat(Y)[t] == hat(beta)[0] + hat(beta)[1]*t + hat(beta)[2]*t^2),
  x = "Año", y = "Nivel (pies)", colour = NULL
) +
theme_bw(base_size = 12) +
theme(legend.position = "top")

Ctrl+Enter

En la tabla de criterios: si el AIC baja de grado 1 a 2 pero apenas cambia de 2 a 3, el principio de parsimonia indica quedarse con grado 2. Agregar más términos cuando AIC no mejora significativamente solo añade complejidad sin beneficio real.

Cómo interpretar la tabla de criterios

Para LakeHuron es típico observar que el AIC cae claramente de grado 1 a grado 2 (la curvatura es real), pero apenas cambia de grado 2 a 3. El BIC, que penaliza más fuerte la complejidad, suele señalar directamente el grado 2. El R² ajustado también debería aumentar al pasar a grado 2 y estabilizarse después.

Recuerda: un R² alto no es suficiente. Un polinomio de grado 8 tendrá R² = 0.99 pero sus residuos mostrarán oscilaciones salvajes (sobreajuste). El diagnóstico de residuos es el criterio definitivo.

▶ Visualización del AIC/BIC por grado — curva de selección(ggplot2)

library(ggplot2)

grados    <- 1:4
lista_mod <- lapply(grados, function(g)
lm(as.numeric(yt) ~ ., data = Mipoly(t, g))
)

criterios <- data.frame(
Grado = grados,
AIC   = round(sapply(lista_mod, AIC), 2),
BIC   = round(sapply(lista_mod, BIC), 2),
R2adj = round(sapply(lista_mod, function(m) summary(m)$adj.r.squared), 4),
RMSE  = round(sapply(lista_mod, function(m) sqrt(mean(residuals(m)^2))), 3)
)
print(criterios)

# AIC por grado: curva de selección
ggplot(criterios, aes(x = Grado, y = AIC)) +
geom_line(colour = "#1d4ed8", linewidth = 1.2) +
geom_point(colour = "#1d4ed8", size = 3) +
geom_point(data = criterios[which.min(criterios$AIC), ],
           colour = "#dc2626", size = 5, shape = 8) +
labs(title    = "AIC según grado del polinomio — LakeHuron",
     subtitle = "Punto rojo = grado óptimo (mínimo AIC)",
     x = "Grado p", y = "AIC") +
theme_bw(base_size = 12)

Ctrl+Enter

Este tipo de gráfico es el que usarás en la práctica cada vez que compares modelos: identifica visualmente dónde la curva 'se aplana' — ese es el punto de parsimonia.

2.5 Ejemplo B — CO₂ Mauna Loa: modelo log-lineal

El dataset co2 (base R) contiene la concentración mensual de CO₂ atmosférico medida en el observatorio de Mauna Loa, Hawaii, entre enero de 1959 y diciembre de 1997 ( $n = 468$ observaciones). La serie tiene crecimiento prácticamente exponencial — el modelo log-lineal captura la tendencia de largo plazo. La estacionalidad se añade en el Módulo 3.

El modelo ajustado sobre $\log(Y_t)$ es:

$\widehat{\log(Y_t)} = \hat{\beta}_0 + \hat{\beta}_1\, t$

Y en la escala original, con corrección por sesgo:

$\hat{Y}_t = \exp\!\left(\hat{\beta}_0 + \hat{\beta}_1\, t\right) \cdot \exp\!\left(\frac{MSE}{2}\right)$

El coeficiente $\hat{\beta}_1$ en escala log se interpreta como la tasa de crecimiento mensual continua del CO₂. La tasa porcentual mensual es $(\exp(\hat{\beta}_1) - 1) \times 100\%$ .

▶ co2 — modelo log-lineal con corrección por sesgo(ggplot2)

# ── Función auxiliar ─────────────────────────────────────
Mipoly <- function(tiempo, grado) {
df <- data.frame(matrix(nrow = length(tiempo), ncol = grado))
for (j in 1:grado) df[, j] <- tiempo^j
names(df) <- paste0("t", 1:grado); df
}

# ── Dataset ───────────────────────────────────────────────
yt2 <- co2                   # base R — mensual 1959–1997
n2  <- length(yt2)           # 468 observaciones
t2  <- 1:n2

# ── Modelo log-lineal (tendencia exponencial) ─────────────
mod_log <- lm(log(as.numeric(yt2)) ~ ., data = Mipoly(t2, grado = 1))
mse_log <- summary(mod_log)$sigma^2        # MSE en escala ln

# Corrección por sesgo de la distribución log-normal
ythat_log <- exp(fitted(mod_log)) * exp(mse_log / 2)

# Interpretación del coeficiente
tasa_mensual <- (exp(coef(mod_log)["t1"]) - 1) * 100
cat(sprintf("Tasa de crecimiento mensual: %.4f%%\n", tasa_mensual))
cat(sprintf("Tasa de crecimiento anual:   %.2f%%\n", tasa_mensual * 12))
cat(sprintf("MSE (escala log):            %.6f\n", mse_log))
cat(sprintf("Factor de corrección Duan:   %.6f\n", exp(mse_log / 2)))

# ── Gráfico ggplot2 ───────────────────────────────────────
library(ggplot2)

df_co2 <- data.frame(
Año    = as.numeric(time(yt2)),
Obs    = as.numeric(yt2),
Ajuste = ythat_log
)

ggplot(df_co2, aes(x = Año)) +
geom_line(aes(y = Obs,    colour = "Observado"),
          linewidth = 0.7, alpha = 0.8) +
geom_line(aes(y = Ajuste, colour = "Tendencia log-lineal"),
          linewidth = 1.4) +
scale_colour_manual(
  values = c("Observado"          = "#78716c",
             "Tendencia log-lineal" = "#dc2626")
) +
labs(
  title    = "CO₂ atmosférico — Mauna Loa (1959–1997)",
  subtitle = "log(Yₜ) = β₀ + β₁t + Eₜ  |  corrección exp(MSE/2) aplicada",
  x = "Año", y = "CO₂ (ppm)", colour = NULL
) +
theme_bw(base_size = 12) +
theme(legend.position = "top")

Ctrl+Enter

Los residuos de este modelo mostrarán un ciclo anual muy claro: la curva roja captura solo la tendencia. Ese ciclo es la estacionalidad, que modelaremos con variables dummy y series de Fourier en el Módulo 3.

¿Qué deberías ver en los resultados?

La tasa de crecimiento mensual del CO₂ estará alrededor de 0.12% mensual, lo que equivale a ~1.4% anual. El factor de corrección de Duan $\exp(MSE/2)$ estará muy cerca de 1.0 —el MSE en escala log es pequeño porque la tendencia es muy limpia. La gráfica mostrará la curva roja siguiendo el centro de la serie, ignorando las oscilaciones anuales —eso es intencional.

2.6 Inferencia y selección del grado

Para determinar el grado $p$ del polinomio, añadimos términos de mayor grado y verificamos si el último coeficiente es estadísticamente significativo.

Hipótesis para el término de grado $p$ :

$H_0: \beta_p = 0 \quad \text{vs} \quad H_1: \beta_p \neq 0$

El estadístico de prueba sigue una distribución $t$ con $n - (p+1)$ grados de libertad bajo $H_0$ . Si el p-valor es grande (por convención, mayor a 0.05), reducimos el grado en uno. Si es pequeño, el término aporta información real.

Error frecuente: elegir el grado por R² solamente

El R² siempre sube cuando añades un término al modelo —eso no significa que el modelo mejore. Un R² = 0.95 con grado 3 puede ser peor que un R² = 0.93 con grado 2 si el término cúbico tiene p-valor = 0.8. Usa R² ajustado y AIC/BIC para comparar, no el R² simple.

2.7 Diagnóstico de residuos

El diagnóstico de residuos es el paso más importante del análisis —y el más frecuentemente omitido. Un buen ajuste en entrenamiento no garantiza nada si los residuos tienen estructura. Verificamos cuatro condiciones sobre $\hat{E}_t = Y_t - \hat{Y}_t$ :

| Supuesto | Herramienta visual | Test formal | |---|---|---| | Media cero | $\hat{E}_t$ vs $t$ | — | | Varianza constante | $\hat{E}_t$ vs $\hat{Y}_t$ | Breusch-Pagan | | Independencia | ACF / PACF de residuos | Durbin-Watson, Ljung-Box | | Normalidad | Q-Q plot | Shapiro-Wilk |

El orden de importancia de los supuestos

En series de tiempo, el supuesto más crítico es la independencia (no autocorrelación de residuos). Si los residuos están autocorrelacionados, los errores estándar de los coeficientes son incorrectos, las pruebas t no son válidas y los intervalos de confianza no tienen el nivel nominal. La normalidad es el supuesto menos crítico —los estimadores MCO son robustos a desviaciones moderadas de normalidad gracias al teorema central del límite.

▶ Diagnóstico completo de residuos — LakeHuron (grado 2)(ggplot2)

library(ggplot2)

Mipoly <- function(tiempo, grado) {
df <- data.frame(matrix(nrow = length(tiempo), ncol = grado))
for (j in 1:grado) df[, j] <- tiempo^j
names(df) <- paste0("t", 1:grado); df
}
yt    <- LakeHuron
n     <- length(yt)
t     <- 1:n
mod_p2 <- lm(as.numeric(yt) ~ ., data = Mipoly(t, grado = 2))
res    <- residuals(mod_p2)

# ── Tests formales ────────────────────────────────────────
sw <- shapiro.test(res)
lb <- Box.test(res, lag = 10, type = "Ljung-Box")
cat(sprintf("Shapiro-Wilk:  W=%.4f  p=%.4f  → %s\n",
  sw$statistic, sw$p.value,
  ifelse(sw$p.value > 0.05, "normalidad OK", "NO normal")))
cat(sprintf("Ljung-Box(10): Q=%.4f  p=%.4f  → %s\n",
  lb$statistic, lb$p.value,
  ifelse(lb$p.value > 0.05,
         "sin autocorrelacion signif.",
         "autocorrelacion detectada → revisar modelo")))

# ── Residuos vs Tiempo ────────────────────────────────────
df_res <- data.frame(t = t, res = res)
ggplot(df_res, aes(x = t, y = res)) +
geom_line(colour = "#3b82f6", linewidth = 0.8) +
geom_hline(yintercept = 0, colour = "#ef4444", linetype = "dashed") +
geom_hline(yintercept = c(-2, 2)*sd(res), colour = "#94a3b8",
           linetype = "dotted") +
labs(title    = "Residuos vs Tiempo — LakeHuron grado 2",
     subtitle = sprintf("SW p=%.3f  |  LB(10) p=%.3f  |  SD residuos=%.3f",
                        sw$p.value, lb$p.value, sd(res)),
     x = "t (período relativo)", y = expression(hat(E)[t])) +
theme_bw(base_size = 12)

Ctrl+Enter

Las líneas punteadas indican ±2σ — aproximadamente el 95% de los residuos debería caer dentro. Si ves un patrón sistemático (curva ascendente, embudo, o ciclos), el modelo no capturó toda la estructura de la serie.

▶ ACF de residuos — ¿hay autocorrelación?(ggplot2)

Ctrl+Enter

Si hay barras rojas en los primeros rezagos, los residuos no son ruido blanco. Para LakeHuron es común encontrar autocorrelación positiva en los primeros rezagos — señal de que hay un ciclo lento (decadal) que el modelo de tendencia no capturó. Esto se resolvería con un componente ARIMA para el error.

Patrón cíclico en los residuos: diagnóstico y solución

Es muy frecuente encontrar ciclos en los residuos de un modelo de tendencia. Esto indica: (1) correlación positiva entre errores consecutivos (viola independencia), o (2) que hay estacionalidad no modelada.

La solución depende del diagnóstico:

Si el ciclo tiene el mismo período que la estacionalidad (por ejemplo, 12 meses para datos mensuales) → añadir componente estacional (Módulo 3).
Si el ciclo es de largo plazo (varios años) → se necesita un componente ARIMA para el error (módulos avanzados).

Figura 2.1 — Residuos con ciclos (izq.) y ACF con barras significativas en rojo (der.) indican correlación serial: el modelo de tendencia solo no es suficiente.

Conexión con el Módulo 3

Los residuos del modelo de tendencia para co2 muestran un ciclo anual perfecto: el modelo no capturó la estacionalidad. En el Módulo 3 añadirás variables dummy estacionales (una por mes) o series de Fourier para modelar ese patrón. El resultado será un modelo con tendencia + estacionalidad cuyos residuos finalmente deberían comportarse como ruido blanco.