Módulo 11·Modelos estocásticos para series de tiempo

Procesos ARMA, SARIMA y Regresión con Errores ARMA

Procesos AR(p), MA(q), ARMA(p,q) y SARIMA(p,d,q)(P,D,Q). Identificación por Box-Jenkins (ACF, PACF, EACF), ajuste por máxima verosimilitud, diagnóstico de residuos y pronóstico. Regresión con errores ARMA (SARIMAX) como modelo unificador. Ejemplo completo: salario real colombiano.

Operador de retardo BProceso AR(p) y condición de estacionariedadProceso MA(q) y condición de invertibilidadProceso ARMA(p,q)Modelo SARIMA(p,d,q)(P,D,Q)[S]Regresión con errores ARMA (SARIMAX)Identificación Box-JenkinsPronóstico con SARIMA

Anterior: Procesos Estocásticos Estacionarios y Funciones de Autocorrelación

SARIMAX: el modelo que engloba todos los anteriores

El modelo SARIMAX $(p,d,q)(P,D,Q)_S \times \mathbf{X}$ es la generalización más amplia de la familia ARMA. Sus casos particulares son:

AR(p):

(p,0,0)(0,0,0)

MA(q):

(0,0,q)(0,0,0)

ARMA(p,q):

(p,0,q)(0,0,0)

ARIMA(p,d,q):

(p,d,q)(0,0,0)

SARMA:

(p,0,q)(P,0,Q)_S

SARIMA:

(p,d,q)(P,D,Q)_S

Al agregar regresores externos $\mathbf{X}$ (componente X de SARIMAX), se obtiene el modelo de regresión con errores SARIMA, que es el marco unificador del análisis de series de tiempo.

Definición formal de los procesos AR(p), MA(q), ARMA(p,q) y SARIMA(p,d,q)(P,D,Q). Identificación mediante ACF, PACF y EACF. Ajuste por máxima verosimilitud, diagnóstico de residuos y pronóstico. Regresión con errores ARMA: el modelo del salario real colombiano.

11.1 Notación con el operador de retardo B

El operador de retardo $B$ se define como:

$B \, Y_t = Y_{t-1}, \quad B^k \, Y_t = Y_{t-k}$

Con esta notación, las operaciones fundamentales se escriben:

Diferencia regular

(1 - B) \, Y_t = Y_t - Y_{t-1} = \nabla Y_t

Diferencia estacional

(1 - B^S) \, Y_t = Y_t - Y_{t-S} = \nabla_S Y_t

El polinomio autorregresivo de orden $p$ :

$\phi(B) = 1 - \phi_1 B - \phi_2 B^2 - \cdots - \phi_p B^p$

El polinomio de media móvil de orden $q$ :

$\theta(B) = 1 + \theta_1 B + \theta_2 B^2 + \cdots + \theta_q B^q$

Y los polinomios estacionales análogos:

$\Phi(B^S) = 1 - \Phi_1 B^S - \cdots - \Phi_P B^{PS}$

$\Theta(B^S) = 1 + \Theta_1 B^S + \cdots + \Theta_Q B^{QS}$

11.2 Proceso autorregresivo AR(p)

Un proceso AR(p) expresa $Y_t$ como una combinación lineal de sus $p$ valores pasados más un error:

$Y_t = \phi_1 Y_{t-1} + \phi_2 Y_{t-2} + \cdots + \phi_p Y_{t-p} + a_t$

En notación de operador:

$\phi(B) \, Y_t = a_t$

Condición de estacionariedad del AR(p)

El proceso AR(p) es estacionario en covarianza si y solo si todas las raíces del polinomio $\phi(z) = 0$ están fuera del círculo unitario (módulo > 1):

$|z_i| > 1 \quad \text{para toda raíz } z_i \text{ de } \phi(z) = 0$

Caso particular: AR(1)

$Y_t = \phi_1 Y_{t-1} + a_t$

La condición de estacionariedad es $|\phi_1| < 1$ . Las propiedades clave son:

Media: $E[Y_t] = 0$ (o $\mu$ si se incluye intercepto)
Varianza: $\gamma(0) = \sigma^2_a / (1 - \phi_1^2)$
ACF: $\rho(k) = \phi_1^k$ (decae exponencialmente)
PACF: se corta en $k = 1$ ( $\phi(1) = \phi_1$ , $\phi(k) = 0$ para $k > 1$ )

AR(1) con raíz unitaria

Si $\phi_1 = 1$ , el proceso es una caminata aleatoria (random walk): $Y_t = Y_{t-1} + a_t$ . No es estacionario (la varianza crece con $t$ ). Este es el caso $d = 1$ en ARIMA: se necesita diferenciar.

