Simulación, calibración y estimación: ejemplos en AR(1) en R

 Las palabras simular, estimar, calibrar e identificar, aparecen de forma recurrente en la economía aplicada. Estas representan una extensión de la forma de pensar y el desarrollo de aplicaciones teóricas.  Es frecuente que la identificación y calibración se tomen como sub conjuntos de la estimación. Suponga observamos una serie y(t), sobre la que suponemos sigue un proceso ar(1). Por ejemplo:

 y(t)=phi_1* y(t-1)+ e(t)  e(t)~N(0,sigma)

Podemos a partir de ciertos supuestos simular una serie que cumpla con estas características, con el objetivo de estudiar sus propiedades o su comportamiento bajo ciertas modificaciones. Un caso distinto es intentar recuperar los coeficientes phi_1 y sigma de un AR(1), aquí entra la estimación, que permite obtenerlos a partir de algún procedimiento estadístico como el de Mínimo Cuadrados Ordinarios (MCO) o máxima verosimilitud (mle). Mientras que la calibración, también pretende recuperar los coeficientes del modelo, pero a partir de suposiciones o conocimientos previos, buscando que el valor calibrado permita obtener una variable y*(t) parecida al valor observado. 

 1.       Simular

 En la siguiente entrada se simula en R un autoregresivo ar(1):

 y(t)=0.8* y(t-1)+ e(t)  e(t)~N(0,1)

 Se muestra que la simulación de series temporales, y entre calibrar y estimar un modelo. “Simular datos implica generar una serie temporal sintética a partir de un modelo AR(1) con parámetros conocidos” (ChatGPT, 2024).

 # Cargar las librerías necesarias

library(ggplot2)
library(dplyr)
library(tidyr)
library(stats4)
 
# Simular ---------------------
# Parámetros del modelo AR(1)
phi <- 0.8
sigma <- 1
n <- 100
 
# Simular la serie temporal
set.seed(2)
y <- numeric(n)
y[1] <- rnorm(1, mean = 0, sd = sigma)
 
for (t in 2:n) {
  y[t] <-  phi * y[t-1] + rnorm(1, mean = 0, sd = sigma)
}
 
# Crear una secuencia de fechas
fechas <- seq.Date(from = as.Date("2000-01-01"), by = "days", length.out = n)
 
# Crear un data frame con la serie temporal simulada
data <- data.frame(fecha = fechas, y = y)
 
ggplot(data, aes(fechas, y)) +
  geom_line() +
  theme_minimal()

2.       Calibrar

“Calibrar un modelo implica ajustar los parámetros del modelo basándose en conocimientos previos o estimaciones empíricas sin realizar un proceso formal de estimación. Es una forma de "ajustar" el modelo para que sea razonable para la situación específica. La calibración se refiere ajustar los parámetros del modelo basándose en conocimientos previos o estimaciones empíricas sin un proceso formal de estimación. Utiliza datos observados pero no sigue un método estadístico forma” (ChatGPT, 2024).

En la práctica, consiste en seleccionar un conjunto de parámetros de preferencia a partir de las referencias disponibles en la teoría económica, con la intención de que los modelos imiten alguna característica histórica de los datos (Hoover, 1995). Es decir, “La calibración se realiza con mayor frecuencia probando distintos valores para los parámetros hasta que el modelo logre predicciones con la menor desviación posible de los datos o reproduzca otras características empíricas” (Stackexchange, 2020). Los modelos calibrados son una herramienta importante de la investigación económica, esta herramienta ha surgido como evolución y extensión de la capacidad de los economistas para teorizar (Cooley, 1997).

# Serie temporal observada (ejemplo)
y_obs <- y
 
# Calibración: Asumir un valor razonable para phi
phi_calibrado <- 0.8
 
phi_calibrado <- 0.8
y_calibrado <- numeric(n)
y_calibrado[1] <- y[1]
 
