#Función de Simulación

El modelo GAS–Poisson considerado se define de la siguiente forma:

\[ y_t \mid \lambda_t \sim \text{Poisson}(\lambda_t), \qquad \lambda_t = \exp(\alpha_t), \]

donde el estado dinámico \(\alpha_t\) evoluciona según la ecuación GAS de primer orden:

\[ \alpha_t = \omega + A s_{t-1} + B \alpha_{t-1}, \]

y el score de la log-verosimilitud condicional Poisson está dado por:

\[ s_{t-1} = \frac{\partial \log f(y_{t-1} \mid \alpha_{t-1})}{\partial \alpha_{t-1}} = y_{t-1} - \lambda_{t-1}. \] :

\[ y_t \mid \lambda_t \sim \text{Poisson}(\lambda_t), \qquad \lambda_t = \exp(\alpha_t), \] La función sim_gas_poisson() implementa la simulación de un modelo GAS–Poisson de primer orden. En ella, el estado dinámico \(\alpha_t\) se actualiza recursivamente a partir de su valor pasado y del score de la verosimilitud Poisson, el cual mide la discrepancia entre la observación y su valor esperado. La intensidad \(\lambda_t\) se obtiene mediante un enlace logarítmico, garantizando su positividad. Adicionalmente, se incorpora una cota al estado para evitar problemas de inestabilidad numérica asociados a la función exponencial. La función retorna la serie observada, el estado dinámico y la intensidad correspondiente.

sim_gas_poisson <- function(n, omega, A, B, alpha0 = 0) {
  
  alpha <- numeric(n)
  lambda <- numeric(n)
  y <- numeric(n)
  score <- numeric(n)
  
  alpha[1] <- alpha0
  lambda[1] <- exp(alpha[1])
  y[1] <- rpois(1, lambda[1])
  score[1] <- y[1] - lambda[1]
  
  for (t in 2:n) {
    alpha[t] <- omega + A * score[t-1] + B * alpha[t-1]
    
    # Protección numérica
    alpha[t] <- max(min(alpha[t], 20), -20)
    
    lambda[t] <- exp(alpha[t])
    y[t] <- rpois(1, lambda[t])
    score[t] <- y[t] - lambda[t]
  }
  
  list(y = y, alpha = alpha, lambda = lambda)
}

Log-verosimilitud (con reparametrización)

La función loglik_gas_poisson() implementa el cálculo de la log-verosimilitud del modelo GAS–Poisson. Los parámetros se reparametrizan para asegurar positividad del coeficiente asociado al score y estabilidad de la dinámica del estado. A partir de una inicialización robusta del estado, la función reconstruye recursivamente la trayectoria de \(\alpha_t\) y evalúa la verosimilitud condicional Poisson en cada período. Se incorporan cotas al estado y verificaciones de finitud en la log-verosimilitud para prevenir inestabilidades numéricas durante la optimización.

loglik_gas_poisson <- function(par, y) {
  
  omega <- par[1]
  A <- exp(par[2])      # A > 0
  B <- tanh(par[3])     # |B| < 1
  
  n <- length(y)
  alpha <- numeric(n)
  lambda <- numeric(n)
  score <- numeric(n)
  
  # Inicialización robusta
  alpha[1] <- log(mean(y) + 1e-6)
  lambda[1] <- exp(alpha[1])
  score[1] <- y[1] - lambda[1]
  
  loglik <- dpois(y[1], lambda[1], log = TRUE)
  
  for (t in 2:n) {
    alpha[t] <- omega + A * score[t-1] + B * alpha[t-1]
    alpha[t] <- max(min(alpha[t], 20), -20)
    
    lambda[t] <- exp(alpha[t])
    score[t] <- y[t] - lambda[t]
    
    ll <- dpois(y[t], lambda[t], log = TRUE)
    if (!is.finite(ll)) return(1e10)
    
    loglik <- loglik + ll
  }
  
  return(-loglik)
}

Estimación por Máxima Verosimilitud

La función estimate_gas_poisson_mc() utiliza esta log-verosimilitud para estimar los parámetros mediante máxima verosimilitud usando el algoritmo BFGS. En el contexto del estudio Monte Carlo, se reportan únicamente las estimaciones puntuales, evitando el uso del Hessiano, el cual puede resultar no invertible en muestras finitas para modelos GAS no lineales.