▶ Ajuste AR(1) — serie simulada Z1t(forecast, lmtest)

library(forecast); library(lmtest)
set.seed(123)
Z1t <- ts(arima.sim(list(ar = 0.7), n = 100), freq = 1)
par(mfrow = c(1, 2))
acf(Z1t,  lag.max = 25, ci.type = "ma", main = "ACF de Z1t", lwd = 2)
pacf(Z1t, lag.max = 25, main = "PACF de Z1t", lwd = 2)
par(mfrow = c(1, 1))
ajuste <- Arima(Z1t, order = c(1, 0, 0), include.mean = FALSE, method = "ML")
coeftest(ajuste, df = length(Z1t) - 1)
par(mfrow = c(2, 3))
plot(residuals(ajuste)); abline(h = c(-2*sqrt(ajuste$sigma2), 0, 2*sqrt(ajuste$sigma2)))
plot(as.numeric(fitted(ajuste)), residuals(ajuste)); abline(h = c(-2*sqrt(ajuste$sigma2), 0, 2*sqrt(ajuste$sigma2)))
acf(residuals(ajuste), ci.type = "ma", lag.max = 25)
pacf(residuals(ajuste), lag.max = 25)
qqnorm(residuals(ajuste)); qqline(residuals(ajuste))
par(mfrow = c(1, 1))
print(Box.test(residuals(ajuste), lag = 25, type = "Ljung-Box"))
print(shapiro.test(residuals(ajuste)))
fc <- forecast(ajuste, h = 5, level = 95)
print(fc)
plot(fc, main = "Pronóstico AR(1) — Z1t")
lines(fitted(ajuste), col = 2)

Ctrl+Enter

Arima() del paquete forecast ajusta por máxima verosimilitud (method='ML'). coeftest() del paquete lmtest da los p-valores.

11.3 Proceso de media móvil MA(q)

Un proceso MA(q) expresa $Y_t$ como una combinación lineal de los $q$ errores pasados:

$Y_t = a_t + \theta_1 a_{t-1} + \theta_2 a_{t-2} + \cdots + \theta_q a_{t-q}$

En notación de operador:

$Y_t = \theta(B) \, a_t$

Condición de invertibilidad del MA(q)

El proceso MA(q) es invertible si todas las raíces de $\theta(z) = 0$ están fuera del círculo unitario. La invertibilidad garantiza que el proceso puede expresarse como un AR de orden infinito, lo que permite estimar los $a_t$ recursivamente a partir de observaciones.

Caso particular: MA(1)

$Y_t = a_t + \theta_1 a_{t-1}$

Propiedades:

Media: $E[Y_t] = 0$
Varianza: $\gamma(0) = \sigma^2_a (1 + \theta_1^2)$
ACF: $\rho(1) = \theta_1 / (1 + \theta_1^2)$ , $\rho(k) = 0$ para $k > 1$ (se corta en 1)
PACF: decae exponencialmente (no se corta)

▶ Ajuste MA(1) — serie simulada Z4t(forecast, lmtest)

Ctrl+Enter

Para MA(1), el pronóstico converge rápidamente a la media (cero) porque los errores futuros no son observables.

11.4 Proceso ARMA(p,q)

El proceso ARMA(p,q) combina la estructura autorregresiva y de media móvil:

$Y_t = \phi_1 Y_{t-1} + \cdots + \phi_p Y_{t-p} + a_t + \theta_1 a_{t-1} + \cdots + \theta_q a_{t-q}$

En notación compacta:

$\phi(B) \, Y_t = \theta(B) \, a_t$

Ventajas del modelo ARMA

Parsimonia: un ARMA(1,1) con 2 parámetros puede modelar lo que un AR(p) o MA(q) necesitaría muchos más.
Flexibilidad: captura tanto la inercia (AR) como los choques persistentes (MA).
Estacionariedad: si $d = 0$ , el proceso es estacionario (raíces de $\phi(z)$ fuera del círculo unitario).

Identificación de órdenes (p,q)

El procedimiento de Box-Jenkins para identificar los órdenes $(p, q)$ es:

Graficar la ACF y PACF muestrales.
Identificar el patrón (corte vs. decaimiento).
Confirmar con la EACF (función de autocorrelación extendida).
Validar con auto.arima() por criterio AIC/BIC.
Ajustar el modelo candidato y verificar residuos.

▶ Identificación y ajuste ARMA(1,2) — serie simulada Z3t(forecast, TSA, lmtest)

