Modelo de Espacio-Estado: Modelo de Nivel Local

En este trabajo se estudia el Modelo de Nivel Local dentro del marco general de los modelos de espacio-estado. Este tipo de modelos permite representar procesos dinámicos en los que una parte de la información es observable y otra corresponde a un estado latente que evoluciona en el tiempo. En particular, se analiza un modelo en el cual las observaciones son la suma de un nivel no observado y un error aleatorio, mientras que dicho nivel sigue un proceso de paseo aleatorio.

Para la estimación del estado latente se utiliza el filtro de Kalman, un algoritmo recursivo que combina predicción y actualización para obtener estimaciones óptimas en el sentido de mínima varianza cuadrática. Además, los parámetros del modelo se estiman mediante el método de máxima verosimilitud, aprovechando que el propio filtro entrega los residuos y varianzas necesarias para calcular la verosimilitud del modelo.

Ecuación de observación

\[ y_t = \alpha_t + \varepsilon_t, \qquad \varepsilon_t \sim N(0, \sigma_\varepsilon^2) \]

Ecuación de estado

\[ \alpha_t = \alpha_{t-1} + \eta_t, \qquad \eta_t \sim N(0, \sigma_\eta^2) \]

Simulación:

1) Librerias, Fijar semilla y definir los parámetros verdaderos del modelo de nivel local

library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.5.1
library(KFAS)
## Warning: package 'KFAS' was built under R version 4.5.2
## Please cite KFAS in publications by using: 
## 
##   Jouni Helske (2017). KFAS: Exponential Family State Space Models in R. Journal of Statistical Software, 78(10), 1-39. doi:10.18637/jss.v078.i10.
set.seed(2025)

sigma_eps <- 1   # √1 = 1
sigma_eta <- 2   # √4 = 2

##Simulacion n =100

Generamos una trayectoria del estado verdadero \(α_t\) y posteriormente las observaciones \(y_t\)

n1 <- 100
alpha_100 <- numeric(n1)  # Estado oculto
y_100 <- numeric(n1)      # Observaciones

# Valor inicial del estado (puede asumirse α_1 = 0 sin pérdida de generalidad)

alpha_100[1] <- 0

# Generamos un paseo aleatorio, que corresponde a la evolución del estado
for (t in 2:n1) {
  alpha_100[t] <- alpha_100[t-1] + rnorm(1, 0, sigma_eta)
}

# Agregamos el ruido de observación
y_100 <- alpha_100 + rnorm(n1, 0, sigma_eps)

##Simulacion n =500 Lo mismo que antes, pero ahora con más observaciones. El objetivo es estudiar cómo cambia la calidad de la estimación al aumentar n.

n2 <- 500
alpha_500 <- numeric(n2)
y_500 <- numeric(n2)

alpha_500[1] <- 0
for (t in 2:n2) {
  alpha_500[t] <- alpha_500[t-1] + rnorm(1, 0, sigma_eta)
}

y_500 <- alpha_500 + rnorm(n2, 0, sigma_eps)

# Convertir a series de tiempo
y_ts_100 <- ts(y_100)
y_ts_500 <- ts(y_500)
alpha_ts_100 <- ts(alpha_100)
alpha_ts_500 <- ts(alpha_500)

Gráficos que permiten comparar visualmente el estado real vs las observaciones.

Se aprecia que las observaciones son más ruidosas que el estado, lo cual es coherente con la estructura del modelo.

# Gráfico para n=100
plot(y_ts_100, type = "l", col = "blue", 
     main = "Simulación n=100", 
     ylab = "Valor", xlab = "Tiempo",
     ylim = range(c(y_100, alpha_100)))
lines(alpha_ts_100, col = "red", lwd = 2)
legend("topleft", 
       legend = c("Observaciones (y)", "Estado real (α)"),
       col = c("blue", "red"), lty = 1)

# Gráfico para n=500
plot(y_ts_500, type = "l", col = "blue",
     main = "Simulación n=500",
     ylab = "Valor", xlab = "Tiempo",
     ylim = range(c(y_500, alpha_500)))
lines(alpha_ts_500, col = "red", lwd = 2)
legend("topleft",
       legend = c("Observaciones (y)", "Estado real (α)"),
       col = c("blue", "red"), lty = 1)

##Estimacion por Maxima Verosimilitud

Aquí estimamos los parámetros del modelo usando fitSSM,que implementa el método de máxima verosimilitud basado en el filtro de Kalman.