estimate_gas_poisson_mc <- function(y) {
  
  init <- c(
    log(mean(y) + 1e-6),
    log(0.1),
    atanh(0.7)
  )
  
  fit <- optim(
    par = init,
    fn = loglik_gas_poisson,
    y = y,
    method = "BFGS",
    control = list(maxit = 300)
  )
  
  c(
    omega = fit$par[1],
    A = exp(fit$par[2]),
    B = tanh(fit$par[3])
  )
}

Filtrado GAS (estado estimado)

Se implementa el filtrado del modelo GAS–Poisson utilizando los parámetros estimados. A partir de la serie observada, se reconstruye recursivamente el estado dinámico mediante la ecuación GAS, obteniendo tanto la trayectoria del estado \(\alpha_t\) como la intensidad \(\lambda_t\). Este procedimiento corresponde al filtrado determinístico propio de los modelos GAS, donde el estado se actualiza en función del score y no mediante un término de ruido latente.

filter_gas_poisson <- function(y, par) {
  
  omega <- par[1]
  A <- par[2]
  B <- par[3]
  
  n <- length(y)
  alpha <- numeric(n)
  lambda <- numeric(n)
  score <- numeric(n)
  
  alpha[1] <- log(mean(y) + 1e-6)
  lambda[1] <- exp(alpha[1])
  score[1] <- y[1] - lambda[1]
  
  for (t in 2:n) {
    alpha[t] <- omega + A * score[t-1] + B * alpha[t-1]
    alpha[t] <- max(min(alpha[t], 20), -20)
    
    lambda[t] <- exp(alpha[t])
    score[t] <- y[t] - lambda[t]
  }
  
  list(alpha = alpha, lambda = lambda)
}

Ejecución completa con n = 100

Los estimadores obtenidos bajo el modelo GAS–Poisson son consistentes con los valores verdaderos utilizados en la simulación y comparables con los resultados obtenidos previamente. En particular, el parámetro de persistencia 𝐵 B se mantiene elevado, confirmando una dinámica altamente dependiente del pasado, mientras que el parámetro 𝐴 A presenta una magnitud similar a la observada en las estimaciones anteriores, reflejando una respuesta estable del estado al score. En conjunto, los resultados muestran que el modelo GAS–Poisson captura adecuadamente la dinámica de la serie y ofrece un ajuste coherente con los enfoques considerados previamente.

set.seed(123)

# Parámetros verdaderos
omega <- 0
A <- 0.5
B <- 0.8

# Simulación
sim <- sim_gas_poisson(100, omega, A, B)

# Estimación
est <- estimate_gas_poisson_mc(sim$y)
est
##      omega          A          B 
## 0.01126156 0.18847294 0.76666193

Gráficos

Interpretación de los resultados

La serie observada \[y_t\] presenta una alta proporción de ceros y algunos episodios aislados de valores positivos, lo cual es coherente con el supuesto Poisson del modelo. El estado latente estimado \[\hat{\alpha}_t\] sigue de cerca al estado verdadero \[\alpha_t\], capturando adecuadamente su dinámica general, aunque con un leve suavizamiento, especialmente en los cambios bruscos. Esto es esperable dado que el estado estimado se basa únicamente en la información observada. Finalmente, la intensidad

\[ \lambda_t = \exp(\alpha_t) \] refleja correctamente los períodos de mayor actividad, mostrando que el modelo GAS–Poisson logra representar de forma adecuada la evolución temporal de la media condicional del proceso.

filtered <- filter_gas_poisson(sim$y, est)

par(mfrow = c(3,1))

plot(sim$y, type = "h", main = "Serie observada y_t", ylab = "y")

plot(sim$alpha, type = "l", col = "red", lwd = 2,
     main = "Estado verdadero vs estimado",
     ylab = expression(alpha))
lines(filtered$alpha, col = "blue", lwd = 2)
legend("topleft", c("Verdadero", "Estimado"),
       col = c("red", "blue"), lwd = 2)

plot(sim$lambda, type = "l", col = "darkgreen", lwd = 2,
     main = expression(lambda[t] == exp(alpha[t])),
     ylab = expression(lambda))

#Base para Monte Carlo

