Análisis Descriptivo Serie ISE Colombia
library(readxl)
library(tidyverse)
library(zoo)
library(xts)
library(TSstudio)
library(car)
library(forecast)
library(astsa)
library(nonlinearTseries)
library(energy)
library(biwavelet)
library(knitr)1. Gráfico de la serie de tiempo
Índice de Seguimiento a la Economía (ISE) de Colombia, mensual, de enero de 2005 a junio de 2026. La serie tiene frecuencia de medición 12: cada 12 periodos se vuelve a ver el mismo mes. Las características a revisar son tendencia, estacionalidad, ciclos y varianza marginal no constante.
data <- read_excel("ISE2005_2026.xlsx")
str(data)## tibble [258 × 2] (S3: tbl_df/tbl/data.frame)
## $ Fecha: POSIXct[1:258], format: "2005-01-01" "2005-02-01" ...
## $ Iset : num [1:258] 60 61.5 62.2 63 62.7 ...head(data)## # A tibble: 6 × 2
## Fecha Iset
## <dttm> <dbl>
## 1 2005-01-01 00:00:00 60.0
## 2 2005-02-01 00:00:00 61.5
## 3 2005-03-01 00:00:00 62.2
## 4 2005-04-01 00:00:00 63.0
## 5 2005-05-01 00:00:00 62.7
## 6 2005-06-01 00:00:00 63.5tail(data)## # A tibble: 6 × 2
## Fecha Iset
## <dttm> <dbl>
## 1 2026-01-01 00:00:00 118.
## 2 2026-02-01 00:00:00 120.
## 3 2026-03-01 00:00:00 129.
## 4 2026-04-01 00:00:00 124.
## 5 2026-05-01 00:00:00 129.
## 6 2026-06-01 00:00:00 127.fechas = as.yearmon(data$Fecha)
ise_xts = xts(x = data$Iset, frequency = 12, order.by = fechas)
ise_ts = ts(data$Iset, start = c(2005, 1), frequency = 12)
plot(ise_xts, main = "ISE Colombia")
ts_plot(ise_xts,
title = "Índice de Seguimiento a la Economía (ISE) - Colombia",
Ytitle = "Índice",
Xtitle = "Año",
Xgrid = TRUE,
Ygrid = TRUE)2. Análisis de varianza marginal: transformación Box-Cox
En ocasiones la serie presenta varianza marginal no constante a lo largo del tiempo, lo cual hace necesario tener en cuenta tal característica. En este caso se sugiere hacer una transformación de potencia para estabilizar la varianza. Esta familia de transformaciones se llama transformaciones Box-Cox:
\[ f_{\lambda}(u_{t})= \begin{cases} \lambda^{-1}(u^{\lambda}_{t}-1), & \text{si } \lambda \neq 0,\\ \ln(u_{t}), & \text{si } \lambda=0, \end{cases} \qquad u_t>0. \]
El parámetro \(\lambda\) se estima
por dos vías. Por máxima verosimilitud
(car::powerTransform), que elige el \(\lambda\) que maximiza la log-verosimilitud
perfilada \(\hat{\lambda}=\arg\max_\lambda
\ell(\lambda)\) suponiendo normalidad de los datos transformados.
Por el método de Guerrero (forecast::BoxCox.lambda), que
divide la serie en subseries de un año y escoge el \(\lambda\) que hace más constante el
cociente \(s_i/\bar{x}_i^{\,1-\lambda}\) entre ellas,
es decir, que desliga la variabilidad del nivel. Valores de referencia:
\(\lambda = 1\) no transforma, \(\lambda = 0.5\) es tipo raíz cuadrada y
\(\lambda = 0\) es el logaritmo.
Cálculo del lambda de Guerrero
El método de Guerrero (1993) tiene la ventaja de no requerir ningún supuesto distribucional sobre la serie (a diferencia de la máxima verosimilitud, que supone normalidad de los datos transformados), por lo que resulta más robusto cuando ese supuesto no se cumple o hay datos atípicos. La idea es la siguiente: la serie se divide en \(n\) subperiodos no traslapados de longitud igual al periodo estacional (12 meses, es decir, un año cada uno), y para cada subperiodo \(i=1,\dots,n\) se calculan su media \(\bar{x}_i\) y su desviación estándar \(s_i\). Si la varianza marginal no es constante, estas desviaciones estándar tienden a crecer o decrecer junto con el nivel de la serie, de modo que \(s_i\) y \(\bar{x}_i\) quedan relacionadas entre los distintos subperiodos. Para cada valor de \(\lambda\) se construye el cociente \(s_i/\bar{x}_i^{\,1-\lambda}\) en cada uno de los \(n\) subperiodos, y luego se calcula el coeficiente de variación de ese conjunto de \(n\) cocientes, es decir, la desviación estándar de esos \(n\) valores dividida entre su media. El lambda de Guerrero es precisamente el que minimiza ese coeficiente de variación: el que hace más estable, entre todos los subperiodos, la relación entre las desviaciones estándar y las medias de la serie: \[\hat{\lambda}_{Guerrero}=\arg\min_{\lambda}\;\frac{\operatorname{sd}_i\!\left(s_i/\bar{x}_i^{\,1-\lambda}\right)}{\operatorname{mean}_i\!\left(s_i/\bar{x}_i^{\,1-\lambda}\right)}.\]
# Perfil de log-verosimilitud: el máximo es el lambda estimado y las líneas punteadas el IC del 95 %
car::boxCox(ise_ts ~ 1)
# Máxima verosimilitud
lambda_verosim_v = car::powerTransform(ise_ts ~ 1)
summary(lambda_verosim_v)## bcPower Transformation to Normality
## Est Power Rounded Pwr Wald Lwr Bnd Wald Upr Bnd
## Y1 0.8647 1 0.2356 1.4939
##
## Likelihood ratio test that transformation parameter is equal to 0
## (log transformation)
## LRT df pval
## LR test, lambda = (0) 7.32938 1 0.0067836
##
## Likelihood ratio test that no transformation is needed
## LRT df pval
## LR test, lambda = (1) 0.1772189 1 0.67377lambda_verosim = as.numeric(lambda_verosim_v$lambda)
# Guerrero
lambda_guerrero = forecast::BoxCox.lambda(ise_ts, method = "guerrero", lower = -1, upper = 3)
sprintf("Lambda máxima verosimilitud: %.4f", lambda_verosim)## [1] "Lambda máxima verosimilitud: 0.8647"sprintf("Lambda Guerrero: %.4f", lambda_guerrero)## [1] "Lambda Guerrero: -0.5582"par(mfrow = c(2, 2), mar = c(3, 3, 3, 1))
plot(ise_ts, main = "ISE sin transformar")
plot(forecast::BoxCox(ise_ts, lambda = lambda_verosim),
main = sprintf("ISE Box-Cox máx. verosimilitud (lambda = %.3f)", lambda_verosim))
plot(forecast::BoxCox(ise_ts, lambda = lambda_guerrero),
main = sprintf("ISE Box-Cox Guerrero (lambda = %.3f)", lambda_guerrero))
plot(log(ise_ts), main = "ISE logaritmo (lambda = 0)")
par(mfrow = c(1, 1))# DECISIÓN: el análisis de tendencia se hace para la serie sin transformar y para la transformada por Guerrero
lambda_final = lambda_guerrero
ise_boxcox = forecast::BoxCox(ise_ts, lambda = lambda_final)
plot(ise_boxcox, main = sprintf("ISE Box-Cox Guerrero (lambda = %.3f)", lambda_final))
3. Análisis de tendencia
El modelo inicial para el análisis de tendencia es \[y_t=\mu_t+a_t,\] donde \(\mu_t\) es la tendencia y \(a_t\) la componente aleatoria. Como \(\mu_t\) no es observable se estima con distintos métodos, y la serie sin tendencia se obtiene como \(\hat{a}_t = y_t - \hat{\mu}_t\). Cada método se aplica a la serie original y a la serie transformada por Box-Cox (Guerrero).
i) Descomposición clásica (decompose)
La descomposición clásica representa la serie como \[y_t=\mu_t+s_t+a_t,\] donde \(s_t\) es la componente estacional. Primero estima la tendencia \(\mu_t\) con un filtro de promedios móviles centrado de orden 12 (el mismo filtro 2×12 de la siguiente sección), luego la componente estacional promediando cada mes sobre la serie sin tendencia, y lo que sobra es la componente aleatoria.
decomp_original = decompose(ise_ts)
plot(decomp_original)
fit_tendencia_original_decompose = decomp_original$trend
ise_original_decompose = ise_ts - fit_tendencia_original_decompose
par(mfrow = c(2, 1), mar = c(3, 3, 3, 1))
plot(ise_ts, main = "ISE original - Tendencia decompose")
lines(fit_tendencia_original_decompose, lwd = 2, col = 4)
plot(ise_original_decompose, main = "ISE original sin tendencia - decompose")
abline(h = 0, lty = 2)
par(mfrow = c(1, 1))decomp_boxcox = decompose(ise_boxcox)
plot(decomp_boxcox)
fit_tendencia_boxcox_decompose = decomp_boxcox$trend
ise_boxcox_decompose = ise_boxcox - fit_tendencia_boxcox_decompose
par(mfrow = c(2, 1), mar = c(3, 3, 3, 1))
plot(ise_boxcox, main = "ISE Box-Cox - Tendencia decompose")
lines(fit_tendencia_boxcox_decompose, lwd = 2, col = 4)
plot(ise_boxcox_decompose, main = "ISE Box-Cox sin tendencia - decompose")
abline(h = 0, lty = 2)
par(mfrow = c(1, 1))ii) Promedio móvil (filtro Boxcar 2×12)
El promedio móvil es útil para descubrir rasgos de una serie como tendencias de largo plazo y componentes estacionales. Si \(x_t\) representa las observaciones, una forma de estimar la tendencia es \[m_t=\sum_{j=-k}^{k}a_jx_{t-j},\] donde, si \(a_j=a_{-j}\ge0\) y \(\sum_{j=-k}^{k}a_j=1\), se conoce como el promedio móvil simétrico de los datos. Con periodo par (12) se usan 13 pesos: \(1/24\) en los extremos y \(1/12\) en los 11 del centro.
# Parámetro importante: número de pesos. Debe cubrir exactamente un año (12 meses, o múltiplos)
# para que cada mes pese igual y la estacionalidad se cancele dentro del promedio.
wgts = c(.5, rep(1, 11), .5) / 12
fit_tendencia_original_ma = stats::filter(ise_ts, sides = 2, filter = wgts)
ise_original_ma = ise_ts - fit_tendencia_original_ma
par(mfrow = c(2, 1), mar = c(3, 3, 3, 1))
plot(ise_ts, main = "ISE original - Tendencia promedio móvil 2x12")
lines(fit_tendencia_original_ma, lwd = 2, col = 4)
plot(ise_original_ma, main = "ISE original sin tendencia - promedio móvil 2x12")
abline(h = 0, lty = 2)
par(mfrow = c(1, 1))fit_tendencia_boxcox_ma = stats::filter(ise_boxcox, sides = 2, filter = wgts)
ise_boxcox_ma = ise_boxcox - fit_tendencia_boxcox_ma
par(mfrow = c(2, 1), mar = c(3, 3, 3, 1))
plot(ise_boxcox, main = "ISE Box-Cox - Tendencia promedio móvil 2x12")
lines(fit_tendencia_boxcox_ma, lwd = 2, col = 4)
plot(ise_boxcox_ma, main = "ISE Box-Cox sin tendencia - promedio móvil 2x12")
abline(h = 0, lty = 2)
par(mfrow = c(1, 1))
# Forma del filtro
nwgts = c(rep(0, 20), wgts, rep(0, 20))
plot(nwgts, type = "l", ylim = c(-.02, .1), xaxt = "n", yaxt = "n", ann = FALSE)
iii) Suavizamiento Kernel (12, 24 y 36 meses)
El suavizamiento kernel es un suavizador de promedio móvil que
utiliza una función de ponderación, o kernel, para promediar las
observaciones: \[m_t=\sum_{i=1}^{n}w_i(t)x_i,\] donde \[w_i(t)=K\left(\frac{t-i}{b}\right)\Big/\sum_{j=1}^{n}K\left(\frac{t-j}{b}\right)\]
son los pesos, \(K(\cdot)\) es una
función kernel y \(b\) el ancho de
banda. Este estimador se llama estimador de Nadaraya-Watson. Se usa un
kernel normal con la función ksmooth.
# Parámetro importante: bandwidth. Como time(ise_ts) está en años, bandwidth = 1 equivale a una
# ventana de 12 meses, 2 a 24 meses y 3 a 36 meses. Más ancho = tendencia más suave.
fit_tendencia_original_kernel12 = ts(ksmooth(time(ise_ts), ise_ts, "normal", bandwidth = 1, x.points = time(ise_ts))$y,
start = start(ise_ts), frequency = 12)
fit_tendencia_original_kernel24 = ts(ksmooth(time(ise_ts), ise_ts, "normal", bandwidth = 2, x.points = time(ise_ts))$y,
start = start(ise_ts), frequency = 12)
fit_tendencia_original_kernel36 = ts(ksmooth(time(ise_ts), ise_ts, "normal", bandwidth = 3, x.points = time(ise_ts))$y,
start = start(ise_ts), frequency = 12)
ise_original_kernel12 = ise_ts - fit_tendencia_original_kernel12
ise_original_kernel24 = ise_ts - fit_tendencia_original_kernel24
ise_original_kernel36 = ise_ts - fit_tendencia_original_kernel36
par(mfrow = c(3, 2), mar = c(3, 3, 3, 1))
plot(ise_ts, main = "ISE original - Tendencia Kernel 12 meses")
lines(fit_tendencia_original_kernel12, lwd = 2, col = 4)
plot(ise_original_kernel12, main = "ISE original sin tendencia - Kernel 12")
abline(h = 0, lty = 2)
plot(ise_ts, main = "ISE original - Tendencia Kernel 24 meses")
lines(fit_tendencia_original_kernel24, lwd = 2, col = 4)
plot(ise_original_kernel24, main = "ISE original sin tendencia - Kernel 24")
abline(h = 0, lty = 2)
plot(ise_ts, main = "ISE original - Tendencia Kernel 36 meses")
lines(fit_tendencia_original_kernel36, lwd = 2, col = 4)
plot(ise_original_kernel36, main = "ISE original sin tendencia - Kernel 36")
abline(h = 0, lty = 2)
par(mfrow = c(1, 1))fit_tendencia_boxcox_kernel12 = ts(ksmooth(time(ise_boxcox), ise_boxcox, "normal", bandwidth = 1,
x.points = time(ise_boxcox))$y,
start = start(ise_boxcox), frequency = 12)
fit_tendencia_boxcox_kernel24 = ts(ksmooth(time(ise_boxcox), ise_boxcox, "normal", bandwidth = 2,
x.points = time(ise_boxcox))$y,
start = start(ise_boxcox), frequency = 12)
fit_tendencia_boxcox_kernel36 = ts(ksmooth(time(ise_boxcox), ise_boxcox, "normal", bandwidth = 3,
x.points = time(ise_boxcox))$y,
start = start(ise_boxcox), frequency = 12)
ise_boxcox_kernel12 = ise_boxcox - fit_tendencia_boxcox_kernel12
ise_boxcox_kernel24 = ise_boxcox - fit_tendencia_boxcox_kernel24
ise_boxcox_kernel36 = ise_boxcox - fit_tendencia_boxcox_kernel36
par(mfrow = c(3, 2), mar = c(3, 3, 3, 1))
plot(ise_boxcox, main = "ISE Box-Cox - Tendencia Kernel 12 meses")
lines(fit_tendencia_boxcox_kernel12, lwd = 2, col = 4)
plot(ise_boxcox_kernel12, main = "ISE Box-Cox sin tendencia - Kernel 12")
abline(h = 0, lty = 2)
plot(ise_boxcox, main = "ISE Box-Cox - Tendencia Kernel 24 meses")
lines(fit_tendencia_boxcox_kernel24, lwd = 2, col = 4)
plot(ise_boxcox_kernel24, main = "ISE Box-Cox sin tendencia - Kernel 24")
abline(h = 0, lty = 2)
plot(ise_boxcox, main = "ISE Box-Cox - Tendencia Kernel 36 meses")
lines(fit_tendencia_boxcox_kernel36, lwd = 2, col = 4)
plot(ise_boxcox_kernel36, main = "ISE Box-Cox sin tendencia - Kernel 36")
abline(h = 0, lty = 2)
par(mfrow = c(1, 1))# Forma del kernel normal
gauss = function(x) { 1/sqrt(2*pi) * exp(-(x^2)/2) }
x = seq(from = -3, to = 3, by = 0.001)
plot(x, gauss(x), type = "l", ylim = c(-.02, .45), xaxt = "n", yaxt = "n", ann = FALSE)
iv) Lowess: vecinos más cercanos (12, 24 y 36)
Otro enfoque para suavizar un gráfico de tiempo es la regresión del
vecino más cercano. La técnica se basa en la regresión de \(k\) vecinos más cercanos, en la que uno usa
solo los datos \(\{x_{t-k/2}, \dots, x_t,
\dots, x_{t+k/2}\}\) para predecir \(x_t\) mediante regresión del tiempo, y
luego establece \(m_t = \hat{x}_t\).
Primero, una cierta proporción de vecinos más cercanos a \(x_t\) para el tiempo \(t\) se incluyen en un esquema de
ponderación (\(\nu_t =
W\left(\frac{|t_i-t|}{\lambda_q(t)}\right)\)); los valores más
cercanos a \(x_t\) en el tiempo
obtienen más peso. Luego se utiliza una regresión (polinomial de grado
\(d\), usualmente 1 o 2, es decir un
ajuste localmente lineal o cuadrático) ponderada robusta entre \(t\) y \(x_t\) para predecir \(x_t\) y obtener los valores suavizados
\(m_t\). Cuanto mayor sea la fracción
de vecinos más cercanos incluidos, más suave será el ajuste. En R se usa
la función loess. La función \(W(\cdot)\) es la función de peso tricúbica
\[
W(u)= \begin{cases}
(1-u^3)^3, & \text{para } 0 \le u < 1,\\
0, & \text{para } u \ge 1.
\end{cases}
\]
# Parámetro importante: f = fracción de observaciones en cada regresión local (vecinos / n).
# Más vecinos = tendencia más suave.
n = length(ise_ts)
f_12 = 12/n
f_24 = 24/n
f_36 = 36/n
fit_tendencia_original_lowess12 = ts(lowess(ise_ts, f = f_12)$y, start = start(ise_ts), frequency = 12)
fit_tendencia_original_lowess24 = ts(lowess(ise_ts, f = f_24)$y, start = start(ise_ts), frequency = 12)
fit_tendencia_original_lowess36 = ts(lowess(ise_ts, f = f_36)$y, start = start(ise_ts), frequency = 12)
ise_original_lowess12 = ise_ts - fit_tendencia_original_lowess12
ise_original_lowess24 = ise_ts - fit_tendencia_original_lowess24
ise_original_lowess36 = ise_ts - fit_tendencia_original_lowess36
par(mfrow = c(3, 2), mar = c(3, 3, 3, 1))
plot(ise_ts, main = "ISE original - Tendencia Lowess 12 vecinos")
lines(fit_tendencia_original_lowess12, lwd = 2, col = 4)
plot(ise_original_lowess12, main = "ISE original sin tendencia - Lowess 12")
abline(h = 0, lty = 2)
plot(ise_ts, main = "ISE original - Tendencia Lowess 24 vecinos")
lines(fit_tendencia_original_lowess24, lwd = 2, col = 4)
plot(ise_original_lowess24, main = "ISE original sin tendencia - Lowess 24")
abline(h = 0, lty = 2)
plot(ise_ts, main = "ISE original - Tendencia Lowess 36 vecinos")
lines(fit_tendencia_original_lowess36, lwd = 2, col = 4)
plot(ise_original_lowess36, main = "ISE original sin tendencia - Lowess 36")
abline(h = 0, lty = 2)
par(mfrow = c(1, 1))fit_tendencia_boxcox_lowess12 = ts(lowess(ise_boxcox, f = f_12)$y, start = start(ise_boxcox), frequency = 12)
fit_tendencia_boxcox_lowess24 = ts(lowess(ise_boxcox, f = f_24)$y, start = start(ise_boxcox), frequency = 12)
fit_tendencia_boxcox_lowess36 = ts(lowess(ise_boxcox, f = f_36)$y, start = start(ise_boxcox), frequency = 12)
ise_boxcox_lowess12 = ise_boxcox - fit_tendencia_boxcox_lowess12
ise_boxcox_lowess24 = ise_boxcox - fit_tendencia_boxcox_lowess24
ise_boxcox_lowess36 = ise_boxcox - fit_tendencia_boxcox_lowess36
par(mfrow = c(3, 2), mar = c(3, 3, 3, 1))
plot(ise_boxcox, main = "ISE Box-Cox - Tendencia Lowess 12 vecinos")
lines(fit_tendencia_boxcox_lowess12, lwd = 2, col = 4)
plot(ise_boxcox_lowess12, main = "ISE Box-Cox sin tendencia - Lowess 12")
abline(h = 0, lty = 2)
plot(ise_boxcox, main = "ISE Box-Cox - Tendencia Lowess 24 vecinos")
lines(fit_tendencia_boxcox_lowess24, lwd = 2, col = 4)
plot(ise_boxcox_lowess24, main = "ISE Box-Cox sin tendencia - Lowess 24")
abline(h = 0, lty = 2)
plot(ise_boxcox, main = "ISE Box-Cox - Tendencia Lowess 36 vecinos")
lines(fit_tendencia_boxcox_lowess36, lwd = 2, col = 4)
plot(ise_boxcox_lowess36, main = "ISE Box-Cox sin tendencia - Lowess 36")
abline(h = 0, lty = 2)
par(mfrow = c(1, 1))v) Tendencia determinística cúbica
Si la tendencia fuera lineal, el modelo sería \(y_t=\beta_0+\beta_1 t +a_t\). Como la tendencia del ISE no se ve lineal, se ajusta un polinomio de grado 3 por mínimos cuadrados: \[y_t=\beta_0+\beta_1 t+\beta_2 t^2+\beta_3 t^3+a_t.\]
# Parámetro importante: grado del polinomio (aquí 3)
t_ise = 1:length(ise_ts)
lm_original_cubica = lm(ise_ts ~ t_ise + I(t_ise^2) + I(t_ise^3))
fit_tendencia_original_cubica = ts(fitted(lm_original_cubica), start = start(ise_ts), frequency = 12)
ise_original_cubica = ise_ts - fit_tendencia_original_cubica
par(mfrow = c(2, 1), mar = c(3, 3, 3, 1))
plot(ise_ts, main = "ISE original - Tendencia cúbica")
lines(fit_tendencia_original_cubica, lwd = 2, col = 2)
plot(ise_original_cubica, main = "ISE original sin tendencia - cúbica")
abline(h = 0, lty = 2)
par(mfrow = c(1, 1))lm_boxcox_cubica = lm(ise_boxcox ~ t_ise + I(t_ise^2) + I(t_ise^3))
fit_tendencia_boxcox_cubica = ts(fitted(lm_boxcox_cubica), start = start(ise_boxcox), frequency = 12)
ise_boxcox_cubica = ise_boxcox - fit_tendencia_boxcox_cubica
par(mfrow = c(2, 1), mar = c(3, 3, 3, 1))
plot(ise_boxcox, main = "ISE Box-Cox - Tendencia cúbica")
lines(fit_tendencia_boxcox_cubica, lwd = 2, col = 2)
plot(ise_boxcox_cubica, main = "ISE Box-Cox sin tendencia - cúbica")
abline(h = 0, lty = 2)
par(mfrow = c(1, 1))vi) Tendencia determinística lineal
Se ajusta la tendencia más simple, una recta en el tiempo, por mínimos cuadrados: \[y_t=\beta_0+\beta_1 t +a_t.\]
lm_original_lineal = lm(ise_ts ~ t_ise)
summary(lm_original_lineal)##
## Call:
## lm(formula = ise_ts ~ t_ise)
##
## Residuals:
## Min 1Q Median 3Q Max
## -27.3748 -2.9384 0.1219 2.6005 16.8571
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 65.007925 0.697954 93.14 <2e-16 ***
## t_ise 0.247204 0.004672 52.91 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 5.589 on 256 degrees of freedom
## Multiple R-squared: 0.9162, Adjusted R-squared: 0.9159
## F-statistic: 2800 on 1 and 256 DF, p-value: < 2.2e-16fit_tendencia_original_lineal = ts(fitted(lm_original_lineal), start = start(ise_ts), frequency = 12)
ise_original_lineal = ise_ts - fit_tendencia_original_lineal
par(mfrow = c(2, 1), mar = c(3, 3, 3, 1))
plot(ise_ts, main = "ISE original - Tendencia lineal")
lines(fit_tendencia_original_lineal, lwd = 2, col = 2)
plot(ise_original_lineal, main = "ISE original sin tendencia - lineal")
abline(h = 0, lty = 2)
par(mfrow = c(1, 1))lm_boxcox_lineal = lm(ise_boxcox ~ t_ise)
summary(lm_boxcox_lineal)##
## Call:
## lm(formula = ise_boxcox ~ t_ise)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0214972 -0.0034461 0.0000991 0.0035711 0.0134376
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 1.622e+00 6.603e-04 2456.93 <2e-16 ***
## t_ise 2.106e-04 4.420e-06 47.64 <2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.005287 on 256 degrees of freedom
## Multiple R-squared: 0.8986, Adjusted R-squared: 0.8983
## F-statistic: 2270 on 1 and 256 DF, p-value: < 2.2e-16fit_tendencia_boxcox_lineal = ts(fitted(lm_boxcox_lineal), start = start(ise_boxcox), frequency = 12)
ise_boxcox_lineal = ise_boxcox - fit_tendencia_boxcox_lineal
par(mfrow = c(2, 1), mar = c(3, 3, 3, 1))
plot(ise_boxcox, main = "ISE Box-Cox - Tendencia lineal")
lines(fit_tendencia_boxcox_lineal, lwd = 2, col = 2)
plot(ise_boxcox_lineal, main = "ISE Box-Cox sin tendencia - lineal")
abline(h = 0, lty = 2)
par(mfrow = c(1, 1))vii) Suavizamiento por splines
Una forma obvia de suavizar los datos sería ajustar una regresión polinomial en términos del tiempo. Por ejemplo, un polinomio cúbico tendría \(x_t = m_t + w_t\), donde \(m_t = \beta_0+\beta_1 t+\beta_2 t^2+\beta_3 t^3\). Entonces podríamos ajustar \(m_t\) mediante mínimos cuadrados ordinarios.
Una extensión de la regresión polinomial es dividir primero el tiempo
\(t=1,\dots,n\) en \(k\) intervalos, \([t_0=1,t_1],[t_1+1,t_2],\dots,[t_{k-1}+1,t_k=n]\);
los valores \(t_0,t_1,\dots,t_k\) se
llaman nodos. Luego, en cada intervalo se ajusta una regresión
polinomial, normalmente de orden 3, y esto se llama splines cúbicos. Un
método relacionado es suavizar splines, que minimiza el compromiso entre
el ajuste y el grado de suavidad dado por \[\sum_{t=1}^{n}[x_t-m_t]^2+\lambda\int
(m_t'')^2\,dt,\] donde \(m_t\) es un spline cúbico con nodos en cada
tiempo \(t\) y el grado de suavidad es
controlado por \(\lambda>0\). El
parámetro de suavizado en R está controlado por el argumento
spar de la función smooth.spline del paquete
stats, que si no se fija explícitamente se escoge por
validación cruzada.
# Parámetro importante: spar (parámetro de suavizado, 0-1). Si se omite, se elige por validación cruzada;
# más spar = tendencia más suave.
sm_original_spline = smooth.spline(as.numeric(time(ise_ts)), ise_ts)
fit_tendencia_original_spline = ts(sm_original_spline$y, start = start(ise_ts), frequency = 12)
ise_original_spline = ise_ts - fit_tendencia_original_spline
par(mfrow = c(2, 1), mar = c(3, 3, 3, 1))
plot(ise_ts, main = "ISE original - Tendencia spline")
lines(fit_tendencia_original_spline, lwd = 2, col = 2)
plot(ise_original_spline, main = "ISE original sin tendencia - spline")
abline(h = 0, lty = 2)
par(mfrow = c(1, 1))sm_boxcox_spline = smooth.spline(as.numeric(time(ise_boxcox)), ise_boxcox)
fit_tendencia_boxcox_spline = ts(sm_boxcox_spline$y, start = start(ise_boxcox), frequency = 12)
ise_boxcox_spline = ise_boxcox - fit_tendencia_boxcox_spline
par(mfrow = c(2, 1), mar = c(3, 3, 3, 1))
plot(ise_boxcox, main = "ISE Box-Cox - Tendencia spline")
lines(fit_tendencia_boxcox_spline, lwd = 2, col = 2)
plot(ise_boxcox_spline, main = "ISE Box-Cox sin tendencia - spline")
abline(h = 0, lty = 2)
par(mfrow = c(1, 1))viii) Eliminación de la tendencia por primeras diferencias
La diferenciación no estima la tendencia, la elimina. La diferencia ordinaria de orden 1 es \[\nabla^1 Y_t=(1-B)^1 Y_t=Y_t-Y_{t-1},\] donde \(B\) es el operador de rezago (\(BY_t = Y_{t-1}\)). La serie resultante es la variación mes a mes y tiene una observación menos que la original.
# Parámetro importante: differences (orden d de la diferencia). Con d = 1 se elimina una tendencia lineal.
ise_original_diff = diff(ise_ts, lag = 1, differences = 1)
par(mfrow = c(2, 1), mar = c(3, 3, 3, 1))
plot(ise_ts, main = "ISE original")
plot(ise_original_diff, main = "ISE original - Primera diferencia")
abline(h = 0, lty = 2)
par(mfrow = c(1, 1))ise_boxcox_diff = diff(ise_boxcox, lag = 1, differences = 1)
par(mfrow = c(2, 1), mar = c(3, 3, 3, 1))
plot(ise_boxcox, main = "ISE Box-Cox")
plot(ise_boxcox_diff, main = "ISE Box-Cox - Primera diferencia")
abline(h = 0, lty = 2)
par(mfrow = c(1, 1))4. Selección de la tendencia
i) Serie sin transformar: tendencia Kernel 12 meses
plot(ise_original_kernel12, main = "ISE original sin tendencia - Kernel 12 meses", ylab = "ISE - tendencia")
abline(h = 0, lty = 2)
Al quitar la tendencia a la serie sin transformar, la amplitud de las oscilaciones crece claramente con el tiempo: la desviación estándar anual pasa de alrededor de 2.6 en 2005 a más de 5 desde 2022. Esa variabilidad que aumenta con el nivel de la serie es varianza marginal no constante, por lo que se justifica usar la serie transformada por Box-Cox, cuya desviación anual se mantiene en un rango estable.
ii) Serie transformada (Guerrero): comparación de métodos
par(mfrow = c(2, 1), mar = c(3, 4, 3, 1))
plot(ise_boxcox, col = "grey50", main = "ISE Box-Cox - Tendencias MA, Kernel 12 y Lowess 12", ylab = "ISE Box-Cox")
lines(fit_tendencia_boxcox_ma, lwd = 2, col = "red")
lines(fit_tendencia_boxcox_kernel12, lwd = 2, col = "blue")
lines(fit_tendencia_boxcox_lowess12, lwd = 2, col = "darkgreen")
legend("topleft", legend = c("MA 2x12", "Kernel 12", "Lowess 12"),
col = c("red", "blue", "darkgreen"), lwd = 2, cex = 0.7)
plot(ise_boxcox, col = "grey50", xlim = c(2019, 2022.5),
ylim = range(window(ise_boxcox, start = c(2019, 1), end = c(2022, 6))),
main = "ISE Box-Cox - Zoom COVID (2019-2022)", ylab = "ISE Box-Cox")
lines(fit_tendencia_boxcox_ma, lwd = 2, col = "red")
lines(fit_tendencia_boxcox_kernel12, lwd = 2, col = "blue")
lines(fit_tendencia_boxcox_lowess12, lwd = 2, col = "darkgreen")
legend("topleft", legend = c("MA 2x12", "Kernel 12", "Lowess 12"),
col = c("red", "blue", "darkgreen"), lwd = 2, cex = 0.7)
par(mfrow = c(1, 1))# Misma escala vertical en todos los paneles para poder compararlos
lim_y = range(c(ise_boxcox_ma, ise_boxcox_kernel12, ise_boxcox_lowess12), na.rm = TRUE)
par(mfrow = c(3, 1), mar = c(3, 3, 3, 1))
plot(ise_boxcox_ma, ylim = lim_y, main = "ISE Box-Cox sin tendencia - MA 2x12", col = "red")
abline(h = 0, lty = 2)
plot(ise_boxcox_kernel12, ylim = lim_y, main = "ISE Box-Cox sin tendencia - Kernel 12", col = "blue")
abline(h = 0, lty = 2)
plot(ise_boxcox_lowess12, ylim = lim_y, main = "ISE Box-Cox sin tendencia - Lowess 12", col = "darkgreen")
abline(h = 0, lty = 2)
par(mfrow = c(1, 1))iii) Tendencia seleccionada: Kernel 12 meses
# DECISIÓN: tendencia final = Kernel 12 meses sobre la serie Box-Cox (Guerrero)
fit_tendencia = fit_tendencia_boxcox_kernel12
ise_sin_tendencia = ise_boxcox_kernel12
par(mfrow = c(2, 1), mar = c(3, 3, 3, 1))
plot(ise_boxcox, main = "ISE Box-Cox - Tendencia seleccionada: Kernel 12 meses")
lines(fit_tendencia, lwd = 2, col = 4)
plot(ise_sin_tendencia, main = "ISE Box-Cox sin tendencia - Kernel 12 meses")
abline(h = 0, lty = 2)
par(mfrow = c(1, 1))5. Correlación y dependencia
Desde aquí se trabaja solo con la serie transformada por Box-Cox
(Guerrero) sin tendencia (Kernel 12 meses),
ise_sin_tendencia, porque con tendencia las
autocorrelaciones quedan dominadas por ella.
i) Gráficos de dispersión de retardos
Se hacen gráficos de dispersión para chequear qué tipos de relaciones hay entre los retardos de la variable de interés: la serie en \(t\) contra la serie en \(t-h\), para \(h = 1, \dots, 12\). Una nube alineada indica relación lineal (positiva o negativa) en ese retardo; una curva indica relación no lineal. Como la serie es mensual, \(h = 12\) compara observaciones separadas por un año.
par(mar = c(2, 2, 2, 2))
plot(ise_sin_tendencia, main = "ISE Box-Cox sin tendencia")
abline(h = 0, lty = 2)
par(mar = c(3, 2, 3, 2))
astsa::lag1.plot(ise_sin_tendencia, 12, corr = F)
ii) Función de autocorrelación simple (ACF)
Cuando el proceso es estacionario, o al menos no presenta tendencia, se puede usar el gráfico ACF para explorar las posibles relaciones lineales a diferentes rezagos: \[\hat{\rho}(k)=\frac{\sum_{t=k+1}^{T}(y_t-\bar{y})(y_{t-k}-\bar{y})}{\sum_{t=1}^{T}(y_t-\bar{y})^2}.\] Periodicidades en las correlaciones que corresponden a valores separados por 12 unidades (y múltiplos 24, 36, 48, \(\dots\)) indican estacionalidad anual. Con tendencia, la ACF decae lentamente.
# DECISIÓN: se compara la ACF con y sin tendencia para ver el efecto de quitarla
par(mfrow = c(2, 1), mar = c(3, 3, 3, 1))
acf(ise_boxcox, 48, main = "ACF ISE Box-Cox")
acf(ise_sin_tendencia, 48, main = "ACF ISE Box-Cox sin tendencia")
par(mfrow = c(1, 1))iii) Función de autocorrelación parcial (PACF)
La PACF mide la relación directa entre \(y_t\) y \(y_{t-k}\), eliminando la relación explicada por los retardos intermedios \(y_{t-1}, \dots, y_{t-k+1}\). Si se ajusta la regresión \[y_t=\beta_0+\phi_{k1}y_{t-1}+\phi_{k2}y_{t-2}+\dots+\phi_{kk}y_{t-k}+a_t,\] la autocorrelación parcial en el retardo \(k\) es el coeficiente del último retardo: \(\text{PACF}(k)=\phi_{kk}\).
par(mfrow = c(2, 1), mar = c(3, 3, 3, 1))
pacf(ise_boxcox, lag.max = 48, main = "PACF ISE Box-Cox")
pacf(ise_sin_tendencia, lag.max = 48, main = "PACF ISE Box-Cox sin tendencia")
par(mfrow = c(1, 1))iv) Información mutua promedio (AMI)
La información mutua promedio mide cuánto nos dice una variable aleatoria sobre otra y se define como \[I(X;Y)=\sum_{i}\sum_{j}p(x_i,y_j)\log_2\left(\frac{p(x_i,y_j)}{p(x_i)p(y_j)}\right).\] En el contexto de series de tiempo, la AMI cuantifica la cantidad de conocimiento obtenido sobre el valor de \(X_{t+d}\) al observar \(X_t\). Equivalentemente, es una medida de qué tanto el conocimiento de \(X\) reduce la incertidumbre acerca de \(Y\); esto implica que \(I(X;Y)=0\) si y solo si \(X\) y \(Y\) son independientes. A diferencia de la ACF, detecta también relaciones no lineales. Si se elige \(d\) alrededor del primer mínimo de la AMI, \(X_t\) y \(X_{t+d}\) son parcialmente, pero no totalmente, independientes.
# Parámetro importante: n.partitions (celdas del histograma con que se estiman las probabilidades).
# Con solo 258 datos, 50 particiones deja muchas celdas casi vacías; se puede bajar para una AMI más estable.
nonlinearTseries::mutualInformation(ise_sin_tendencia, lag.max = 48,
n.partitions = 50, units = "Bits", do.plot = TRUE)
## $time.lag
## [1] 0 1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24
## [26] 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48
##
## $mutual.information
## [1] 4.452384 1.886142 1.828053 1.603754 1.674828 1.575223 1.637594 1.576680
## [9] 1.577816 1.627761 1.826505 1.824755 2.260776 1.803471 1.763743 1.750441
## [17] 1.724847 1.674984 1.632309 1.625029 1.588814 1.711779 1.772566 1.885754
## [25] 2.016940 1.837835 1.865643 1.798491 1.829846 1.738063 1.664696 1.710037
## [33] 1.659673 1.719609 1.883828 1.836784 2.180442 1.861083 1.874075 1.798664
## [41] 1.787046 1.715643 1.711766 1.822754 1.779157 1.771955 1.917859 2.041906
## [49] 2.315325
##
## $units
## [1] "Bits"
##
## $n.partitions
## [1] 50
##
## attr(,"class")
## [1] "mutualInf"v) Autocorrelación de distancia (dACF)
La autocorrelación de distancia también permite detectar relaciones no lineales. Para cada retardo \(k\) se calcula la correlación de distancia entre \(y_t\) y \(y_{t-k}\): \[\text{dACF}(k)=\text{dCor}(y_t,y_{t-k})=\frac{\text{dCov}(y_t,y_{t-k})}{\sqrt{\text{dVar}(y_t)\,\text{dVar}(y_{t-k})}},\] donde la covarianza de distancia se construye a partir de las matrices de distancias \(|y_i-y_j|\) doblemente centradas. Toma valores entre 0 y 1, y vale 0 si y solo si las variables son independientes.
calcular_dacf_interna <- function(serie, max_lag = 20) {
sapply(1:max_lag, function(k) {
n_k <- length(serie)
x_t <- serie[(k + 1):n_k] # y_t
x_tk <- serie[1:(n_k - k)] # y_(t-k)
return(dcor(x_t, x_tk))
})
}
lags <- 1:48
dacf_ise <- calcular_dacf_interna(ise_sin_tendencia, max_lag = 48)
plot(lags, dacf_ise, type = "h", main = "dACF ISE Box-Cox sin tendencia",
xlab = "Retardo k", ylab = "Correlación de distancia", ylim = c(0, 1))
6. Ciclos y estacionalidad
Es necesario eliminar la tendencia de la serie para pasar a detectar
la estacionalidad, por eso todas las herramientas se aplican a
ise_sin_tendencia; solo en el mapa de calor y en el primer
periodograma se muestra también la serie con tendencia, para ver cómo
esta tapa el patrón estacional.
i) Mapa de calor
El mapa de calor ubica los meses en el eje vertical y los años en el horizontal, y colorea cada celda según el valor de la serie. Si el color cambia principalmente en dirección horizontal (los mismos meses siempre más oscuros o más claros) hay estacionalidad; si cambia en franjas verticales (años completos más claros u oscuros) hay ciclos o tendencia. Con la serie con tendencia dominan las franjas verticales; al quitarla debería quedar a la vista el patrón por meses.
ts_heatmap(ise_boxcox, title = "Mapa de calor - ISE Box-Cox (con tendencia)")ts_heatmap(ise_sin_tendencia, title = "Mapa de calor - ISE Box-Cox sin tendencia", color = "Reds")ii) Estadísticas descriptivas por mes
Se agrupan las observaciones por mes del año y se calcula la media y la desviación estándar de cada grupo. Si las medias son claramente distintas entre meses, hay un patrón estacional; la desviación estándar indica qué tan estable es ese mes a lo largo de los años.
ise_df_mes = data.frame(anio = floor(time(ise_sin_tendencia)),
mes = factor(month.abb[cycle(ise_sin_tendencia)], levels = month.abb),
ise = as.numeric(ise_sin_tendencia))
ise_resumen_mes = ise_df_mes %>%
group_by(mes) %>%
summarise(n = n(), media = mean(ise), sd = sd(ise))
kable(ise_resumen_mes, digits = 4, caption = "Media y desviación estándar por mes - ISE Box-Cox sin tendencia")| mes | n | media | sd |
|---|---|---|---|
| Jan | 22 | -0.0053 | 0.0013 |
| Feb | 22 | -0.0027 | 0.0015 |
| Mar | 22 | -0.0004 | 0.0013 |
| Apr | 22 | -0.0033 | 0.0030 |
| May | 22 | -0.0011 | 0.0023 |
| Jun | 22 | -0.0008 | 0.0012 |
| Jul | 21 | 0.0005 | 0.0010 |
| Aug | 21 | 0.0006 | 0.0009 |
| Sep | 21 | 0.0005 | 0.0009 |
| Oct | 21 | 0.0005 | 0.0008 |
| Nov | 21 | 0.0037 | 0.0010 |
| Dec | 21 | 0.0081 | 0.0011 |
ggplot(data = ise_resumen_mes, aes(x = mes, y = media)) +
geom_col(fill = "steelblue") +
geom_errorbar(aes(ymin = media - sd, ymax = media + sd), width = 0.3) +
geom_hline(yintercept = 0, linetype = 2) +
labs(title = "ISE Box-Cox sin tendencia - Media por mes (± 1 desviación estándar)", x = "Mes", y = "Media")
iii) Monthplot
El monthplot agrupa las observaciones según el mes del
año (enero con enero, febrero con febrero, etc.) y dibuja cada subserie
a lo largo de los años, con una línea horizontal en su media. Si las
medias son distintas para algunos meses, hay un ciclo estacional en la
serie.
monthplot(ise_sin_tendencia, main = "ISE Box-Cox sin tendencia - Comportamiento mensual",
ylab = "ISE Box-Cox sin tendencia", xlab = "Mes")
iv) Wavelet
La transformada wavelet descompone la serie en oscilaciones de distintos periodos de forma local en el tiempo: el eje horizontal es el tiempo, el vertical el periodo y el color la potencia. Una banda continua en un periodo indica un ciclo presente durante toda la serie; manchas aisladas indican ciclos o choques que aparecen solo en ciertos años. La zona sombreada (cono de influencia) es poco confiable por el efecto de los bordes.
series_wavelets <- cbind(seq(1:length(ise_sin_tendencia)), as.numeric(ise_sin_tendencia))
wt_ise <- biwavelet::wt(series_wavelets)
plot(wt_ise, type = "power.corr.norm", main = "ISE Box-Cox sin tendencia - Bias-corrected wavelet power",
ylab = "Periodo (meses)", plot.cb = TRUE)
v) Boxplot por mes
Otro enfoque para analizar patrones estacionales es trazar la distribución de las unidades de frecuencia. Esto permite examinar si cada unidad de frecuencia (cada mes) tiene una distribución única que la distingue del resto. La línea dentro de cada caja es la mediana, la caja va de \(Q_1\) a \(Q_3\) y los puntos fuera de los bigotes son posibles atípicos.
# cycle() da la posición de cada dato dentro del año (1 = enero, ..., 12 = diciembre)
ise_boxplot <- data.frame(mes = cycle(ise_sin_tendencia), ise = as.numeric(ise_sin_tendencia))
ise_boxplot$mes <- factor(month.abb[ise_boxplot$mes], levels = month.abb) # levels mantiene el orden cronológico
ggplot(data = ise_boxplot, aes(x = mes, y = ise)) +
geom_boxplot() +
labs(title = "ISE Box-Cox sin tendencia - Distribución por mes", x = "Mes", y = "ISE Box-Cox sin tendencia")
vi) Periodograma
Dada una serie de longitud \(T\),
las frecuencias básicas o de Fourier son \[f_j=\frac{j}{T},\quad j=1,2,\dots,T/2,\]
así que \(1/T\leq f_j\leq 1/2\); el
valor máximo \(f_j=0.5\) se conoce como
la frecuencia de Nyquist. Una serie periódica \(Z_t\) se puede representar como \[Z_t=\mu+\sum_{j=1}^{T/2}A_j\sin(\omega_j t)+
\sum_{j=1}^{T/2}B_j\cos(\omega_j t)+a_t,\qquad \omega_j=2\pi
f_j,\] con estimadores \[\hat{A}_j=\frac{2}{T}\sum_{t=1}^{T}Z_t\sin(\omega_j
t),\qquad \hat{B}_j=\frac{2}{T}\sum_{t=1}^{T}Z_t\cos(\omega_j t),\qquad
\hat{R}_j^2=\hat{A}_j^2+\hat{B}_j^2.\] La contribución de cada
onda a la varianza es \(\hat{R}_j^2/2\), por lo que las ondas con
amplitudes grandes son las importantes en la explicación de la
variabilidad de la serie. Se le da el nombre de periodograma a la
representación de \[I(f_j)=\frac{T\hat{R}_j^2}{2},\qquad 1/T\leq
f_j\leq 1/2,\] en función de la frecuencia. En una serie mensual
estacional de periodo \(s=12\) se
espera un valor alto del periodograma en \(f=1/12\), pero también en \(f=j/12\), es decir, \(1/6, 1/4, 1/3\), que son armónicos del
periodo estacional. Como la serie es un objeto ts de
frecuencia 12, spectrum reporta la frecuencia en ciclos por
año: \(f = 1\) es el ciclo de 12 meses,
\(f = 2\) el de 6 meses, \(f = 3\) el de 4 meses, etc. Primero se
compara el periodograma con y sin tendencia: con tendencia, casi toda la
potencia se concentra cerca de la frecuencia 0 y tapa los picos
estacionales.
par(mfrow = c(2, 1), mar = c(4, 4, 3, 1))
spectrum(ise_boxcox, log = "no", main = "Periodograma - ISE Box-Cox con tendencia") # la tendencia concentra la potencia cerca de 0
spectrum(ise_sin_tendencia, log = "no", main = "Periodograma - ISE Box-Cox sin tendencia")
par(mfrow = c(1, 1))periodograma_ise = spectrum(ise_sin_tendencia, log = "no", main = "Periodograma - ISE Box-Cox sin tendencia")
abline(v = 1:6, lty = 2, col = 2) # frecuencia anual (f = 1) y sus armónicos
# Ubicación del máximo
ubicacion_ise = which.max(periodograma_ise$spec)
sprintf("Frecuencia donde se maximiza el periodograma: %.4f ciclos por año", periodograma_ise$freq[ubicacion_ise])## [1] "Frecuencia donde se maximiza el periodograma: 4.0000 ciclos por año"sprintf("Periodo correspondiente: %.2f meses", 12 / periodograma_ise$freq[ubicacion_ise])## [1] "Periodo correspondiente: 3.00 meses"# Las 8 frecuencias con mayor potencia y su periodo en meses
orden_frec = order(periodograma_ise$spec, decreasing = TRUE)[1:8]
tabla_frecuencias = data.frame(frecuencia_ciclos_anio = periodograma_ise$freq[orden_frec],
periodo_meses = 12 / periodograma_ise$freq[orden_frec],
potencia = periodograma_ise$spec[orden_frec])
kable(tabla_frecuencias, digits = 4, caption = "Frecuencias con mayor potencia - ISE Box-Cox sin tendencia")| frecuencia_ciclos_anio | periodo_meses | potencia |
|---|---|---|
| 4.0000 | 3.0000 | 0 |
| 0.9778 | 12.2727 | 0 |
| 2.0000 | 6.0000 | 0 |
| 1.0222 | 11.7391 | 0 |
| 3.0222 | 3.9706 | 0 |
| 2.9778 | 4.0299 | 0 |
| 4.9778 | 2.4107 | 0 |
| 5.0222 | 2.3894 | 0 |
vii) Periodograma suavizado
Como el periodograma sin suavizar puede tener muchos picos que en verdad no son grandes, se suaviza con pesos simétricos \(p_i\) en una ventana \(q\): \[\hat{I}(f_j)=\sum_{i=-q}^{q}p_iI(f_{j+i}).\]
# Parámetro importante: span (anchos de los suavizadores de Daniell modificados). Más grande = más suave,
# pero picos cercanos se pueden fundir en uno solo.
par(mfrow = c(3, 1), mar = c(4, 4, 3, 1))
spectrum(ise_sin_tendencia, log = "no", span = 5, main = "Periodograma - ISE Box-Cox sin tendencia - span = 5")
spectrum(ise_sin_tendencia, log = "no", span = c(5, 5), main = "Periodograma - ISE Box-Cox sin tendencia - span = c(5,5)")
spectrum(ise_sin_tendencia, log = "no", span = c(2, 2), main = "Periodograma - ISE Box-Cox sin tendencia - span = c(2,2)")
par(mfrow = c(1, 1))7. Modelamiento de la estacionalidad
Una vez estimada y quitada la tendencia, el modelo para la serie
Box-Cox sin tendencia es \[y_t-\hat{\mu}_t =
s_t + a_t,\] donde \(s_t\) es la
componente estacional de periodo \(s =
12\) (\(s_t = s_{t+12}\)). La
estacionalidad se modela con dos alternativas determinísticas, ambas
ajustadas por mínimos cuadrados sobre ise_sin_tendencia. En
las dos, lo que ajusta el modelo es la componente estacional estimada
\(\hat{s}_t\), y lo que sobra, \(\hat{a}_t = y_t-\hat{\mu}_t-\hat{s}_t\), es
la serie sin tendencia ni estacionalidad.
i) Variables dummy estacionales
Se crea una variable indicadora por cada mes y se ajusta una regresión lineal de la serie sobre ellas: \[y_t-\hat{\mu}_t = \delta_1 \gamma_{1,t} + \delta_2 \gamma_{2,t} + \dots + \delta_{12}\gamma_{12,t} + a_t,\] donde \(\gamma_{i,t}=1\) si la observación \(t\) corresponde al mes \(i\) y \(0\) en otro caso. Cada coeficiente \(\delta_i\) es el efecto promedio de ese mes, así que el patrón estacional puede tener cualquier forma. Usa 12 parámetros; en la práctica se ajustan 11 dummies más el intercepto (el mes base es enero) para evitar colinealidad perfecta.
# Como factor, lm crea automáticamente las 11 dummies (enero es el mes base, recogido por el intercepto)
mes_ise = factor(month.abb[cycle(ise_sin_tendencia)], levels = month.abb)
lm_estacional_dummy = lm(ise_sin_tendencia ~ mes_ise)
summary(lm_estacional_dummy)##
## Call:
## lm(formula = ise_sin_tendencia ~ mes_ise)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0127174 -0.0006434 0.0000375 0.0006984 0.0058051
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) -0.0052701 0.0003205 -16.446 < 2e-16 ***
## mes_iseFeb 0.0025575 0.0004532 5.643 4.58e-08 ***
## mes_iseMar 0.0048714 0.0004532 10.749 < 2e-16 ***
## mes_iseApr 0.0019513 0.0004532 4.306 2.40e-05 ***
## mes_iseMay 0.0041925 0.0004532 9.251 < 2e-16 ***
## mes_iseJun 0.0045134 0.0004532 9.959 < 2e-16 ***
## mes_iseJul 0.0058161 0.0004586 12.684 < 2e-16 ***
## mes_iseAug 0.0058929 0.0004586 12.851 < 2e-16 ***
## mes_iseSep 0.0057841 0.0004586 12.614 < 2e-16 ***
## mes_iseOct 0.0057971 0.0004586 12.642 < 2e-16 ***
## mes_iseNov 0.0089347 0.0004586 19.485 < 2e-16 ***
## mes_iseDec 0.0134022 0.0004586 29.227 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.001503 on 246 degrees of freedom
## Multiple R-squared: 0.8321, Adjusted R-squared: 0.8245
## F-statistic: 110.8 on 11 and 246 DF, p-value: < 2.2e-16fit_estacional_dummy = ts(fitted(lm_estacional_dummy), start = start(ise_sin_tendencia), frequency = 12)
ise_sin_estacional_dummy = ise_sin_tendencia - fit_estacional_dummy
par(mfrow = c(2, 1), mar = c(3, 3, 3, 1))
plot(ise_sin_tendencia, col = "grey50", main = "ISE Box-Cox sin tendencia - Estacionalidad con variables dummy")
lines(fit_estacional_dummy, lwd = 2, col = 2)
plot(ise_sin_estacional_dummy, main = "ISE Box-Cox sin tendencia ni estacionalidad - dummy")
abline(h = 0, lty = 2)
par(mfrow = c(1, 1))
# Patrón estacional estimado: efecto de cada mes en un año
patron_dummy = tapply(fit_estacional_dummy, cycle(fit_estacional_dummy), mean)
plot(1:12, patron_dummy, type = "b", pch = 19, col = 2, xaxt = "n",
main = "Patrón estacional estimado - dummy", xlab = "Mes", ylab = "Efecto estacional")
axis(1, at = 1:12, labels = month.abb)
abline(h = 0, lty = 2)
ii) Regresión armónica (series de Fourier)
La estacionalidad se representa como una suma de ondas seno y coseno de la frecuencia anual y sus armónicos: \[s_t=\sum_{k=1}^{K}\left[\beta_{1k}\cos\left(\frac{2\pi k t}{12}\right)+\beta_{2k}\sin\left(\frac{2\pi k t}{12}\right)\right].\] Gracias a la identidad \(A\cos(\omega t+\varphi)=\beta_1\cos(\omega t)+\beta_2\sin(\omega t)\), con amplitud \(A=\sqrt{\beta_1^2+\beta_2^2}\), el modelo es lineal en los parámetros y se estima por mínimos cuadrados. El parámetro a variar es \(K\), el número de armónicos (\(1 \leq K \leq s/2 = 6\)): con \(K\) pequeño el patrón es suave y usa pocos parámetros (\(2K\)); con \(K = 6\) reproduce exactamente el modelo de dummies. Los armónicos \(k = 1, 2, 3, \dots\) corresponden a las frecuencias \(f = 1, 2, 3, \dots\) ciclos por año que aparecen como picos en el periodograma.
Primero se ajusta un solo armónico, la onda de periodo 12 meses:
# t_ise = 1, ..., T ya se definió en la sección de tendencia determinística
z1_ise = cos(2 * pi * t_ise / 12)
z2_ise = sin(2 * pi * t_ise / 12)
lm_estacional_armonico = lm(ise_sin_tendencia ~ z1_ise + z2_ise)
summary(lm_estacional_armonico)##
## Call:
## lm(formula = ise_sin_tendencia ~ z1_ise + z2_ise)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0135894 -0.0016603 0.0000183 0.0015929 0.0092758
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 2.128e-05 1.958e-04 0.109 0.914
## z1_ise 1.317e-03 2.769e-04 4.757 3.29e-06 ***
## z2_ise -2.089e-03 2.770e-04 -7.543 8.10e-13 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.003145 on 255 degrees of freedom
## Multiple R-squared: 0.2377, Adjusted R-squared: 0.2317
## F-statistic: 39.76 on 2 and 255 DF, p-value: 9.34e-16sprintf("Amplitud estimada de la onda anual: %.4f",
sqrt(coef(lm_estacional_armonico)[2]^2 + coef(lm_estacional_armonico)[3]^2))## [1] "Amplitud estimada de la onda anual: 0.0025"fit_estacional_armonico = ts(fitted(lm_estacional_armonico), start = start(ise_sin_tendencia), frequency = 12)
plot(ise_sin_tendencia, col = "grey50", main = "ISE Box-Cox sin tendencia - Un armónico (periodo 12 meses)")
lines(fit_estacional_armonico, lwd = 2, col = 4)
Luego se compara el ajuste según el número de armónicos \(K\), para ver cuánto mejora con cada
armónico adicional. La función fourier del paquete
forecast construye las columnas seno y coseno de los \(K\) armónicos:
# Comparación del número de armónicos (K = 6 equivale a las dummies)
tabla_K = data.frame(K = 1:6, parametros = NA, R2 = NA, R2_ajustado = NA, AIC = NA, BIC = NA)
for (K in 1:6) {
X_fourier = fourier(ise_sin_tendencia, K = K)
lm_K = lm(ise_sin_tendencia ~ X_fourier)
tabla_K$parametros[K] = length(coef(lm_K))
tabla_K$R2[K] = summary(lm_K)$r.squared
tabla_K$R2_ajustado[K] = summary(lm_K)$adj.r.squared
tabla_K$AIC[K] = AIC(lm_K)
tabla_K$BIC[K] = BIC(lm_K)
}
kable(tabla_K, digits = 4, caption = "Regresión armónica según el número de armónicos K")| K | parametros | R2 | R2_ajustado | AIC | BIC |
|---|---|---|---|---|---|
| 1 | 3 | 0.2377 | 0.2317 | -2235.977 | -2221.765 |
| 2 | 5 | 0.3504 | 0.3402 | -2273.265 | -2251.947 |
| 3 | 7 | 0.5275 | 0.5162 | -2351.387 | -2322.964 |
| 4 | 9 | 0.7116 | 0.7023 | -2474.712 | -2439.183 |
| 5 | 11 | 0.8210 | 0.8138 | -2593.830 | -2551.194 |
| 6 | 12 | 0.8321 | 0.8245 | -2608.247 | -2562.059 |
Para ver cómo cambia el patrón estacional al agregar armónicos, se grafican los ajustes con \(K = 3, 4\) y \(5\) antes de pasar a \(K = 6\). Con pocos armónicos la curva es una onda suave que no alcanza a seguir los cambios bruscos entre meses; cada armónico adicional \(k\) suma una onda de periodo más corto (\(12/k\) meses: 4, 3 y 2.4 meses para \(k = 3, 4, 5\)), que permite capturar picos y caídas más abruptas de un mes a otro.
# Ajustes con K = 3, 4 y 5 armónicos para ver la transición antes de K = 6
X_fourier3_ise = fourier(ise_sin_tendencia, K = 3)
X_fourier4_ise = fourier(ise_sin_tendencia, K = 4)
X_fourier5_ise = fourier(ise_sin_tendencia, K = 5)
lm_estacional_fourier3 = lm(ise_sin_tendencia ~ X_fourier3_ise)
lm_estacional_fourier4 = lm(ise_sin_tendencia ~ X_fourier4_ise)
lm_estacional_fourier5 = lm(ise_sin_tendencia ~ X_fourier5_ise)
fit_estacional_fourier3 = ts(fitted(lm_estacional_fourier3), start = start(ise_sin_tendencia), frequency = 12)
fit_estacional_fourier4 = ts(fitted(lm_estacional_fourier4), start = start(ise_sin_tendencia), frequency = 12)
fit_estacional_fourier5 = ts(fitted(lm_estacional_fourier5), start = start(ise_sin_tendencia), frequency = 12)
ise_sin_estacional_fourier3 = ise_sin_tendencia - fit_estacional_fourier3
ise_sin_estacional_fourier4 = ise_sin_tendencia - fit_estacional_fourier4
ise_sin_estacional_fourier5 = ise_sin_tendencia - fit_estacional_fourier5
par(mfrow = c(3, 2), mar = c(3, 3, 3, 1))
plot(ise_sin_tendencia, col = "grey50", main = "ISE Box-Cox sin tendencia - Fourier K = 3")
lines(fit_estacional_fourier3, lwd = 2, col = 4)
plot(ise_sin_estacional_fourier3, main = "ISE sin tendencia ni estacionalidad - Fourier K = 3")
abline(h = 0, lty = 2)
plot(ise_sin_tendencia, col = "grey50", main = "ISE Box-Cox sin tendencia - Fourier K = 4")
lines(fit_estacional_fourier4, lwd = 2, col = 4)
plot(ise_sin_estacional_fourier4, main = "ISE sin tendencia ni estacionalidad - Fourier K = 4")
abline(h = 0, lty = 2)
plot(ise_sin_tendencia, col = "grey50", main = "ISE Box-Cox sin tendencia - Fourier K = 5")
lines(fit_estacional_fourier5, lwd = 2, col = 4)
plot(ise_sin_estacional_fourier5, main = "ISE sin tendencia ni estacionalidad - Fourier K = 5")
abline(h = 0, lty = 2)
par(mfrow = c(1, 1))# Patrón estacional de un año con K = 3, 4 y 5, frente al de dummies (al que llega K = 6)
patron_fourier3 = tapply(fit_estacional_fourier3, cycle(fit_estacional_fourier3), mean)
patron_fourier4 = tapply(fit_estacional_fourier4, cycle(fit_estacional_fourier4), mean)
patron_fourier5 = tapply(fit_estacional_fourier5, cycle(fit_estacional_fourier5), mean)
plot(1:12, patron_dummy, type = "b", pch = 19, lwd = 2, col = 1, xaxt = "n",
ylim = range(c(patron_dummy, patron_fourier3, patron_fourier4, patron_fourier5)),
main = "Patrón estacional: transición de Fourier K = 3, 4, 5 hacia dummies", xlab = "Mes", ylab = "Efecto estacional")
lines(1:12, patron_fourier3, type = "b", pch = 15, col = 2)
lines(1:12, patron_fourier4, type = "b", pch = 17, col = 3)
lines(1:12, patron_fourier5, type = "b", pch = 18, col = 4)
axis(1, at = 1:12, labels = month.abb)
abline(h = 0, lty = 2)
legend("topleft", c("Dummy", "Fourier K = 3", "Fourier K = 4", "Fourier K = 5"),
col = c(1, 2, 3, 4), pch = c(19, 15, 17, 18), lwd = c(2, 1, 1, 1), cex = 0.8)
Finalmente se usa \(K = 6\), el máximo posible para \(s = 12\).
# Parámetro importante: K (número de armónicos, de 1 a 6).
# DECISIÓN: K = 6, el máximo para s = 12 (equivale al modelo de dummies)
K_fourier = 6
X_fourier_ise = fourier(ise_sin_tendencia, K = K_fourier)
lm_estacional_fourier = lm(ise_sin_tendencia ~ X_fourier_ise)
summary(lm_estacional_fourier)##
## Call:
## lm(formula = ise_sin_tendencia ~ X_fourier_ise)
##
## Residuals:
## Min 1Q Median 3Q Max
## -0.0127174 -0.0006434 0.0000375 0.0006984 0.0058051
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 3.933e-05 9.360e-05 0.420 0.675
## X_fourier_iseS1-12 -2.068e-03 1.324e-04 -15.626 < 2e-16 ***
## X_fourier_iseC1-12 1.369e-03 1.324e-04 10.342 < 2e-16 ***
## X_fourier_iseS2-12 -9.540e-04 1.324e-04 -7.207 7.01e-12 ***
## X_fourier_iseC2-12 1.439e-03 1.324e-04 10.869 < 2e-16 ***
## X_fourier_iseS3-12 -1.608e-03 1.324e-04 -12.145 < 2e-16 ***
## X_fourier_iseC3-12 1.396e-03 1.324e-04 10.549 < 2e-16 ***
## X_fourier_iseS4-12 -1.157e-03 1.324e-04 -8.737 3.76e-16 ***
## X_fourier_iseC4-12 1.833e-03 1.324e-04 13.850 < 2e-16 ***
## X_fourier_iseS5-12 4.549e-06 1.324e-04 0.034 0.973
## X_fourier_iseC5-12 1.679e-03 1.324e-04 12.684 < 2e-16 ***
## X_fourier_iseC6-12 3.763e-04 9.360e-05 4.020 7.74e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.001503 on 246 degrees of freedom
## Multiple R-squared: 0.8321, Adjusted R-squared: 0.8245
## F-statistic: 110.8 on 11 and 246 DF, p-value: < 2.2e-16fit_estacional_fourier = ts(fitted(lm_estacional_fourier), start = start(ise_sin_tendencia), frequency = 12)
ise_sin_estacional_fourier = ise_sin_tendencia - fit_estacional_fourier
par(mfrow = c(2, 1), mar = c(3, 3, 3, 1))
plot(ise_sin_tendencia, col = "grey50",
main = sprintf("ISE Box-Cox sin tendencia - Estacionalidad con Fourier (K = %d)", K_fourier))
lines(fit_estacional_fourier, lwd = 2, col = 4)
plot(ise_sin_estacional_fourier, main = sprintf("ISE Box-Cox sin tendencia ni estacionalidad - Fourier (K = %d)", K_fourier))
abline(h = 0, lty = 2)
par(mfrow = c(1, 1))
# Verificación: con K = 6 los valores ajustados coinciden con los del modelo de dummies (debe dar TRUE)
ncol(X_fourier_ise) # 11 columnas: el seno del armónico 6 es idénticamente 0 y se omite## [1] 11all.equal(as.numeric(fit_estacional_fourier), as.numeric(fit_estacional_dummy))## [1] TRUEiii) Comparación de las dos alternativas
Con \(K = 6\) ambos modelos dan exactamente los mismos valores ajustados, así que en las gráficas las curvas de dummy y Fourier se superponen y los residuos son idénticos; la comparación sirve para verificar esa equivalencia y para revisar sobre los residuos si la estacionalidad quedó bien modelada. Se comparan el patrón estacional estimado, el ajuste en los últimos años, los residuos con su ACF y el periodograma de los residuos. Si la estacionalidad quedó bien modelada, la ACF de los residuos ya no debería mostrar picos en los rezagos 12, 24, 36 y el periodograma no debería tener picos en las frecuencias \(f = 1, 2, 3, \dots\) ciclos por año.
tabla_estacional = data.frame(modelo = c("Dummy", sprintf("Fourier (K = %d)", K_fourier)),
parametros = c(length(coef(lm_estacional_dummy)), length(coef(lm_estacional_fourier))),
R2_ajustado = c(summary(lm_estacional_dummy)$adj.r.squared,
summary(lm_estacional_fourier)$adj.r.squared),
AIC = c(AIC(lm_estacional_dummy), AIC(lm_estacional_fourier)),
BIC = c(BIC(lm_estacional_dummy), BIC(lm_estacional_fourier)))
kable(tabla_estacional, digits = 4, caption = "Comparación de los modelos de estacionalidad")| modelo | parametros | R2_ajustado | AIC | BIC |
|---|---|---|---|---|
| Dummy | 12 | 0.8245 | -2608.247 | -2562.059 |
| Fourier (K = 6) | 12 | 0.8245 | -2608.247 | -2562.059 |
# Patrón estacional de un año con cada alternativa
patron_fourier = tapply(fit_estacional_fourier, cycle(fit_estacional_fourier), mean)
plot(1:12, patron_dummy, type = "b", pch = 19, col = 2, xaxt = "n",
ylim = range(c(patron_dummy, patron_fourier)),
main = "Patrón estacional estimado: dummy vs Fourier", xlab = "Mes", ylab = "Efecto estacional")
lines(1:12, patron_fourier, type = "b", pch = 17, col = 4)
axis(1, at = 1:12, labels = month.abb)
abline(h = 0, lty = 2)
legend("topleft", c("Dummy", sprintf("Fourier (K = %d)", K_fourier)), col = c(2, 4),
pch = c(19, 17), lwd = 1, cex = 0.8)
# Ajuste en los últimos años
plot(window(ise_sin_tendencia, start = c(2022, 1)), col = "grey50", lwd = 2,
main = "ISE Box-Cox sin tendencia - Ajuste estacional 2022-2026", ylab = "ISE Box-Cox sin tendencia")
lines(window(fit_estacional_dummy, start = c(2022, 1)), col = 2, lwd = 2)
lines(window(fit_estacional_fourier, start = c(2022, 1)), col = 4, lwd = 2, lty = 2)
abline(h = 0, lty = 2)
legend("topleft", c("Serie", "Dummy", "Fourier"), col = c("grey50", 2, 4),
lwd = 2, lty = c(1, 1, 2), cex = 0.8)
par(mfrow = c(2, 2), mar = c(3, 3, 3, 1))
plot(ise_sin_estacional_dummy, main = "Residuo - dummy")
abline(h = 0, lty = 2)
plot(ise_sin_estacional_fourier, main = sprintf("Residuo - Fourier (K = %d)", K_fourier))
abline(h = 0, lty = 2)
acf(ise_sin_estacional_dummy, lag.max = 48, main = "ACF residuo - dummy")
acf(ise_sin_estacional_fourier, lag.max = 48, main = "ACF residuo - Fourier")
par(mfrow = c(1, 1))# Si la estacionalidad quedó bien modelada, los picos en f = 1, 2, 3, ... desaparecen
par(mfrow = c(2, 1), mar = c(4, 4, 3, 1))
spectrum(ise_sin_estacional_dummy, log = "no", main = "Periodograma residuo - dummy")
abline(v = 1:6, lty = 2, col = 2)
spectrum(ise_sin_estacional_fourier, log = "no", main = sprintf("Periodograma residuo - Fourier (K = %d)", K_fourier))
abline(v = 1:6, lty = 2, col = 2)
par(mfrow = c(1, 1))