for (t in 2:n) {
  y_calibrado[t] <- phi_calibrado * y_calibrado[t-1] + rnorm(1, mean = 0, sd = sigma)
}
 
data$y_calibrado <- y_calibrado
 
ggplot(data) +
  geom_line(aes(fechas, y, color = "Obs")) +
  geom_line(aes(fechas, y_calibrado, color = "Calibrado")) +
  theme_minimal()

 

Algunas recomendaciones identificadas en la literatura:

 ·         No justificar el uso solo porque otros autores lo han utilizado (Cooley, 1997).

 3.       Estimar

 La estimación toma los datos como dados y los utilizan para determinar los parámetros del modelo a partir de un estadístico (que es una función de la muestra), buscando que este sea insesgado, consistente y eficiente. En este procedimiento, el valor de un estadístico se utiliza para estimar un parámetro, por tanto, la estimación “utiliza métodos estadísticos formales (como MLE) para determinar los parámetros del modelo a partir de los datos observados. Es un enfoque más riguroso que la calibración y proporciona estimaciones con propiedades estadísticas deseables” (ChatGPT, 2024). En otras palabras, se utiliza un procedimiento estadístico formal para la estimación de los coeficientes del modelo.

 # Estimación usando MLE
library(stats4)
 
# Log-verosimilitud de un modelo AR(1)
logLik_AR1 <- function(phi, sigma) {
  n <- length(y_obs)
  eps <- numeric(n)
  eps[1] <- y_obs[1] # inicializamos el primer error
 
  for (t in 2:n) {
    eps[t] <- y_obs[t] - phi * y_obs[t-1]
  }
 
  logLik <- -0.5 * n * log(2 * pi * sigma^2) - sum(eps^2) / (2 * sigma^2)
  return(logLik)
}
 
# Función de verosimilitud negativa para optimización
negLogLik_AR1 <- function(phi, sigma) {
  -logLik_AR1(phi, sigma)
}
 
# Optimizar para encontrar los parámetros MLE
mle_fit <- mle(negLogLik_AR1, start = list(phi = 0.5, sigma = 1))
 
# Resultados de la estimación
summary(mle_fit)
 
Maximum likelihood estimation
 
Call:
mle(minuslogl = negLogLik_AR1, start = list(phi = 0.5, sigma = 1))
 
Coefficients:
       Estimate Std. Error
phi   0.7218984 0.06979479
sigma 1.1476188 0.08114858
 
-2 log L: 311.3256

# -------------------------------------------------------------------------------------

# Obtener los parámetros estimados
phi_estimado <- coef(mle_fit)["phi"]
sigma_estimado <- coef(mle_fit)["sigma"]
 
# Generar la serie temporal con los parámetros estimados
y_estimado <- numeric(n)
y_estimado[1] <- y[1]
 
for (t in 2:n) {
  y_estimado[t] <- phi_estimado * y_estimado[t-1] + rnorm(1, mean = 0, sd = sigma_estimado)
}
 
data$y_estimado <- y_estimado
 
# Graficar los datos simulados, calibrados y estimados
data_long <- data %>%
  pivot_longer(cols = c(y, y_calibrado, y_estimado), names_to = "serie", values_to = "valor")
 
ggplot(data_long, aes(x = fecha, y = valor, color = serie)) +
  geom_line(size = 0.7) +
  theme_minimal() +
  labs(title = "Simulación, Calibración y Estimación de un Modelo AR(1)",
       x = "Fecha", y = "Valor",
       color = "Serie") +
scale_color_manual(values = c("y" = "black", "y_calibrado" = "green", "y_estimado" = "red")) +
  theme(plot.title = element_text(hjust = 0.5),
        legend.position = "bottom")

Bibliography

Cooley, T. (1997). Calibrated Models. University of Rochester.

Hoover, K. (1995). Facts and artifacts: calibration and the empirical assessment of real-business-cycle models. Oxford Economic Papers 47 (1), 24-44.

 

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

 