## Función simple para estimar un modelo
estimar_modelo <- function(serie) {
  # Crear modelo
  mod <- SSModel(serie ~ SSMtrend(1, Q = NA), H = NA)
  
  # Estimar parámetros
  fit <- fitSSM(mod, inits = log(c(1, 1)))
  
  # Obtener varianzas estimadas
  # CUIDADO: fit$optim.out$par[1] es log(σ²_η)
  #          fit$optim.out$par[2] es log(σ²_ε)
  var_eta_est <- exp(fit$optim.out$par[1])  # σ²_η
  var_eps_est <- exp(fit$optim.out$par[2])  # σ²_ε
  
  # Aplicar filtro de Kalman
  kalman <- KFS(fit$model)
  
  # Estado filtrado
  estado_filtrado <- kalman$att[, 1]
  
  # Devolver resultados
  return(list(
    var_eps = var_eps_est,
    var_eta = var_eta_est,
    estado = estado_filtrado,
    modelo = fit$model
  ))
}

resultados_100 <- estimar_modelo(y_ts_100)
resultados_500 <- estimar_modelo(y_ts_500)

Resultados

cat("=== RESULTADOS PARA n=100 ===\n")
## === RESULTADOS PARA n=100 ===
cat("σ²_ε estimada:", round(resultados_100$var_eps, 4), "\n")
## σ²_ε estimada: 0.3914
cat("σ_ε estimada:", round(sqrt(resultados_100$var_eps), 4), "\n")
## σ_ε estimada: 0.6257
cat("σ²_η estimada:", round(resultados_100$var_eta, 4), "\n")
## σ²_η estimada: 4.6603
cat("σ_η estimada:", round(sqrt(resultados_100$var_eta), 4), "\n")
## σ_η estimada: 2.1588
cat("Valores verdaderos: σ²_ε=1, σ²_η=4\n\n")
## Valores verdaderos: σ²_ε=1, σ²_η=4
cat("=== RESULTADOS PARA n=500 ===\n")
## === RESULTADOS PARA n=500 ===
cat("σ²_ε estimada:", round(resultados_500$var_eps, 4), "\n")
## σ²_ε estimada: 1.1
cat("σ_ε estimada:", round(sqrt(resultados_500$var_eps), 4), "\n")
## σ_ε estimada: 1.0488
cat("σ²_η estimada:", round(resultados_500$var_eta, 4), "\n")
## σ²_η estimada: 4.0057
cat("σ_η estimada:", round(sqrt(resultados_500$var_eta), 4), "\n")
## σ_η estimada: 2.0014
cat("Valores verdaderos: σ²_ε=1, σ²_η=4\n")
## Valores verdaderos: σ²_ε=1, σ²_η=4

Gráfica del estado estimado

El objetivo es comparar visualmente cuán cerca está el estado filtrado del verdadero estado oculto.

Como podemos observar entre más datos, el filtro es más estable. Mostrando la mayor eficiencia del filtro cuando hay más informacion disponible

# Para n=100
plot(y_ts_100, type = "l", col = "gray",
     main = "Estado estimado vs Real (n=100)",
     ylab = "Valor", xlab = "Tiempo")
lines(alpha_ts_100, col = "green", lwd = 2)        # Estado real
lines(resultados_100$estado, col = "red", lwd = 2) # Estado estimado
legend("topleft",
       legend = c("Observado", "Estado real", "Estado estimado"),
       col = c("gray", "green", "red"), lty = 1)

# Para n=500
plot(y_ts_500, type = "l", col = "gray",
     main = "Estado estimado vs Real (n=500)",
     ylab = "Valor", xlab = "Tiempo")
lines(alpha_ts_500, col = "green", lwd = 2)        # Estado real
lines(resultados_500$estado, col = "red", lwd = 2) # Estado estimado
legend("topleft",
       legend = c("Observado", "Estado real", "Estado estimado"),
       col = c("gray", "green", "red"), lty = 1)

Comentario:

Note que efectivamente a mayor numero de datos, la estimacion mejor.

# Note que:
cat("1. Para n=100:\n")
## 1. Para n=100:
cat("   - σ²_ε estimada:", round(resultados_100$var_eps, 4), "(verdadera: 1)\n")
##    - σ²_ε estimada: 0.3914 (verdadera: 1)
cat("   - σ²_η estimada:", round(resultados_100$var_eta, 4), "(verdadera: 4)\n")
##    - σ²_η estimada: 4.6603 (verdadera: 4)
cat("   - Error en σ²_ε:", round(abs(1 - resultados_100$var_eps), 4), "\n")
##    - Error en σ²_ε: 0.6086
cat("   - Error en σ²_η:", round(abs(4 - resultados_100$var_eta), 4), "\n\n")
##    - Error en σ²_η: 0.6603
cat("2. Para n=500:\n")
## 2. Para n=500:
cat("   - σ²_ε estimada:", round(resultados_500$var_eps, 4), "(verdadera: 1)\n")
##    - σ²_ε estimada: 1.1 (verdadera: 1)
cat("   - σ²_η estimada:", round(resultados_500$var_eta, 4), "(verdadera: 4)\n")
##    - σ²_η estimada: 4.0057 (verdadera: 4)
cat("   - Error en σ²_ε:", round(abs(1 - resultados_500$var_eps), 4), "\n")
##    - Error en σ²_ε: 0.1
cat("   - Error en σ²_η:", round(abs(4 - resultados_500$var_eta), 4), "\n\n")
##    - Error en σ²_η: 0.0057