library(forecast); library(TSA); library(lmtest)
set.seed(123)
Z3t <- ts(arima.sim(list(ar = 0.5, ma = c(0.4, 0.3)), n = 100), freq = 1)
par(mfrow = c(1, 2))
acf(Z3t,  lag.max = 25, ci.type = "ma", main = "ACF de Z3t", lwd = 2)
pacf(Z3t, lag.max = 25, main = "PACF de Z3t", lwd = 2)
par(mfrow = c(1, 1))
eacf(Z3t, ar.max = 13, ma.max = 13)
auto.arima(Z3t, ic = "bic")
ajuste <- Arima(Z3t, order = c(1, 0, 2), include.mean = FALSE, method = "ML")
coeftest(ajuste, df = length(Z3t) - 3)
par(mfrow = c(2, 3))
plot(residuals(ajuste)); abline(h = c(-2*sqrt(ajuste$sigma2), 0, 2*sqrt(ajuste$sigma2)))
plot(as.numeric(fitted(ajuste)), residuals(ajuste)); abline(h = c(-2*sqrt(ajuste$sigma2), 0, 2*sqrt(ajuste$sigma2)))
acf(residuals(ajuste), ci.type = "ma", lag.max = 25)
pacf(residuals(ajuste), lag.max = 25)
qqnorm(residuals(ajuste)); qqline(residuals(ajuste))
par(mfrow = c(1, 1))
Box.test(residuals(ajuste), lag = 25, type = "Ljung-Box")
shapiro.test(residuals(ajuste))
fc <- forecast(ajuste, h = 5, level = 95)
print(fc)
plot(fc, main = "Pronóstico ARMA(1,2) — Z3t")
lines(fitted(ajuste), col = 2)

Ctrl+Enter

El procedimiento Box-Jenkins: identificar → estimar → diagnosticar. Si los residuos no son ruido blanco, iterar.

11.5 Modelo SARIMA(p,d,q)(P,D,Q)[S]

El modelo SARIMA extiende ARIMA incorporando estructura estacional. La especificación completa es:

$\phi(B) \, \Phi(B^S) \, (1-B)^d \, (1-B^S)^D \, Y_t = \theta(B) \, \Theta(B^S) \, a_t$

donde:

$(p, d, q)$ : órdenes no estacionales (AR, diferenciación, MA)
$(P, D, Q)$ : órdenes estacionales
$S$ : período estacional (12 mensual, 4 trimestral, 52 semanal)
$\phi(B), \theta(B)$ : polinomios AR y MA no estacionales
$\Phi(B^S), \Theta(B^S)$ : polinomios AR y MA estacionales
$(1-B)^d$ : diferencia regular de orden $d$
$(1-B^S)^D$ : diferencia estacional de orden $D$

Casos particulares de SARIMA

| Modelo | Especificación | Condición | |---|---|---| | AR(p) | $(p,0,0)(0,0,0)$ | Sin MA, sin diferenciación, sin estacionalidad | | MA(q) | $(0,0,q)(0,0,0)$ | Sin AR, sin diferenciación, sin estacionalidad | | ARMA(p,q) | $(p,0,q)(0,0,0)$ | Estacionario, sin estacionalidad | | ARIMA(p,d,q) | $(p,d,q)(0,0,0)$ | Con diferenciación, sin estacionalidad | | SARMA | $(p,0,q)(P,0,Q)_S$ | Estacionario con estacionalidad | | SARIMA | $(p,d,q)(P,D,Q)_S$ | Caso general |

Cuándo usar diferenciación estacional

La diferencia estacional $(1-B^S)^D$ es necesaria cuando la ACF tiene picos significativos en los rezagos múltiplos de $S$ (12, 24, 36, ...) que no decaen. Esto indica una raíz unitaria estacional: el patrón estacional es no estacionario (cambia de nivel a lo largo del tiempo).

▶ SARIMA automático con auto.arima()(forecast)

Ctrl+Enter

auto.arima() con seasonal=TRUE evalúa automáticamente modelos SARIMA. La diferencia estacional diff(yt, 12) remueve el patrón estacional no estacionario.

11.6 Regresión con errores ARMA

El modelo de regresión con errores ARMA combina una componente determinística (tendencia + estacionalidad) con una estructura estocástica en los errores:

$Y_t = \mathbf{x}_t^\top \boldsymbol{\beta} + Z_t$

donde $\mathbf{x}_t$ son las variables explicativas (polinomios de tiempo, dummies estacionales, armónicos de Fourier, regresores externos) y $Z_t$ sigue un proceso ARMA:

$\phi(B) \, Z_t = \theta(B) \, a_t$

