Estimar un modelo de regresión lineal univariada con muestreo de Gibbs (estadística bayesiana)

 1.1.           Simulando datos

Primero simulamos los datos, asumiendo un beta1 real de 0.5. Puntualmente se generan dos variables aleatorias x e y, con una relación de 0.5, más un término de error. El código se desarrolló en el marco de un curso sobre “Modelos para el análisis de políticas macroeconómicas”. En esta sección se puede indicar que BETA y SIG2 son los parámetros a estimar en el procedimiento. Como créditos se aclara que el código original en Matlab fue elaborado por el profesor Daniel Parra.

rm(list=ls())
set.seed(1)
 

# Generar datos
# Modelo
# y_{t} = BETA * x_{t} + v_{t}
# v_{t} ~ N(0,SIG2)
 
BETA <- 0.5
SIG2 <- 0.006 # Varianza del error
k <- 1  # Número de variables explicativas
t <- 500  # Longitud de la muestra
e1 <- rnorm(t, mean = 0, sd = sqrt(SIG2))  # Generar errores distribuidos N(0, SIG2)
x <- rnorm(t, mean = 0, sd = 1)  # Variable explicativa: extracciones i.i.d. de N(0,1)
 
y <- x * BETA + e1 # Calcular la variable dependiente

1.2. Enfoque frecuentista

Utiizando el enfoque tradicional del modelo de regresión podemos recuperar, a partir de las variables simuladas, los coeficientes a estimar –tenga pendiente que estos procesos están automatizados en R mediante el comando lm tradicional (lm(y~x))-:

 # Estimación de OLS (para comparación)
bols <- solve(t(x) %*% x) %*% t(x) %*% y
sigols <- t(y - x * bols) %*% (y - x * bols) / (t - k)

Observando los resultados se puede apreciar como el enfoque logra aproximar bien los datos.

        coef  real         mco
1      beta 0.500 0.496877122
2 sig.error 0.006 0.006136154

1.3. Enfoque bayesiano

En primer lugar, es importante recordar que el enfoque bayesiano está interesado en estimar una posterior que resulta del producto de la prior y una función de verosimilitud. La prior se estima o supone basado en el conocimiento previo del analista, mientras que la mle debe estimarse mediante algún algoritmo. En el siguiente ejemplo se utiliza el algoritmo de Gibbs. En esta línea se asume una distribución normal del beta y una gamma del error del modelo y se asumen parámetros de estas definiciones, como la media del set de parámetros que se estimaran mediante la estimación bayesiana.

b0 <- 0.2  # Media de la distribución a priori para beta
P0 <- 10   # Varianza de la distribución a priori para beta
 
t0 <- 0  # Parámetro de forma a priori para sigma
R0 <- 0  # Parámetro de escala a priori para sigma: cero implica un prior "plano"
# es decir, sin información a priori

Posteriormente se establecen los parámetros para el rastreador de Gibbs. La idea de estos es tomar x cantidad de simulaciones, pero omitir las primeras dado que estas pueden estar incididas por las condiciones de inicio establecidas en el algoritmo. Puntualmente, se generan 100,000 simulaciones (n), pero de estas, posteriormente se desechan (Período de quema: las extracciones se descartan) las primeras 90,000. Por tanto, es sobre las ultimas 10 extracciones sobre las que se basara la inferencia. Note el condicional if que se encuentra en el bucle (for) de abajo, este guarda solo las últimas 10,000 iteracciones.

nburn <- 90000  # Período de quema: las extracciones se descartan
nkeep <- 10000  # Número de extracciones en las que se basará la inferencia
n <- nkeep + nburn  # Número total de iteraciones del muestreador de Gibbs

Luego se generan los objetos donde se guardará la distribución de betas y errores (parámetros) que serán obtenidas a partir de Gibbs.

sig2 <- 1  # Valor inicial para el muestreador de Gibbs
 
Beta <- matrix(0, nkeep, k)
Sig2 <- numeric(nkeep)