Usar la función summarise+apply en dplyr

En la siguiente entrada, se simula una data, que contienen diferentes tipos de series, con el objetivo de mostrar, como al combinar las funciones summarise (del paquete dplyr) y la función apply (base de R), podemos obtener resúmenes estadísticos de un conjunto de variables, incluyendo test formales como el Dickey Fuller, o estadísticos como el coeficiente de autocorrelación. Puntualmente, tentemos una base de datos con 30 variables y deseamos obtener una tabla donde a cada variable le calculemos el promedio y la desviación estándar.

library(dplyr)
library(moments)
library(tseries)

Generamos de forma aleatoria una data llamada my_data, esta es un ejemplo de datos aleatorios para ilustrar el proceso.

set.seed(123)
my_data <- data.frame(
  serie1 = rnorm(100),
  serie2 = rnorm(100),
  fecha = as.Date("2024-01-01") + 0:99,
  serie3 = rnorm(100),
  categoria = sample(c("A", "B", "C"), 100, replace = TRUE)
)
 
          serie1      serie2      fecha      serie3 categoria
1   -0.560475647 -0.71040656 2024-01-01  2.19881035         B
2   -0.230177489  0.25688371 2024-01-02  1.31241298         A
3    1.558708314 -0.24669188 2024-01-03 -0.26514506         C
4    0.070508391 -0.34754260 2024-01-04  0.54319406         C
5    0.129287735 -0.95161857 2024-01-05 -0.41433995         C
…

Ahora creamos un filtro para seleccionar únicamente las series numéricas dentro de mi base de datos, dato que puede contener fechas y caracteres. Es decir, las series 1,2 y 3.

numeric_columns <- my_data %>%
  select_if(is.numeric) %>%
  colnames()
 
umeric_columns
[1] "serie1" "serie2" "serie3"

Ahora, usamos select para recuperar las variables numéricas y luego a cada columna de aplicamos la opción deseada combinando summarise y apply como se muestra en los siguientes ejemplos.

my_data %>%
  dplyr::select(all_of(numeric_columns)) %>%
  summarise(
    varianza = apply(., 2, var),
    media = apply(., 2, mean)
  )
 
varianza       media
1 0.8332328  0.09040591
2 0.9350631 -0.10754680
3 0.9022700  0.12046511

Finalmente, calculamos para cada columna (numérica, porque solo va a trabajar con las condiciones impuestas por alguna operación lógica), se estimarán los estadísticos indicados, y los organizara en un cuadro. Usamos la función select del paquete dplyr, esta se indica el nombre del paquete, dado que puede coincidir con funciones de otros paquetes y generar errores. Además, usamos mutate, para colocar el nombre en cada fila, permitiendo la lectura de los datos.

Error in select(., all_of(numeric_columns)) :
  unused argument (all_of(numeric_columns))

Posteriormente, se solicita un resumen de la data, pero realizada a cada una de las columnas usando apply. En el caso de funciones menos directas, como el p valor del test DF, se inserta una fusión que permita extraerlo.  

stats <- my_data %>%
  dplyr::select(all_of(numeric_columns)) %>%
  summarise(
    varianza = apply(., 2, var),
    asimetria = apply(., 2, skewness),
    curtosis = apply(., 2, kurtosis),
    autocorr_lag1 = apply(., 2, function(x) {
      acf(x, plot = FALSE)$acf[2]
    }),
    p_valor_dickey_fuller = apply(., 2, function(x) {
      adf.test(x)$p.value
    })
  )
 
# Imprimir las estadísticas
print(stats)

Creando variables por grupos en dplyr (group_by + mutate)

 Simulemos una base de hogares, donde se identifica el hogar, el sexo (1 mujer) y provincia y edad para cada miembro.
 