Esto es equivalente a un modelo SARIMAX con regresores externos. La función Arima() del paquete forecast lo implementa mediante el argumento xreg.

Por qué modelar los errores como ARMA

Cuando los residuos de una regresión por MCO están autocorrelacionados, los errores estándar de los coeficientes son incorrectos (subestimados), las pruebas t y F son inválidas, y los pronósticos son ineficientes. Modelar los errores como ARMA:

Corrige la inferencia (errores estándar válidos)
Mejora la eficiencia del pronóstico
Proporciona intervalos de predicción bien calibrados

Ejemplo completo: Salario real colombiano

La serie del índice de salario real de Colombia (enero 1990 – junio 2012, mensual) tiene tendencia creciente y estacionalidad. Se ajusta un modelo polinomial de grado 4 con dummies estacionales y se comparan cinco especificaciones de errores: sin modelar (MCO), AR(13), ARMA(2,0)(1,1)[12], ARMA(6,4), y ARMA(9,12) con restricciones.

Regresión con errores ARMA — Salario real Colombia

1library(forecast); library(lmtest)
2source("https://raw.githubusercontent.com/NelfiGonzalez/Funciones-de-Usuario-Estadistica-III/main/Funcion-Mipoly.R")
3source("https://raw.githubusercontent.com/NelfiGonzalez/Funciones-de-Usuario-Estadistica-III/main/Funciones-Criterios.Informacion-Calidad.Intervalos.R")
4source("https://raw.githubusercontent.com/NelfiGonzalez/Funciones-de-Usuario-Estadistica-III/main/Funciones-BP.LB.test-pruebaDW1.R")
5salario <- read.table("DATOSSALARIOREALENER1990-JUNIO2012.txt", header = TRUE)
6salario <- ts(salario, freq = 12, start = c(1990, 1))
7n <- length(salario) - 12; t <- 1:n
8yt <- ts(salario[1:n], freq = 12, start = c(1990, 1))
9mes  <- seasonaldummy(yt)
10poli <- Mipoly(tiempo = t, grado = 4)
11X    <- data.frame(poli, mes)
12modelo1 <- lm(yt ~ ., data = X)
13k1 <- length(coef(modelo1)[coef(modelo1) != 0])
14par(mfrow = c(1, 2))
15acf(as.numeric(residuals(modelo1)), ci.type = "ma", lag.max = 36, main = "ACF residuos MCO")
16pacf(as.numeric(residuals(modelo1)), lag.max = 36, main = "PACF residuos MCO")
17par(mfrow = c(1, 1))
18auto.arima(ts(residuals(modelo1), freq = 12), ic = "aic")
19modelo2 <- Arima(yt, order = c(13, 0, 0), xreg = as.matrix(X), method = "ML")
20modelo3 <- Arima(yt, order = c(2, 0, 0), seasonal = list(order = c(1, 0, 1), period = 12), xreg = as.matrix(X), method = "ML")
21modelo4 <- Arima(yt, order = c(6, 0, 4), xreg = as.matrix(X), method = "ML")
22criterios <- data.frame(
23Modelo = 1:4,
24AIC = c(AIC(modelo1), AIC(modelo2), AIC(modelo3), AIC(modelo4)),
25BIC = c(BIC(modelo1), BIC(modelo2), BIC(modelo3), BIC(modelo4)))
26print(criterios)
27for (i in 1:4) {
28m <- get(paste0("modelo", i))
29e <- residuals(m)
30lb <- Box.test(as.numeric(e), lag = 36, type = "Ljung-Box")
31cat(sprintf("Modelo %d: Ljung-Box Q=%.2f, p=%.4f → %s\n",
32            i, lb$statistic, lb$p.value,
33            ifelse(lb$p.value > 0.05, "ruido blanco ✓", "autocorrelación ✗")))
34}
35tnuevo  <- (n + 1):length(salario)
36Xnuevo  <- data.frame(Mipoly(tnuevo, 4), seasonaldummy(yt, h = 12))
37fc3 <- forecast(modelo3, xreg = as.matrix(Xnuevo), level = 95)
38plot(fc3, main = "Pronóstico salario real — Modelo 3: ARMA(2,0)(1,1)[12]")

El argumento xreg de Arima() incorpora regresores externos. Los errores se modelan como ARMA. La identificación del orden ARMA se hace sobre los residuos del MCO inicial.

11.7 SARIMAX: el modelo unificador

Combinando regresión externa con SARIMA, obtenemos el modelo SARIMAX completo:

$\phi(B) \, \Phi(B^S) \, (1-B)^d \, (1-B^S)^D \, (Y_t - \mathbf{x}_t^\top \boldsymbol{\beta}) = \theta(B) \, \Theta(B^S) \, a_t$

