Sea \(\{y_t\}_{t=1}^n\) una serie temporal que representa la concentración promedio de dióxido de carbono (CO₂) medida en partes por millón (ppm). El objetivo del modelo es describir la dinámica temporal de la media condicional de la serie, permitiendo que esta evolucione de forma flexible en el tiempo.
Se asume que, condicional a la información pasada \(\mathcal{F}_{t-1}\), la variable \(y_t\) sigue una distribución Normal:
\[ y_t \mid \mathcal{F}_{t-1} \sim \mathcal{N}(u_t, \sigma_\varepsilon^2), \]
donde: - \(u_t\) es la media condicional dinámica, - \(\sigma_\varepsilon^2\) es la varianza del término de error, asumida constante.
Los modelos GAS (Generalized Autoregressive Score) se caracterizan porque los parámetros dinámicos se actualizan utilizando el score, es decir, el gradiente de la log-verosimilitud condicional con respecto al parámetro dinámico.
En este caso, el parámetro dinámico es la media \(u_t\). Para una distribución Normal, el score con respecto a la media es:
\[ s_t = \frac{\partial \log f(y_t \mid u_t)}{\partial u_t} = \frac{y_t - u_t}{\sigma_\varepsilon^2}. \]
Este término mide el error de predicción estandarizado y es el mecanismo mediante el cual la información nueva \(y_t\) actualiza la media condicional.
La dinámica del modelo GAS para la media condicional se define como:
\[ u_{t+1} = w + \alpha s_t + \beta u_t, \]
o, de forma explícita,
\[ u_{t+1} = w + \alpha \frac{(y_t - u_t)}{\sigma_\varepsilon^2} + \beta u_t, \]
donde: - \(w\) es un término constante, - \(\alpha\) controla la intensidad con la que el score afecta la actualización, - \(\beta\) captura la persistencia temporal de la media condicional.
datos_co2 = read.csv2("C:/Users/maria/OneDrive/Escritorio/SERIES/co.csv", header = TRUE, fileEncoding = "UTF-8")
datos_co2$promedio <- as.numeric(datos_co2$promedio)
n <- length(datos_co2$promedio)
y <- datos_co2$promedio
##################################
#### Modelo GAS para CO2
##################################
GAS_MODEL_CO2 <- function(w, alpha, beta, sigma_epsilon){
# Inicializar vectores
u <- rep(0, n)
# Condición inicial
u[1] <- y[1]
# Modelo GAS
for(i in 1:(n-1)){
u[i+1] <- w + alpha*(y[i] - u[i])/sigma_epsilon^2 + beta*u[i]
}
# Graficar resultados
plot(y, type="l", main="Modelo GAS para CO2",
col="blue", xlab="Meses desde Marzo 1958",
ylab="ppm", lwd=1.5)
lines(u, col="red", lwd=2)
legend("topleft",
legend=c("Datos CO2", "Media estimada (u)"),
col=c("blue", "red"), lty=1, lwd=2, cex=0.8)
return(list(y = y, u = u))
}
res <- GAS_MODEL_CO2(
w = 0.1,
alpha = 0.5,
beta = 0.9,
sigma_epsilon = sd(y, na.rm = TRUE)
)Los parámetros del modelo GAS se estiman mediante Máxima Verosimilitud (MV). Bajo el supuesto de normalidad condicional,
\[ y_t \mid \mathcal{F}_{t-1} \sim \mathcal{N}(u_t, \sigma_\varepsilon^2), \]
la log-verosimilitud condicional de una observación \(y_t\) es
\[ \ell_t = -\frac{1}{2}\log(2\pi) - \log(\sigma_\varepsilon) - \frac{(y_t - u_t)^2}{2\sigma_\varepsilon^2}. \]
La log-verosimilitud total del modelo se obtiene sumando estas contribuciones para \(t = 2, \dots, n\), dado que el primer valor se utiliza como condición inicial:
\[ \mathcal{L}(w,\alpha,\beta,\sigma_\varepsilon) = \sum_{t=2}^{n} \ell_t. \]
El vector de parámetros a estimar es
\[ \theta = (w,\alpha,\beta,\sigma_\varepsilon). \]
Para garantizar la positividad de la desviación estándar, se reparametriza:
\[ \sigma_\varepsilon = \exp(\eta), \]
donde \(\eta \in \mathbb{R}\) es el parámetro que se optimiza numéricamente.
La maximización de la log-verosimilitud se realiza mediante el algoritmo BFGS, un método cuasi–Newton que utiliza aproximaciones de segundo orden. El Hessiano estimado en el óptimo permite aproximar la matriz de varianzas y covarianzas de los estimadores:
\[ \widehat{\mathrm{Var}}(\hat{\theta}) = \left(-\nabla^2 \mathcal{L}(\hat{\theta})\right)^{-1}. \]
Los errores estándar se obtienen como la raíz cuadrada de los elementos diagonales de esta matriz. En el caso de \(\sigma_\varepsilon\), el error estándar se ajusta usando la regla delta debido a la transformación exponencial.
Para cada parámetro se calcula el estadístico t:
\[ t_j = \frac{\hat{\theta}_j}{\mathrm{se}(\hat{\theta}_j)}, \]
y se evalúa su significancia estadística utilizando el criterio asintótico:
\[ |t_j| > 1.96, \]
correspondiente a un nivel de significancia del 5 %.
El procedimiento produce estimaciones para: - \(w\): nivel base de la media condicional, - \(\alpha\): intensidad con que los errores de predicción afectan la actualización de la media, - \(\beta\): persistencia temporal de la media condicional, - \(\sigma_\varepsilon\): desviación estándar del error.
La log-verosimilitud máxima alcanzada es:
\[ \mathcal{L}_{\max} = -1083.064, \]
lo que resume el ajuste global del modelo a la serie de CO₂ bajo los supuestos planteados.
Los residuos se definen como:
\[ \hat{\varepsilon}_t = y_t - \hat{u}_t. \]
Las estadísticas descriptivas indican que: - la media de los residuos es aproximadamente cero, - la desviación estándar es cercana a \(0.92\).
Esto sugiere que el modelo captura adecuadamente el nivel medio de la serie, sin sesgo sistemático.
El análisis gráfico de los residuos incluye: 1. Serie temporal de residuos, que permite evaluar autocorrelación o patrones sistemáticos no capturados por el modelo. 2. Histograma con densidad normal superpuesta, que permite evaluar visualmente el supuesto de normalidad.
En conjunto, estos diagnósticos permiten verificar si el modelo GAS para la media condicional proporciona una descripción razonable de la dinámica de corto plazo de la serie de CO₂.
##################################
#### Estimación por Máxima Verosimilitud
##################################
GAS_NORMAL_CO2 <- function(q){
# Parámetros a estimar
w <- q[1]
alpha <- q[2]
beta <- q[3]
sigma_epsilon <- exp(q[4]) # Usar exponencial para asegurar positividad
n <- length(y)
u <- rep(0, n)
# Condición inicial
u[1] <- y[1]
# Calcular media condicional
for (i in 1:(n-1)){
u[i+1] <- w + alpha*(y[i] - u[i])/sigma_epsilon^2 + beta*u[i]
}
# Log-verosimilitud (excluyendo el primer punto)
sum_ll <- 0
for (i in 2:n){
sum_ll <- sum_ll + dnorm(y[i], mean = u[i], sd = sigma_epsilon, log = TRUE)
}
return(sum_ll)
}
##################################
#### Estimación de parámetros
##################################
# Valor inicial para sigma_epsilon (estimación basada en varianza de la serie)
sigma_init <- log(sd(y, na.rm = TRUE))
# Valores iniciales para la optimización
initial_params <- c(
w = mean(y), # Valor inicial para w (intercepto)
alpha = 0.1, # Sensibilidad a los errores
beta = 0.9, # Persistencia
log_sigma = sigma_init # Log de la desviación estándar
)
# Optimización por máxima verosimilitud
result <- optim(
par = initial_params,
fn = GAS_NORMAL_CO2,
method = "BFGS",
control = list(fnscale = -1, maxit = 1000), # Maximizar
hessian = TRUE
)
# Extraer parámetros estimados
param_est <- result$par
param_est[4] <- exp(param_est[4]) # Transformar sigma de vuelta
# Calcular errores estándar
if (!is.null(result$hessian)) {
hessian_inv <- solve(-result$hessian)
stderrors <- sqrt(diag(hessian_inv))
stderrors[4] <- stderrors[4] * param_est[4] # Ajustar para sigma
} else {
stderrors <- rep(NA, length(param_est))
}
# Crear tabla de resultados
estim <- data.frame(
Parámetro = c("w", "alpha", "beta", "sigma_epsilon"),
Estimación = round(param_est, 4),
Error_Std = round(stderrors, 4),
t_stat = round(param_est/stderrors, 4),
Significativo = abs(param_est/stderrors) > 1.96
)
print("Resultados de la estimación:")## [1] "Resultados de la estimación:"
## Parámetro Estimación Error_Std t_stat Significativo
## w w 0.1841 0.5930 0.3105 FALSE
## alpha alpha 1.4059 0.0720 19.5165 TRUE
## beta beta 0.9999 0.0016 609.9360 TRUE
## log_sigma sigma_epsilon 0.9230 0.0230 40.2151 TRUE
##
## Log-verosimilitud máxima: -1083.064
# Ejecutar modelo con parámetros estimados
cat("\nEjecutando modelo GAS con parámetros estimados...\n")##
## Ejecutando modelo GAS con parámetros estimados...
sigma_epsilon_est <- param_est[4]
modelo_final <- GAS_MODEL_CO2(
w = param_est[1],
alpha = param_est[2],
beta = param_est[3],
sigma_epsilon = sigma_epsilon_est
)# Calcular y mostrar residuos
residuos <- y - modelo_final$u
cat("\nEstadísticas de los residuos:\n")##
## Estadísticas de los residuos:
## Media: -0.0026
## Desviación estándar: 0.923
## null device
## 1
par(mfrow = c(1,2)) # Ventana 1x2
# Gráfico de residuos
plot(residuos, type="l",
main="Residuos del modelo GAS",
xlab="Meses", ylab="Residuos",
col="darkgreen")
abline(h=0, col="red", lty=2)
# Histograma de residuos
hist(residuos, breaks=30, freq=FALSE,
main="Distribución de residuos",
xlab="Residuos",
col="lightgreen", border="darkgreen")
curve(dnorm(x, mean=mean(residuos), sd=sd(residuos)),
add=TRUE, col="red", lwd=2)
par(mfrow = c(1,1)) # Volver a normalAunque el modelo ajusta adecuadamente la media local y produce residuos centrados, la serie de CO₂ presenta una tendencia estructural de largo plazo. Por ello, este enfoque GAS describe principalmente la dinámica de corto plazo, y podría ampliarse incorporando componentes de tendencia o modelos más generales de espacio de estados para capturar la evolución secular de la concentración de CO₂.