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.
# 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")
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")
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")
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")
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")
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")
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")
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")
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")
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")
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 se usan
solo los datos \(\{x_{t-k/2}, \dots, x_t,
\dots, x_{t+k/2}\}\) para predecir \(x_t\) mediante regresión, y luego se
establece \(m_t = \hat{x}_t\). Primero,
una cierta proporción de vecinos más cercanos a \(x_t\) se incluye en un esquema de
ponderación; los valores más cercanos a \(x_t\) en el tiempo obtienen más peso. Luego
se utiliza una regresión ponderada robusta 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 lowess.
# 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")
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")
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")
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")
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")
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")
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")
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")
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")
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")
par(mfrow = c(1, 1))vii) 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")
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.
i) 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")
ii) 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)
iii) 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")
iv) Periodograma suavizado
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 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}).\]
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, etc.
# 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. Eliminación de la componente estacional
Al igual que con la tendencia, hay al menos dos formas de eliminar la
componente estacional: estimándola y luego quitándola, o usando
diferencias estacionales \((1-B^s)^D\),
donde \(s\) es el periodo del ciclo
estacional. Para \(s = 12\) y \(D = 1\): \[(1-B^{12})Y_t = Y_t - Y_{t-12},\] es
decir, a cada observación se le resta la del mismo mes del año anterior.
En R se usa diff(series, lag = s, differences = D).
i) Diferencia estacional
# DECISIÓN: s = 12 (periodo anual identificado) y D = 1
ise_diff_estacional = diff(ise_sin_tendencia, lag = 12, differences = 1)
par(mfrow = c(2, 1), mar = c(4, 4, 3, 1))
plot(ise_diff_estacional, main = "ISE Box-Cox - Diferencia estacional", ylab = "Diferencia estacional", xlab = "Tiempo")
spectrum(ise_diff_estacional, log = "no", main = "Periodograma - ISE Box-Cox con diferencia estacional")
par(mfrow = c(1, 1))