Este modelo engloba todos los casos particulares vistos en el curso:

| Modelo del curso | Especificación SARIMAX | Módulo | |---|---|---| | Regresión polinomial (MCO) | $(0,0,0)(0,0,0) + \mathbf{X}$ | 2–4 | | Regresión con errores AR(p) | $(p,0,0)(0,0,0) + \mathbf{X}$ | 11 | | Regresión con errores ARMA | $(p,0,q)(0,0,0) + \mathbf{X}$ | 11 | | Regresión con errores SARIMA | $(p,d,q)(P,D,Q)_S + \mathbf{X}$ | 11 | | ARIMA puro (sin regresores) | $(p,d,q)(P,D,Q)_S$ | 11 |

Conexión con TimeSight

TimeSight implementa el flujo SARIMAX completo en su paso de Modelado. Al seleccionar la familia ARIMA, puedes:

Usar auto.arima() para selección automática por AIC
Especificar manualmente los órdenes $(p,d,q)(P,D,Q)[S]$
El back-transform se aplica automáticamente si la serie fue transformada

Los pasos de Diagnóstico y Pronóstico verifican que los residuos sean ruido blanco y generan intervalos de predicción con el modelo SARIMA ajustado.

11.8 Resumen del procedimiento completo

Flujo de análisis SARIMAX

Paso 1 — Exploración: Graficar la serie, ACF, PACF, tests ADF/KPSS.
Paso 2 — Transformación: log o Box-Cox si la varianza no es constante.
Paso 3 — Diferenciación: Aplicar $\nabla^d$ y/o $\nabla_S^D$ hasta lograr estacionariedad.
Paso 4 — Identificación: ACF/PACF de la serie diferenciada → identificar $(p,q)$ y $(P,Q)$ . Confirmar con EACF y auto.arima.
Paso 5 — Estimación: Ajustar SARIMA por ML con Arima(). Si hay regresores, usar xreg.
Paso 6 — Diagnóstico: ACF/PACF de residuos, Ljung-Box, Shapiro-Wilk. Si no son ruido blanco, volver al paso 4.
Paso 7 — Pronóstico: forecast() con horizonte $h$ y nivel de confianza.
Paso 8 — Validación: Comparar con datos held-out. Calcular RMSE, MAPE, cobertura de IP.

▶ Flujo SARIMAX completo — AirPassengers(forecast)

library(forecast)
yt <- AirPassengers
# Paso 1: Exploración
par(mfrow = c(2, 2))
plot(yt, main = "AirPassengers")
acf(yt,  lag.max = 48, main = "ACF")
pacf(yt, lag.max = 48, main = "PACF")
plot(log(yt), main = "log(AirPassengers)")
par(mfrow = c(1, 1))
# Paso 2: Transformación log
yt_log <- log(yt)
# Paso 3: Diferenciación
yt_d <- diff(diff(yt_log, 12))
par(mfrow = c(1, 2))
acf(yt_d,  lag.max = 48, main = "ACF diff(diff(log(Y),12))")
pacf(yt_d, lag.max = 48, main = "PACF diff(diff(log(Y),12))")
par(mfrow = c(1, 1))
# Paso 4-5: Identificación + Estimación
m <- Arima(yt_log, order = c(0, 1, 1),
         seasonal = list(order = c(0, 1, 1), period = 12), method = "ML")
summary(m)
# Paso 6: Diagnóstico
e <- residuals(m)
par(mfrow = c(1, 3))
acf(e, lag.max = 48, main = "ACF residuos")
pacf(e, lag.max = 48, main = "PACF residuos")
qqnorm(e); qqline(e, col = 2)
par(mfrow = c(1, 1))
cat("Ljung-Box:", Box.test(e, lag = 36, type = "Ljung-Box")$p.value, "\n")
cat("Shapiro:",  shapiro.test(e)$p.value, "\n")
# Paso 7: Pronóstico
fc <- forecast(m, h = 24, level = c(80, 95))
fc_orig <- exp(fc$mean)
fc_lo80 <- exp(fc$lower[, 2])
fc_hi80 <- exp(fc$upper[, 2])
fc_lo95 <- exp(fc$lower[, 1])
fc_hi95 <- exp(fc$upper[, 1])
plot(fc, main = "Pronóstico SARIMA(0,1,1)(0,1,1)[12] — log(AirPassengers)")

Ctrl+Enter

Este es el flujo que TimeSight automatiza: exploración → transformación → modelado SARIMA → diagnóstico → pronóstico con back-transform.