# Definir la lista de relaciones familiares por hogar
hogar <- c(1,1,1,1,2,2,2,3,4,4,4)
 prov <- c("a","a","a","a","c","c","c","a","b","b","b")
   id <- 1:length(hogar)
mujer <- c(1,0,1,1,0,0,2,0,1,0,0)
 edad <- c(28,29,10,8,33,27,13,24,32,11,6)
 
hogar_data <- data.frame(id, hogar, prov, mujer, edad)
hogar_data
  id hogar prov mujer edad
1   1     1    a     1   28
2   2     1    a     0   29
3   3     1    a     1   10
4   4     1    a     1    8
5   5     2    c     0   33
6   6     2    c     0   27
7   7     2    c     2   13
8   8     3    a     0   24
9   9     4    b     1   32
10 10     4    b     0   11
11 11     4    b     0    6
 
Si combinanos group_by con summarise, obtendríamos un collapse de la base de datos, con la edad promedio del grupo indicado en group_by. Es decir, pasamos a tener tantas observaciones como hogares en nuestra base de datos. Noten que en este caso solo quedan aquellas variables a las que le hemos indicado la forma de agregación.
 
hogar_data |>
  group_by(hogar) |>
  summarise(mean(edad))
 
# A tibble: 4 x 2
  hogar `mean(edad)`
  <dbl>        <dbl>
1     1         18.8
2     2         24.3
3     3         24 
4     4         16.3
 
Ahora, en lugar de obtener el resumen anterior, quisiéramos agregar una variable que indique la edad promedio por hogar. Podemos combinar group_by + mutate. Ahora agregamos una variable con la edad promedio de cada hogar.
 
hogar_data |>
  group_by(hogar) |>
  mutate(edad_mhogar = mean(edad))
 
# A tibble: 11 x 6
# Groups:   hogar [4]
      id hogar prov  mujer  edad edad_mhogar
   <int> <dbl> <chr> <dbl> <dbl>       <dbl>
 1     1     1 a         1    28        18.8
 2     2     1 a         0    29        18.8
 3     3     1 a         1    10        18.8
 4     4     1 a         1     8        18.8
 5     5     2 c         0    33        24.3
 6     6     2 c         0    27        24.3
 7     7     2 c         2    13        24.3
 8     8     3 a         0    24        24 
 9     9     4 b         1    32        16.3
10    10     4 b         0    11        16.3
11    11     4 b         0     6        16.3
 
Contar número de menores de edad en cada hogar. Contar la cantidad de menores de edad en cada hogar.
 
hogar_data |>
  group_by(hogar) |>
  mutate(num_menores = sum(edad < 18))
 
# A tibble: 11 x 6
# Groups:   hogar [4]
      id hogar prov  mujer  edad num_menores
   <int> <dbl> <chr> <dbl> <dbl>       <int>
 1     1     1 a         1    28           2
 2     2     1 a         0    29           2
 3     3     1 a         1    10           2
 4     4     1 a         1     8           2
 5     5     2 c         0    33           1
 6     6     2 c         0    27           1
 7     7     2 c         2    13           1
 8     8     3 a         0    24           0
 9     9     4 b         1    32           2
10    10     4 b         0    11           2
11    11     4 b         0     6           2
 
Contar el numero de menores de edad que son mujeres.
 
hogar_data |>
  group_by(hogar) |>
  mutate(num_menores = sum(edad < 18 & mujer ==1))
 
# A tibble: 11 x 6
# Groups:   hogar [4]
      id hogar prov  mujer  edad num_menores
   <int> <dbl> <chr> <dbl> <dbl>       <int>
 1     1     1 a         1    28           2
 2     2     1 a         0    29           2
 3     3     1 a         1    10           2
 4     4     1 a         1     8           2
 5     5     2 c         0    33           0
 6     6     2 c         0    27           0
 7     7     2 c         2    13           0
 8     8     3 a         0    24           0
 9     9     4 b         1    32           0
10    10     4 b         0    11           0
11    11     4 b         0     6           0