Resultados del experimento Monte Carlo

El experimento Monte Carlo evidencia comportamientos diferenciados entre los parámetros del modelo GAS–Poisson. Las estimaciones de los parámetros dinámicos \[A\] y \[B\] presentan una dispersión moderada y valores medios relativamente cercanos a los parámetros verdaderos, lo que indica que la dinámica del estado es razonablemente bien identificada incluso con muestras de tamaño \[n = 100\].

En contraste, el parámetro \[\omega\] muestra una alta variabilidad y una distribución claramente asimétrica, con valores extremos negativos que se reflejan tanto en el boxplot como en su elevada desviación estándar. Esto sugiere que \[\omega\] está débilmente identificado en muestras pequeñas y que su estimación es particularmente sensible a realizaciones con muchos ceros en la serie observada, un fenómeno típico en modelos Poisson.

monte_carlo_gas <- function(R, n, omega, A, B) {
  
  est_mat <- matrix(NA, R, 3)
  colnames(est_mat) <- c("omega", "A", "B")
  
  for (r in 1:R) {
    sim <- sim_gas_poisson(n, omega, A, B)
    est_mat[r, ] <- estimate_gas_poisson_mc(sim$y)
  }
  
  est_mat
}

Ejemplo

Resultados del experimento Monte Carlo

En conjunto, los resultados indican que el modelo GAS–Poisson estima de forma más estable los parámetros que controlan la dinámica del estado latente (\[A\] y \[B\]), mientras que el término constante \[\omega\] requiere tamaños muestrales mayores para lograr una estimación precisa.

set.seed(123)

mc <- monte_carlo_gas(100, 100, omega, A, B)

apply(mc, 2, mean)
##         omega             A             B 
## -1382.3614115     0.2572194     0.4962236
apply(mc, 2, sd)
##        omega            A            B 
## 6185.7984236    0.2171619    0.6206364
boxplot(mc)

# Monte Carlo para n = 500
mc_500 <- monte_carlo_gas(100, 500, omega, A, B)

apply(mc_500, 2, mean)
##         omega             A             B 
## -4.846944e+03  5.728693e-02  2.832043e-01
apply(mc_500, 2, sd)
##        omega            A            B 
## 2.988339e+04 1.025278e-01 8.972898e-01
boxplot(mc_500, main = "Monte Carlo GAS–Poisson (n = 500)")

# Monte Carlo para n = 1000
mc_1000 <- monte_carlo_gas(100, 1000, omega, A, B)

apply(mc_1000, 2, mean)
##         omega             A             B 
## -3587.0778012     0.1194341     0.5001461
apply(mc_1000, 2, sd)
##        omega            A            B 
## 1.691320e+04 5.099234e-01 8.146807e-01
boxplot(mc_1000, main = "Monte Carlo GAS–Poisson (n = 1000)")

Comparación Monte Carlo para distintos tamaños muestrales

Al aumentar el tamaño muestral de \(n = 100\) a \(n = 500\) y \(n = 1000\), las estimaciones de los parámetros dinámicos \(A\) y \(B\) muestran una mejora progresiva en su comportamiento, especialmente en el caso de \(B\), cuya media se aproxima al valor verdadero y cuya dispersión se reduce parcialmente. Esto indica que la persistencia del estado latente es mejor identificada a medida que aumenta la información disponible en la serie.

En contraste, el parámetro \(\omega\) presenta una elevada variabilidad y valores extremos incluso para tamaños muestrales grandes. Este comportamiento se debe a que \(\omega\) actúa como un término de nivel en la ecuación del estado y está débilmente identificado en modelos GAS–Poisson cuando la serie observada contiene una alta proporción de ceros, lo que es característico de este tipo de procesos.

Los boxplots confirman que, aunque el aumento del tamaño muestral mejora parcialmente la estabilidad de las estimaciones, el parámetro \(\omega\) sigue siendo el más problemático, mientras que los parámetros \(A\) y \(B\), asociados a la dinámica del score y a la persistencia, presentan un comportamiento más robusto y estable.

par(mfrow = c(1,3))
boxplot(mc, main = "n = 100")
boxplot(mc_500, main = "n = 500")
boxplot(mc_1000, main = "n = 1000")

par(mfrow = c(1,1))