Finalmente se utiliza el algoritmo de Gibbs para obtener la distribución de parámetros, y en función de la media de esta distribución el valor estimado. Puntualmente, aquí se define un algoritmo iterativo de las simulaciones, desde 1 hasta n. En el primer paso se utiliza el prior para estimar el beta, más una corrección conocida como factor de incertidumbre. Luego con este beta, estimado en la primera iteración se utiliza para estimar la varianza. En tal sentido, sig2 guarda un posible candidato aleatorio para la varianza. Como la varianza no puede ser negativa se usa una distribución gamma. posterior usa la información del prior, con la información que recibo de los datos.

for (iter in 1:n) {
 
  # Generar una extracción para beta condicional a sig2
  P1 <- solve(solve(P0) + (1/sig2) * t(x) %*% x)
  b1 <- P1 %*% (solve(P0) %*% b0 + (1/sig2) * t(x) %*% y)
  beta <- b1 + chol(P1) %*% rnorm(k, mean = 0, sd = 1)
 
  # Generar una extracción para sig2 condicional a beta
  t1 <- t0 + t
  R1 <- R0 + t(y - x * beta) %*% (y - x * beta)
  shpe <- t1
  scle <- 1 / R1
  sig2inv <- rgamma(1, shape = shpe, rate = scle)
  sig2 <- 1 / sig2inv
 
  if (iter > nburn) {
    Beta[iter - nburn,] <- beta
    Sig2[iter - nburn] <- sig2
  }
}

Estos pasos están relacionados con la actualización de parámetros en el muestreo de Gibbs.

 -          P0 es la matriz de varianza-covarianza a priori para el parámetro beta.

-          sig2 es la varianza del error en el modelo.

-        P1 es la inversa de la suma de la inversa de la matriz prior P0 y el producto (1/sig2) * (x'*x). se utiliza para actualizar la matriz de covarianza-condicional de la distribución posterior de los parámetros beta dados los valores actuales de otros parámetros y los datos observados. La actualización se realiza en la fase de muestreo de Gibbs, donde se obtienen nuevas muestras de los parámetros del modelo. La matriz resultante P1 es la varianza-covarianza condicional de la distribución posterior de beta dada la varianza del error sig2 y los datos observados.

 - b0 es el valor medio a priori para el parámetro beta. Esta línea actualiza el vector de medias condicionales de la distribución posterior de beta dada la varianza del error sig2 y los datos observados.

-y es el vector de observaciones de la variable dependiente.

-          (1/sig2) * (x'*y) es la contribución de los datos observados al término de actualización.

En resumen, estas líneas de código están realizando la actualización de parámetros en el contexto del muestreo de Gibbs para obtener una nueva muestra de la distribución conjunta de beta condicional en los valores actuales de los demás parámetros y los datos observados.

 1.4.           Resultados

En primer lugar, se observa que pese a tener un prior algo alejado del verdadero valor del beta, el algoritmo es capaz de empujar sus resultados hacia el valor verdadero.

hist(Beta, main = expression(beta), col = "lightblue", xlim = c(0.32, 0.6))
abline(v = BETA, col = "red", lwd = 2)
abline(v = b0, col = "black", lwd = 2)
legend("topright", legend = c("Estimado", "Verdadero", "Prior"), col = c("blue", "red", "black"), lty = 1, cex = 0.8)

Posteriormente, podemos comparar la distribución obtenida, vs. El valor verdadero del parámetro, y las estimaciones mco.

 hist(Beta, main = expression(beta), col = "lightblue")
abline(v = BETA, col = "red", lwd = 2)
abline(v = mean(Beta), col = "green", lwd = 3)
abline(v = bols, col = "black", lwd = 2)
legend("topright", legend = c("Estimado", "Verdadero", "Media posterior", "MCO"), col = c("blue", "red", "green"), lty = 1, cex = 0.8)

 Y las tablas de resultados:

        coef  real         mco       gibbs
1      beta 0.500 0.499524587 0.499512621
2 sig.error 0.006 0.006116056 0.006568336