Función para realizar una simulacion y luego estimar el modelo:

# Función para una simulación y estimación
una_simulacion <- function(n_obs) {
  # Simular
  alpha <- cumsum(c(0, rnorm(n_obs - 1, 0, 2)))  # σ_η = 2
  y <- alpha + rnorm(n_obs, 0, 1)                # σ_ε = 1
  
  # Estimar
  mod <- SSModel(y ~ SSMtrend(1, Q = NA), H = NA)
  fit <- fitSSM(mod, inits = log(c(1, 1)))
  
  # Devolver estimaciones (CORREGIDO: [1]=σ²_η, [2]=σ²_ε)
  return(c(
    sigma2_eps = exp(fit$optim.out$par[2]),  # σ²_ε
    sigma2_eta = exp(fit$optim.out$par[1])   # σ²_η
  ))
}

Método de Monte Carlo: Repetir la estimacion muchas veces:

Este procedimiento permite estudiar el comportamiento estadístico de los estimadores del modelo: - sesgo - variabilidad - convergencia al aumentar el número de réplicas

# Hacer 100 simulaciones
cat("Haciendo 100 simulaciones...\n")
## Haciendo 100 simulaciones...
sim_100 <- matrix(NA, nrow = 100, ncol = 2)
for (i in 1:100) {
  sim_100[i, ] <- una_simulacion(100)
  if (i %% 20 == 0) cat("  ", i, "de 100\n")
}
##    20 de 100
##    40 de 100
##    60 de 100
##    80 de 100
##    100 de 100
# Hacer 500 simulaciones
cat("\nHaciendo 500 simulaciones...\n")
## 
## Haciendo 500 simulaciones...
sim_500 <- matrix(NA, nrow = 500, ncol = 2)
for (i in 1:500) {
  sim_500[i, ] <- una_simulacion(100)
  if (i %% 100 == 0) cat("  ", i, "de 500\n")
}
##    100 de 500
##    200 de 500
##    300 de 500
##    400 de 500
##    500 de 500
# Hacer 1000 simulaciones
cat("\nHaciendo 1000 simulaciones...\n")
## 
## Haciendo 1000 simulaciones...
sim_1000 <- matrix(NA, nrow = 1000, ncol = 2)
for (i in 1:1000) {
  sim_1000[i, ] <- una_simulacion(100)
  if (i %% 200 == 0) cat("  ", i, "de 1000\n")
}
##    200 de 1000
##    400 de 1000
##    600 de 1000
##    800 de 1000
##    1000 de 1000

Box Plot

Los boxplots permiten visualizar la dispersión de las estimaciones y cómo esta cambia al aumentar el número de simulaciones.

# Configurar área de gráficos
par(mfrow = c(1, 2))

# Boxplot para σ²_ε
boxplot(sim_100[, 1], sim_500[, 1], sim_1000[, 1],
        names = c("100 sim", "500 sim", "1000 sim"),
        main = "Estimaciones de σ²_ε",
        ylab = "σ²_ε",
        col = "lightblue")
abline(h = 1, col = "red", lty = 2)  # Valor verdadero

# Boxplot para σ²_η
boxplot(sim_100[, 2], sim_500[, 2], sim_1000[, 2],
        names = c("100 sim", "500 sim", "1000 sim"),
        main = "Estimaciones de σ²_η",
        ylab = "σ²_η",
        col = "lightgreen")
abline(h = 4, col = "red", lty = 2)  # Valor verdadero

# Volver a un solo gráfico simple
par(mfrow = c(1, 1))

Estimaciones de \(σ²_ε\):

  • La dispersión (altura de las cajas y longitud de los bigotes) es relativamente moderada.

  • Al comparar 100, 500 y 1000 simulaciones: Se pueden observar menores puntos atípicos, lo que indica mayor estabilidad estadística en la distribucion estimada cuadno aumnetadmos el bnúmero de réplicas.

La variabilidad es menor que en \(σ²_η\), lo cual concuerda con la teoría: el ruido de observación suele ser más fácil de estimar que el ruido del estado.

Estimaciones de \(σ²_η\):

  • Se observan valores atípicos más grandes, indicando una mayor sensibilidad del estimador a la variabilidad de la simulación.

  • Al aumentar el número de simulaciones (100 → 500 → 1000):Se observan más valores atíícos, indicando una mayor sensibilidad del estimador a la variabilidad de la simulación.

  • Con pocas simulaciones, la varianza empírica puede subestimarse o sobreestimarse. Al aumentar el número de réplicas, las varianzas de los estimadores es más estable y más confiable

En resumen, con más datos se obtienen mejores estiamciones. Con más similaciones los resultados son más estables, Cuando \(σ²_ε\) es grande, es más dificil estimarse.