---
title: null
output:
  html_document:
    toc: true
    toc_float: true
    collapsed: true
    smooth_scroll: false
    code_folding: hide
    theme: cerulean
    highlight: tango
    number_sections: true
    df_print: kable
    fig_height: 5
    fig_width: 7
    css: estilos.css
---

::: {style="text-align:center; margin-bottom:20px;"}
<img src="portada.png" alt="Portada" style="max-width:100%; height:auto;"/>
:::

::: {style="text-align:left; margin-left:40px; margin-bottom:40px;"}
<h1>Fijación ordinaria de comisiones máximas del sistema de tarjetas de pago 2025</h1>

<h4>División Económica y División Sistemas de Pago, BCCR</h4>

<h5>Noviembre 2025</h5>
:::

<hr>

```{r setup, include=FALSE}
knitr::opts_chunk$set(echo = TRUE, message = FALSE, warning = FALSE)

library(AER)#Pruebas sobre endogeneidad y estimación 
library(bootStepAIC) #Bootstrap
library(car)#Para condición de multicolinealidad
library(data.table) #Crear tabla de datos
library(dplyr) # Manipulación de datos
library(furrr)
library(future)
library(ggplot2) #Gráfico
library(glmulti)#Para selección de variables
library(leaps) # Para selección de variables
library(lme4) # modelo efectos aleatorios
library(lmtest)#Para pruebas de heterocedasticidad y autocorrelación
library(MASS) # regresión robusta
library(moments)#Para pruebas de normalidad de los residuos
library(MuMIn) #Selección AICC
library(openxlsx) #Permite la lectura de archivos tipo "xlsx"
library(plm) #Para datos de panel
library(readxl) # lectura archivo en excel
library(sandwich)#Para matriz de varianzas y covarianzas robusta
library(skedastic)#Para pruebas de heterocedasticidad
library(stargazer)#Para mejor apariencia de tablas
library(strucchange)#Para prueba estabilidad
library(tidyr)
library(truncreg)#Para regresión truncada 
```

# Introducción

::: {style="text-align: justify"}
Este código permite replicar las estimaciones del documento “Fijación ordinaria de comisiones máximas del sistema de tarjetas de pago”.

Como modelo principal, mediante una regresión lineal múltiple (MRL) se estiman los valores esperados para Costa Rica, dadas las características observables. En la selección del modelo se consideraron distintas combinaciones de variables de control, y se utilizó el criterio de información de Akaike con una corrección para tamaños de muestra pequeños (AICc) y el criterio de información bayesiana (BIC).

Posteriormente, se presenta un análisis de robustez a estimaciones alternativas, que considera los siguientes modelos:

-   Modelos de regresión truncados (MRT)
-   Modelo con transformación logarítmica de la variable dependiente (MTL)
-   Modelo con inclusión de efectos no lineales en la regresión (MENL)
-   Método bootstrap (MB)
-   Método de regresión robusta (MRR)
-   Modelo de regresión lineal omitiendo al Espacio Económico Europeo (MOEEE)
-   Modelo de regresión lineal restringiendo la muestra a países regulados (MPR)
-   Modelo de regresión lineal que amplía la definición de regulación (MADR)
-   Modelo de regresión lineal que utiliza las características observadas más recientes para cada país (MCR)
-   Modelo de regresión lineal con variables propuestas por el FMI (MFMI)
-   Modelo de regresión lineal sin variables relacionadas al uso de tarjetas (MSVR)
-   Modelo de regresión lineal con variables instrumentales (MVI)
-   Modelo de panel de datos con efectos aleatorios (MPDEA)
:::

# Datos

```{r, include=FALSE}
datos <- read_excel("D:/Base_datos_adquirencia_intercambio.xlsx")
```

::: {style="text-align: justify"}
El conjunto de datos contiene **`r nrow(datos)`** observaciones y **`r ncol(datos)`** variables.
:::

```{r, include=FALSE}
#Cantidad de países en la muestra (distintos de C.R.)

n_paises <- datos %>%
  filter(Codigo_Pais != "CRI") %>%
  distinct(Codigo_Pais) %>%
  nrow()

```

::: {style="text-align: justify"}
El total de países en la muestra sin considerar a Costa Rica es **`r n_paises`**.

Características de las variables:
:::

```{r, include=TRUE}
str(datos)
```

```{r, include=FALSE}
#Tablas de frecuencia

datos_tablas <- dplyr::select(datos, c("Codigo_Pais","OCDE","Regulacion_intercambio","Regulacion_adquirencia", "Dato_mas_reciente"))
datos_tablas <- datos_tablas %>% 
  filter(datos_tablas$Dato_mas_reciente == 1) %>% #Selecciona el dato más reciente
  filter(Codigo_Pais != "CRI") #Se elimina Costa Rica

#Cantidad de países miembros de la OCDE en la muestra
tabla_ocde <- as.data.frame(cbind( 
  c(0,1, "Total"), #Col_1
  
  c(sum(datos_tablas$OCDE==0), sum(datos_tablas$OCDE==1), nrow(datos_tablas)), #Col_2
  
  round(c(100*sum(datos_tablas$OCDE==0)/nrow(datos_tablas),
          100*sum(datos_tablas$OCDE==1)/nrow(datos_tablas),
          100*sum(datos_tablas$OCDE %in% c(0, 1))/nrow(datos_tablas)),2), #Col_3
  
  c(round(100*sum(datos_tablas$OCDE==0)/nrow(datos_tablas),2),
    round(100*sum(datos_tablas$OCDE %in% c(0, 1))/nrow(datos_tablas),2),""))) #Col_4

names(tabla_ocde) <- c("País pertenece a la OCDE", "Frecuencia", "Porcentaje", "Acumulado")

#Cantidad de países con regulación en tasa de intercambio
tabla_int <- as.data.frame(cbind( 
  c(0,1, "Total"), #Col_1
  
  c(sum(datos_tablas$Regulacion_intercambio==0),
  sum(datos_tablas$Regulacion_intercambio==1),
  nrow(datos_tablas)), #Col_2
 
   round(c(100*sum(datos_tablas$Regulacion_intercambio==0)/nrow(datos_tablas),                      
           100*sum(datos_tablas$Regulacion_intercambio==1)/nrow(datos_tablas),
           100*sum(datos_tablas$Regulacion_intercambio %in% c(0, 1))/nrow(datos_tablas)),2), #Col_3                         
         c(round(100*sum(datos_tablas$Regulacion_intercambio==0)/nrow(datos_tablas),2),
    round(100*sum(datos_tablas$Regulacion_intercambio %in% c(0, 1))/nrow(datos_tablas),2), ""))) #Col_4 

names(tabla_int) <- c("Regulación en intercambio", "Frecuencia", "Porcentaje", "Acumulado")

#Cantidad de países con regulación en tasa de adquirencia
tabla_adq <- as.data.frame(cbind( 
  c(0,1, "Total"), #Col_1
  
  c(sum(datos_tablas$Regulacion_adquirencia==0),
  sum(datos_tablas$Regulacion_adquirencia==1),
  nrow(datos_tablas)), #Col_2
 
   round(c(100*sum(datos_tablas$Regulacion_adquirencia==0)/nrow(datos_tablas),                      
           100*sum(datos_tablas$Regulacion_adquirencia==1)/nrow(datos_tablas),
           100*sum(datos_tablas$Regulacion_adquirencia %in% c(0, 1))/nrow(datos_tablas)),2), #Col_3                         
         c(round(100*sum(datos_tablas$Regulacion_adquirencia==0)/nrow(datos_tablas),2),
    round(100*sum(datos_tablas$Regulacion_adquirencia %in% c(0, 1))/nrow(datos_tablas),2), ""))) #Col_4 

names(tabla_adq) <- c("Regulación en adquirencia", "Frecuencia", "Porcentaje", "Acumulado")

```

::: {style="text-align: justify"}
Países miembros de la OCDE dentro de la muestra:
:::

```{r, include=TRUE}
tabla_ocde
```

::: {style="text-align: justify"}
Países dentro de la muestra con regulación de la comisión de intercambio:
:::

```{r, include=TRUE}
tabla_int
```

::: {style="text-align: justify"}
Países dentro de la muestra con regulación de la comisión de adquirencia:
:::

```{r, include=TRUE}
tabla_adq
```

```{r, include=FALSE}
rm(datos_tablas, tabla_ocde, tabla_int, tabla_adq, n_paises)
```

# Comisión de intercambio

## Estimación principal

::: {style="text-align: justify"}
Se estima una regresión lineal múltiple (MRL). A partir de los coeficientes estimados y las características observables para Costa Rica, se calcula la comisión de intercambio esperada para el país. Se calculan 5 modelos:

-   MRL_sc: estima una regresión lineal sin variables de control; únicamente se incluye la variable dicotómica que indica si el país cuenta con regulación.
-   MRL_aicc: corresponde a los resultados obtenidos con el criterio de información de Akaike (AICc) con una corrección para tamaños de muestra pequeñas.
-   MRL_bic: incluye las estimaciones al controlar por las variables seleccionadas mediante el criterio de información bayesiano (BIC).
-   MRL_aicc_et: incluyen efectos fijos de tiempo en las regresiones, para el modelo MRL_aicc.
-   MRL_bic_et: incluyen efectos fijos de tiempo en las regresiones, para el modelo MRL_bic.
:::

```{r, include=FALSE}
#Países con información sobre intercambio en la muestra (distintos de C.R.)

n_paises_int <- datos[!is.na(datos$Intercambio), ] %>%
  filter(Codigo_Pais != "CRI") %>%
  distinct(Codigo_Pais) %>%
  nrow()

```

::: {style="text-align: justify"}
El total de países con información sobre intercambio es **`r n_paises_int`**.
:::

### Definiciones

```{r, include=TRUE}
# Variables de control
controles <- c("Cajero_Autom", "Crimen", "Cuenta", "Credito", "Debito", "Densidad_Pob", "Densidad_Neg", "PIB_pc_PPA", "ICH", "Pago_digital", "Valor_prom", "Europa", "Asia_Este", "Am_lat")

# Efectos fijos de tiempo
ef_tiempo <- c("Año_2013", "Año_2014", "Año_2015", "Año_2016", "Año_2017", "Año_2018", "Año_2019", "Año_2020", "Año_2021", "Año_2022", "Año_2023", "Año_2024", "Año_2025")

# Fórmula que incluye las variables de control
formula_con_cont <- as.formula(
  paste("Intercambio ~ Regulacion_intercambio +", paste(controles, collapse = " + "))
)

# Fórmula que incluye las variables de control y los efectos fijos de tiempo
formula_con_cont_t <- as.formula(
  paste("Intercambio ~ Regulacion_intercambio +", paste(controles, collapse = " + ")," + ", paste(ef_tiempo, collapse = " + "))
)

# Variables fijas para los modelos con efectos de tiempo
variables_fijas <- c("Regulacion_intercambio", paste0("Año_", 2010:2025))

# Definir el conjuntos de datos para los modelos con variables de control
variables_constantes <- c("Pais", "Codigo_Pais", "Año", "Intercambio", "Regulacion_intercambio")
variables_interes <- c(variables_constantes, controles, ef_tiempo)

datos_int_con_cont <- datos %>%
  select(all_of(variables_interes)) %>%
  na.omit()

```

::: {style="text-align: justify"}
Variables de control: `r paste(capture.output(controles), collapse = "<br>")`.

Efectos fijos de tiempo: `r paste(capture.output(ef_tiempo), collapse = "<br>")`.
:::

### MRL_sc

::: {style="text-align: justify"}
Regresión lineal sin variables de control.
:::

#### Estimación

```{r, include=TRUE, results='asis'}
#Remover NA's de la columna "Intercambio" y dejar los demás.
datos_int_sin_cont <- datos %>% 
  filter(!is.na(Intercambio))

MRL_sc <- lm(Intercambio ~ Regulacion_intercambio, datos_int_sin_cont)

# Errores robustos
MRL_sc_corr_hc <- coeftest(MRL_sc, vcov. = vcovHC(MRL_sc, "HC1", adjust = TRUE, sandwich = TRUE))

MRL_sc_corr_hac <- coeftest(MRL_sc, vcov. = vcovHAC(MRL_sc, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MRL_sc, MRL_sc_corr_hc, type = "html", 
          title = "Estimación principal - Modelo sin variables de control", 
          column.labels = c("MCO","MCO HC"), model.names = F, 
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")
```

#### Pruebas estadísticas

```{r, include=TRUE}
#Prueba de normalidad
residuos_MRL_sc <- resid(MRL_sc) #Se guardan los residuos
jarque.test(residuos_MRL_sc) #Se computa la prueba con la función Nula: distribución normal residuos
car::qqPlot(MRL_sc)

#Test de White sobre el modelo original. Hipótesis nula: homocedasticidad
white(MRL_sc, interactions = TRUE)

```

#### Predicción

```{r, include=TRUE}
#Pronóstico para CR
data_CR <- datos
data_CR[ which(datos$Codigo_Pais=="CRI"), which(names(datos)=="Regulacion_intercambio") ]<-1
pred_int_MRL_sc<-predict(MRL_sc, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MRL_sc`**.
:::

```{r, include=FALSE}
#Crea el data frame para guardar la predicción
modelo <- c("Estimación principal - Modelo sin variables de control")
pronostico <- c(pred_int_MRL_sc)  
Pronosticos <- data.frame(Modelo = modelo, Pronostico = pronostico)
```

```{r, include=FALSE}
rm(datos_int_sin_cont, MRL_sc, MRL_sc_corr_hac, MRL_sc_corr_hc, residuos_MRL_sc, pred_int_MRL_sc, modelo, pronostico)
```

### MRL_aicc

::: {style="text-align: justify"}
Regresión lineal con variables de control, sin efectos fijos de tiempo. Utiliza el criterio de información de Akaike (AICc) con una corrección para tamaños de muestra pequeñas.
:::

#### Estimación

```{r, include=TRUE, results='asis'}
# Ajustar el modelo
modelo_completo <- lm(formula_con_cont, data = datos_int_con_cont)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRL_aicc <- dredge(modelo_completo, rank = "AICc", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MRL_aicc <- get.models(MRL_aicc, 1)[[1]]
MRL_aicc

# Errores robustos
MRL_aicc_corr_hc <- coeftest(MRL_aicc, vcov. = vcovHC(MRL_aicc, "HC1", adjust = TRUE, sandwich = TRUE))

MRL_aicc_corr_hac <- coeftest(MRL_aicc, vcov. = vcovHAC(MRL_aicc, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MRL_aicc, MRL_aicc_corr_hc,
          type = "html", 
          title = "Estimación principal - Modelo AICc sin efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")
```

#### Pruebas estadísticas

```{r, include=TRUE}
#Prueba de normalidad

residuos_MRL_aicc <- resid(MRL_aicc) #Se guardan los residuos
jarque.test(residuos_MRL_aicc) #Se computa la prueba con la función Nula: distribución normal residuos
car::qqPlot(MRL_aicc)

#Se obtiene el resultado del índice de condición de multicolinealidad
#create vector of VIF values
print("El VIF del modelo es:")
vif(MRL_aicc)

#Se ejecuta el test de White sobre el modelo original Nula: homocedasticidad
#white_lm(modelo_int_con_cont_aicc, interactions = TRUE)
white(MRL_aicc, interactions = TRUE)
plot(MRL_aicc$fit,MRL_aicc$res,xlab="Valores ajustados",ylab="Residuales")

#Prueba de especificación lineal RESET Nula: especificación correcta lineal
resettest(MRL_aicc, power = 2:3, type = c("fitted"), data = datos_int_con_cont, vcov = vcovHC) #Corregido por heterocedasticidad)

```

#### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MRL_aicc <- predict(MRL_aicc, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MRL_aicc`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación principal - Modelo AICc sin efectos fijos de tiempo")
pronostico <- c(pred_int_MRL_aicc)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRL_aicc, modelo, MRL_aicc_corr_hac, MRL_aicc_corr_hc, pred_int_MRL_aicc, pronostico, residuos_MRL_aicc)
```

### MRL_bic

::: {style="text-align: justify"}
Regresión lineal con variables de control, sin efectos fijos de tiempo. Utiliza el criterio de información bayesiana (BIC).
:::

#### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRL_bic <- dredge(modelo_completo, rank = "BIC", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MRL_bic <- get.models(MRL_bic, 1)[[1]]
MRL_bic

# Errores robustos
MRL_bic_corr_hc<-coeftest(MRL_bic, vcov. = vcovHC(MRL_bic, "HC1", adjust = TRUE, sandwich = TRUE))

MRL_bic_corr_hac<-coeftest(MRL_bic, vcov. = vcovHAC(MRL_bic, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MRL_bic, MRL_bic_corr_hc, type = "html",
          title = "Estimación principal - Modelo BIC sin efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F, 
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")


```

#### Pruebas estadísticas

```{r, include=TRUE}
#Prueba de normalidad
#Se guardan los residuos
residuos_MRL_bic<-resid(MRL_bic)
#Se computa la prueba con la función Nula: distribución normal residuos
jarque.test(residuos_MRL_bic)

#Se obtiene el resultado del índice de condición de multicolinealidad
#create vector of VIF values
print("El VIF del modelo es:")
vif(MRL_bic)

#Se ejecuta el test de White sobre el modelo original Nula: homocedasticidad
#white_lm(modelo_int_con_cont_bic, interactions = TRUE)
white(MRL_bic, interactions = TRUE)
#Prueba de especificación lineal RESET Nula: especificación correcta lineal
resettest(MRL_bic, power = 2:3, type = c("fitted"), data = datos_int_con_cont, vcov = vcovHC)#Corregido por heterocedasticidad)

```

#### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MRL_bic <- predict(MRL_bic, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MRL_bic`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación principal - Modelo BIC sin efectos fijos de tiempo")
pronostico <- c(pred_int_MRL_bic)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRL_bic, modelo, MRL_bic_corr_hac, MRL_bic_corr_hc, pred_int_MRL_bic, pronostico, residuos_MRL_bic)
```

### MRL_aicc_et

::: {style="text-align: justify"}
Regresión lineal con variables de control, con efectos fijos de tiempo. Utiliza el criterio de información de Akaike (AICc) con una corrección para tamaños de muestra pequeñas.
:::

#### Estimación

```{r, include=TRUE, results='asis'}
# Ajustar el modelo
modelo_completo_et <- lm(formula_con_cont_t, data = datos_int_con_cont)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRL_aicc_et <- dredge(modelo_completo_et, rank = "AICc", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MRL_aicc_et <- get.models(MRL_aicc_et, 1)[[1]]
MRL_aicc_et

# Errores robustos
MRL_aicc_et_corr_hc <- coeftest(MRL_aicc_et, vcov. = vcovHC(MRL_aicc_et, "HC1", adjust = TRUE, sandwich = TRUE))

MRL_aicc_et_corr_hac <- coeftest(MRL_aicc_et, vcov. = vcovHAC(MRL_aicc_et, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MRL_aicc_et, MRL_aicc_et_corr_hc,
          type = "html", 
          title = "Estimación principal - Modelo AICc con efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")

MRL_aicc_et_corr_hc
MRL_aicc_et_corr_hac

```

#### Pruebas estadísticas

```{r, include=TRUE}

#Prueba de normalidad
#Se guardan los residuos
residuos_MRL_aicc_et<-resid(MRL_aicc_et)
#Se computa la prueba con la función Nula: distribución normal residuos
jarque.test(residuos_MRL_aicc_et)
car::qqPlot(residuos_MRL_aicc_et)

#Se ejecuta el test de White sobre el modelo original Nula: homocedasticidad
#white_lm(modelo_int_con_cont_ef_aicc, interactions = TRUE)
white(MRL_aicc_et, interactions = TRUE)

```

#### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MRL_aicc_et<-predict(MRL_aicc_et, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MRL_aicc_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación principal - Modelo AICc con efectos fijos de tiempo")
pronostico <- c(pred_int_MRL_aicc_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRL_aicc_et, modelo, MRL_aicc_et_corr_hac, MRL_aicc_et_corr_hc, pred_int_MRL_aicc_et, pronostico, residuos_MRL_aicc_et)
```

### MRL_bic_et

::: {style="text-align: justify"}
Regresión lineal con variables de control, con efectos fijos de tiempo. Utiliza el criterio de información bayesiana (BIC).
:::

#### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRL_bic_et <- dredge(modelo_completo_et, rank = "BIC", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MRL_bic_et <- get.models(MRL_bic_et, 1)[[1]]
MRL_bic_et

# Errores robustos
MRL_bic_et_corr_hc <- coeftest(MRL_bic_et, vcov. = vcovHC(MRL_bic_et, "HC1", adjust = TRUE, sandwich = TRUE))

MRL_bic_et_corr_hac <- coeftest(MRL_bic_et, vcov. = vcovHAC(MRL_bic_et, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MRL_bic_et, MRL_bic_et_corr_hc, type = "html",
          title = "Estimación principal - Modelo BIC con efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F, 
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")

```

#### Pruebas estadísticas

```{r, include=TRUE}
#Prueba de normalidad
#Se guardan los residuos
residuos_MRL_bic_et <- resid(MRL_bic_et)
#Se computa la prueba con la función Nula: distribución normal residuos
jarque.test(residuos_MRL_bic_et)
white(MRL_bic_et, interactions = TRUE)
```

#### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MRL_bic_et <- predict(MRL_bic_et, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MRL_bic_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación principal - Modelo BIC con efectos fijos de tiempo")
pronostico <- c(pred_int_MRL_bic_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRL_bic_et, modelo, MRL_bic_et_corr_hac, MRL_bic_et_corr_hc, pred_int_MRL_bic_et, pronostico, residuos_MRL_bic_et)
```

## Pruebas de robustez

::: {style="text-align: justify"}
A continuación, se evalúa cómo cambian los resultados anteriores al estimar las comisiones con modelos alternativos.
:::

### Modelos de regresión truncados

::: {style="text-align: justify"}
Mientras que el modelo de regresión lineal no impone ninguna restricción a las comisiones, es de esperar que las comisiones sean necesariamente mayores a cero como mínimo, o mayores a los costos asociados en particular. Para tratar estos casos, se pueden utilizar modelos de variable dependiente truncada.

Se estiman estos modelos bajo dos supuestos sobre el punto de truncamiento. Un primer punto de partida es que el punto de truncamiento es en cero. El segundo supuesto es que el punto de truncamiento sucede en el mínimo de la comisión de intercambio observada en la muestra.

Además, se estima cada modelo con las variables seleccionadas por los criterios de información de AICc y BIC.
:::

#### MRT_aicc_0

##### Estimación

```{r, include=TRUE}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRT_aicc_0 <- dredge(modelo_completo, rank = "AICc", fixed = "Regulación_intercambio")

# Extraer el mejor modelo según AICc
MRT_aicc_0 <- names(coef(get.models(MRT_aicc_0, 1)[[1]]))
MRT_aicc_0 <- MRT_aicc_0[MRT_aicc_0 != "(Intercept)"]  # Excluir intercepto

# Ajustar modelo truncado con límite inferior 0
formula_trunc <- as.formula(paste("Intercambio ~", paste(MRT_aicc_0, collapse = " + ")))
MRT_aicc_0 <- truncreg(formula_trunc, data = datos_int_con_cont, point = 0, direction = "left")

summary(MRT_aicc_0)
```

##### Predicción

```{r, include=TRUE, warning=FALSE}
#Pronóstico para CR
pred_int_MRT_aicc_0 <- predict(MRT_aicc_0, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MRT_aicc_0`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación truncada en 0 - Modelo AICc sin efectos fijos de tiempo")
pronostico <- c(pred_int_MRT_aicc_0)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRT_aicc_0, modelo, pred_int_MRT_aicc_0, pronostico, lastFuncGrad, lastFuncParam)
```

#### MRT_bic_0

##### Estimación

```{r, include=TRUE}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRT_bic_0 <- dredge(modelo_completo, rank = "BIC", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MRT_bic_0 <- names(coef(get.models(MRT_bic_0, 1)[[1]]))
MRT_bic_0 <- MRT_bic_0[MRT_bic_0 !="(Intercept)"]  # Excluir intercepto

# Ajustar modelo truncado con límite inferior 0
formula_trunc <- as.formula(paste("Intercambio ~", paste(MRT_bic_0, collapse = " + ")))
MRT_bic_0 <- truncreg(formula_trunc, data = datos_int_con_cont, point = 0, direction = "left")
summary(MRT_bic_0)
```

##### Predicción

```{r truncado 0 bic int , include=TRUE, warning=FALSE}
#Pronóstico para CR
pred_int_MRT_bic_0 <- predict(MRT_bic_0, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MRT_bic_0`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación truncada en 0 - Modelo BIC sin efectos fijos de tiempo")
pronostico <- c(pred_int_MRT_bic_0)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRT_bic_0, modelo, pred_int_MRT_bic_0, pronostico, lastFuncGrad, lastFuncParam)
```

#### MRT_aicc_et_0

##### Estimación

```{r, include=TRUE}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRT_aicc_et_0 <- dredge(modelo_completo_et, rank = "AICc", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MRT_aicc_et_0 <- get.models(MRT_aicc_et_0, 1)[[1]]

# Obtener coeficientes
coeficientes <- coef(MRT_aicc_et_0)

# Filtrar nombres de coeficientes que no sean NA
coeficientes_validos <- names(coeficientes[!is.na(coeficientes)])

# Excluir el intercepto
variables_truncadas <- coeficientes_validos[coeficientes_validos != "(Intercept)"]

# Construir fórmula para modelo truncado
formula_trunc <- as.formula(paste("Intercambio ~", paste(variables_truncadas, collapse = " + ")))

# Ajustar modelo truncado con límite inferior 0
MRT_aicc_et_0 <- truncreg(formula_trunc, data = datos_int_con_cont, point = 0, direction = "left")

summary(MRT_aicc_et_0)

```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MRT_aicc_et_0<-predict(MRT_aicc_et_0, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MRT_aicc_et_0`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación truncada en 0 - Modelo AICc con efectos fijos de tiempo")
pronostico <- c(pred_int_MRT_aicc_et_0)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRT_aicc_et_0, modelo, pred_int_MRT_aicc_et_0, pronostico, lastFuncGrad, lastFuncParam, coeficientes, coeficientes_validos)
```

#### MRT_bic_et_0

##### Estimación

```{r, include=TRUE}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRT_bic_et_0 <- dredge(modelo_completo_et, rank = "BIC", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MRT_bic_et_0 <- get.models(MRT_bic_et_0, 1)[[1]]

# Obtener coeficientes
coeficientes <- coef(MRT_bic_et_0)

# Filtrar nombres de coeficientes que no sean NA
coeficientes_validos <- names(coeficientes[!is.na(coeficientes)])

# Excluir el intercepto
variables_truncadas <- coeficientes_validos[coeficientes_validos != "(Intercept)"]

# Construir fórmula para modelo truncado
formula_trunc <- as.formula(paste("Intercambio ~", paste(variables_truncadas, collapse = " + ")))

# Ajustar modelo truncado con límite inferior 0
MRT_bic_et_0 <- truncreg(formula_trunc, data = datos_int_con_cont, point = 0, direction = "left")

summary(MRT_bic_et_0)
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MRT_bic_et_0 <- predict(MRT_bic_et_0, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MRT_bic_et_0`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación truncada en 0 - Modelo BIC con efectos fijos de tiempo")
pronostico <- c(pred_int_MRT_bic_et_0)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRT_bic_et_0, modelo, pred_int_MRT_bic_et_0, pronostico, lastFuncGrad, lastFuncParam, coeficientes, coeficientes_validos)
```

#### MRT_aicc_026

##### Estimación

```{r, include=TRUE}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRT_aicc_026 <- dredge(modelo_completo, rank = "AICc", fixed = "Regulación_intercambio")

# Extraer el mejor modelo según AICc
MRT_aicc_026 <- names(coef(get.models(MRT_aicc_026, 1)[[1]]))
MRT_aicc_026 <- MRT_aicc_026[MRT_aicc_026 != "(Intercept)"]  # Excluir intercepto

# Ajustar modelo truncado con límite inferior 0.26
formula_trunc <- as.formula(paste("Intercambio ~", paste(MRT_aicc_026, collapse = " + ")))
MRT_aicc_026 <- truncreg(formula_trunc, data = datos_int_con_cont, point = 0.26, direction = "left")

summary(MRT_aicc_026)
```

##### Predicción

```{r, include=TRUE, warning=FALSE}
#Pronóstico para CR
pred_int_MRT_aicc_026 <- predict(MRT_aicc_026, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MRT_aicc_026`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación truncada en 0.26 - Modelo AICc sin efectos fijos de tiempo")
pronostico <- c(pred_int_MRT_aicc_026)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRT_aicc_026, modelo, pred_int_MRT_aicc_026, pronostico, lastFuncGrad, lastFuncParam)
```

#### MRT_bic_026

##### Estimación

```{r, include=TRUE}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRT_bic_026 <- dredge(modelo_completo, rank = "BIC", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MRT_bic_026 <- names(coef(get.models(MRT_bic_026, 1)[[1]]))
MRT_bic_026 <- MRT_bic_026[MRT_bic_026 !="(Intercept)"]  # Excluir intercepto

# Ajustar modelo truncado con límite inferior 0.26
formula_trunc <- as.formula(paste("Intercambio ~", paste(MRT_bic_026, collapse = " + ")))
MRT_bic_026 <- truncreg(formula_trunc, data = datos_int_con_cont, point = 0.26, direction = "left")

summary(MRT_bic_026)
```

##### Predicción

```{r, include=TRUE, warning=FALSE}
#Pronóstico para CR
pred_int_MRT_bic_026 <- predict(MRT_bic_026, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MRT_bic_026`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación truncada en 0.26 - Modelo BIC sin efectos fijos de tiempo")
pronostico <- c(pred_int_MRT_bic_026)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRT_bic_026, modelo, pred_int_MRT_bic_026, pronostico, lastFuncGrad, lastFuncParam)
```

#### MRT_aicc_et_026

##### Estimación

```{r, include=TRUE}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRT_aicc_et_026 <- dredge(modelo_completo_et, rank = "AICc", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MRT_aicc_et_026 <- get.models(MRT_aicc_et_026, 1)[[1]]

# Obtener coeficientes
coeficientes <- coef(MRT_aicc_et_026)

# Filtrar nombres de coeficientes que no sean NA
coeficientes_validos <- names(coeficientes[!is.na(coeficientes)])

# Excluir el intercepto
variables_truncadas <- coeficientes_validos[coeficientes_validos != "(Intercept)"]

# Construir fórmula para modelo truncado
formula_trunc <- as.formula(paste("Intercambio ~", paste(variables_truncadas, collapse = " + ")))

# Ajustar modelo truncado con límite inferior 0.26
MRT_aicc_et_026 <- truncreg(formula_trunc, data = datos_int_con_cont, point = 0.26, direction = "left")

summary(MRT_aicc_et_026)
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MRT_aicc_et_026 <- predict(MRT_aicc_et_026, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MRT_aicc_et_026`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación truncada en 0.26 - Modelo AICc con efectos fijos de tiempo")
pronostico <- c(pred_int_MRT_aicc_et_026)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRT_aicc_et_026, modelo, pred_int_MRT_aicc_et_026, pronostico, lastFuncGrad, lastFuncParam, coeficientes, coeficientes_validos)
```

#### MRT_bic_et_026

##### Estimación

```{r, include=TRUE}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRT_bic_et_026 <- dredge(modelo_completo_et, rank = "BIC", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MRT_bic_et_026 <- get.models(MRT_bic_et_026, 1)[[1]]

# Obtener coeficientes
coeficientes <- coef(MRT_bic_et_026)

# Filtrar nombres de coeficientes que no sean NA
coeficientes_validos <- names(coeficientes[!is.na(coeficientes)])

# Excluir el intercepto
variables_truncadas <- coeficientes_validos[coeficientes_validos != "(Intercept)"]

# Construir fórmula para modelo truncado
formula_trunc <- as.formula(paste("Intercambio ~", paste(variables_truncadas, collapse = " + ")))

# Ajustar modelo truncado con límite inferior 0.26
MRT_bic_et_026 <- truncreg(formula_trunc, data = datos_int_con_cont, point = 0.26, direction = "left")

summary(MRT_bic_et_026)
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MRT_bic_et_026 <- predict(MRT_bic_et_026, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MRT_bic_et_026`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación truncada en 0.26 - Modelo BIC con efectos fijos de tiempo")
pronostico <- c(pred_int_MRT_bic_et_026)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRT_bic_et_026, modelo, pred_int_MRT_bic_et_026, pronostico, lastFuncGrad, lastFuncParam, coeficientes, coeficientes_validos, variables_truncadas, formula_trunc)
```

### Modelos de regresión con transformación logarítmica de la variable dependiente

::: {style="text-align: justify"}
Al no existir ningún tipo de restricción en la estimación de los coeficientes, cabe la posibilidad de que un modelo de regresión lineal genere predicciones negativas. Sin embargo, ninguna de las predicciones obtenidas para el intercambio o adquirencia son negativas.

A pesar de no significar un potencial problema, una manera de garantizar que las predicciones sean positivas es recurrir a una transformación logarítmica de la variable dependiente.
:::

```{r, include=FALSE}
## Transformación a logaritmos de las comisiones
datos_int_con_cont_log <- datos_int_con_cont
datos_int_con_cont_log$Intercambio <- log(datos_int_con_cont_log$Intercambio)
```

#### MTL_aicc

##### Estimación

```{r, include=TRUE}
# Ajustar el modelo
modelo_completo_log <- lm(formula_con_cont, data = datos_int_con_cont_log)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MTL_aicc <- dredge(modelo_completo_log, rank = "AICc", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MTL_aicc <- get.models(MTL_aicc, 1)[[1]]
MTL_aicc

# Errores robustos
MTL_aicc_corr_hc <- coeftest(MTL_aicc, vcov. = vcovHC(MTL_aicc, "HC1", adjust = TRUE, sandwich = TRUE))

MTL_aicc_corr_hac <- coeftest(MTL_aicc, vcov. = vcovHAC(MTL_aicc, adjust = TRUE, sandwich = TRUE))

summary(MTL_aicc)
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MTL_aicc <- predict(MTL_aicc, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MTL_aicc`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación con transformación logarítmica - Modelo AICc sin efectos fijos de tiempo")
pronostico <- c(pred_int_MTL_aicc)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MTL_aicc, modelo, pred_int_MTL_aicc, pronostico, MTL_aicc_corr_hc, MTL_aicc_corr_hac)
```

#### MTL_bic

##### Estimación

```{r, include=TRUE}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MTL_bic <- dredge(modelo_completo_log, rank = "BIC", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MTL_bic <- get.models(MTL_bic, 1)[[1]]
MTL_bic

# Errores robustos
MTL_bic_corr_hc <- coeftest(MTL_bic, vcov. = vcovHC(MTL_bic, "HC1", adjust = TRUE, sandwich = TRUE))

MTL_bic_corr_hac <- coeftest(MTL_bic, vcov. = vcovHAC(MTL_bic, adjust = TRUE, sandwich = TRUE))

summary(MTL_bic)
```

##### Predicción

```{r, include=TRUE}

#Pronóstico para CR
pred_int_MTL_bic <- predict(MTL_bic, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MTL_bic`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación con transformación logarítmica - Modelo BIC sin efectos fijos de tiempo")
pronostico <- c(pred_int_MTL_bic)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MTL_bic, modelo, pred_int_MTL_bic, pronostico, MTL_bic_corr_hc, MTL_bic_corr_hac, modelo_completo_log)
```

#### MTL_aicc_et

##### Estimación

```{r, include=TRUE}
# Ajustar el modelo
modelo_completo_log_et <- lm(formula_con_cont_t, data = datos_int_con_cont_log)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MTL_aicc_et <- dredge(modelo_completo_log_et, rank = "AICc", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MTL_aicc_et <- get.models(MTL_aicc_et, 1)[[1]]
MTL_aicc_et

# Errores robustos
MTL_aicc_et_corr_hc <- coeftest(MTL_aicc_et, vcov. = vcovHC(MTL_aicc_et, "HC1", adjust = TRUE, sandwich = TRUE))

MTL_aicc_et_corr_hac <- coeftest(MTL_aicc_et, vcov. = vcovHAC(MTL_aicc_et, adjust = TRUE, sandwich = TRUE))

summary(MTL_aicc_et)
```

##### Predicción

```{r, include= TRUE}
#Pronóstico para CR
pred_int_MTL_aicc_et <- predict(MTL_aicc_et, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MTL_aicc_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación con transformación logarítmica - Modelo AICc con efectos fijos de tiempo")
pronostico <- c(pred_int_MTL_aicc_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MTL_aicc_et, modelo, pred_int_MTL_aicc_et, pronostico, MTL_aicc_et_corr_hc, MTL_aicc_et_corr_hac)
```

#### MTL_bic_et

##### Estimación

```{r, include=TRUE}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MTL_bic_et <- dredge(modelo_completo_log_et, rank = "BIC", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MTL_bic_et <- get.models(MTL_bic_et, 1)[[1]]
MTL_bic_et

# Errores robustos
MTL_bic_et_corr_hc <- coeftest(MTL_bic_et, vcov. = vcovHC(MTL_bic_et, "HC1", adjust = TRUE, sandwich = TRUE))

MTL_bic_et_corr_hac <- coeftest(MTL_bic_et, vcov. = vcovHAC(MTL_bic_et, adjust = TRUE, sandwich = TRUE))

summary(MTL_bic_et)
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MTL_bic_et <- predict(MTL_bic_et, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **\`r pred_int_MTL_bic_et**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación con transformación logarítmica - Modelo BIC con efectos fijos de tiempo")
pronostico <- c(pred_int_MTL_bic_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MTL_bic_et, modelo, pred_int_MTL_bic_et, pronostico, MTL_bic_et_corr_hc, MTL_bic_et_corr_hac, datos_int_con_cont_log, modelo_completo_log_et)
```

### Modelo de regresión con efectos no lineales

::: {style="text-align: justify"}
Para permitir la presencia de efectos no lineales en la regresión, se estima una especificación donde se incluyen como potenciales candidatos en la selección del modelo los términos cuadráticos de las variables de control.
:::

```{r, include=FALSE}
## Arreglo de los datos
data_CR_enl <- data_CR
datos_int_con_cont_enl <- datos_int_con_cont

controles_numericos <- setdiff(controles, c("Europa", "Asia_Este", "Am_lat"))

for (var in controles_numericos) {
  datos_int_con_cont_enl[[paste0(var, "_sqrt")]] <- datos_int_con_cont_enl[[var]]^2
}

for (var in controles_numericos) {
  data_CR_enl[[paste0(var, "_sqrt")]] <- data_CR_enl[[var]]^2
}

controles_sqrt <- c(controles, paste0(controles_numericos, "_sqrt"))


# Fórmula para incluir las variables de control
formula_con_cont_enl <- as.formula(
  paste("Intercambio ~ Regulacion_intercambio +", paste(controles_sqrt, collapse = " + "))
)

# Fórmula para incluir las variables de control con efectos fijos de tiempo
formula_con_cont_t_enl <- as.formula(
  paste("Intercambio ~ Regulacion_intercambio +", paste(controles_sqrt, collapse = " + ")," + ", paste(ef_tiempo, collapse = " + "))
)

```

#### MENL_aicc

##### Estimación

::: {style="text-align: justify"}
Se proporciona el código completo para su ejecución en R. No obstante, debido al tiempo excesivo que toma el procesamiento automático, se han incorporado manualmente las variables que optimizan el modelo. Estas variables fueron identificadas previamente mediante un proceso de selección realizado en Stata, con el fin de reducir significativamente el tiempo de ejecución en R.
:::

```{r, include=TRUE}
# Ajustar el modelo (automático R)
#modelo_completo_enl <- lm(formula_con_cont_enl, data = datos_int_con_cont_enl)

# Ajustar el modelo (manual, selección previa de variables en Stata)
modelo_completo_enl <- lm(Intercambio ~ Regulacion_intercambio + 
                                     Crimen + Cuenta + Credito + Debito + 
                                     Densidad_Pob + ICH + Pago_digital + 
                                     Europa + Cajero_Autom_sqrt + Crimen_sqrt +      
                                     Cuenta_sqrt + Credito_sqrt + ICH_sqrt +
                                     Pago_digital_sqrt, data = datos_int_con_cont_enl) 

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MENL_aicc <- dredge(modelo_completo_enl, rank = "AICc", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MENL_aicc <- get.models(MENL_aicc, 1)[[1]]
MENL_aicc

# Errores robustos
MENL_aicc_corr_hc <- coeftest(MENL_aicc, vcov. = vcovHC(MENL_aicc, "HC1", adjust = TRUE, sandwich = TRUE))

MENL_aicc_corr_hac <- coeftest(MENL_aicc, vcov. = vcovHAC(MENL_aicc, adjust = TRUE, sandwich = TRUE))

summary(MENL_aicc)
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MENL_aicc <- predict(MENL_aicc, newdata = data_CR_enl, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MENL_aicc`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación con efectos no lineales - Modelo AICc sin efectos fijos de tiempo")
pronostico <- c(pred_int_MENL_aicc)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MENL_aicc, modelo, pred_int_MENL_aicc, pronostico, MENL_aicc_corr_hc, MENL_aicc_corr_hac, var, modelo_completo_enl)
```

#### MENL_bic

##### Estimación

::: {style="text-align: justify"}
Se proporciona el código completo para su ejecución en R. No obstante, debido al tiempo excesivo que toma el procesamiento automático, se han incorporado manualmente las variables que optimizan el modelo. Estas variables fueron identificadas previamente mediante un proceso de selección realizado en Stata, con el fin de reducir significativamente el tiempo de ejecución en R.
:::

```{r, include=TRUE}
# Ajustar el modelo (automático R)
#modelo_completo_enl <- lm(formula_con_cont_enl, data = datos_int_con_cont_enl)

# Ajustar el modelo (manual, selección previa de variables en Stata)
modelo_completo_enl <- lm(Intercambio ~ Regulacion_intercambio +
                                     Cuenta + Credito + Densidad_Pob + 
                                     Europa + Credito_sqrt, data = datos_int_con_cont_enl) 

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MENL_bic <- dredge(modelo_completo_enl, rank = "BIC", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MENL_bic <- get.models(MENL_bic, 1)[[1]]
MENL_bic

# Errores robustos
MENL_bic_corr_hc <- coeftest(MENL_bic, vcov. = vcovHC(MENL_bic, "HC1", adjust = TRUE, sandwich = TRUE))

MENL_bic_corr_hac <- coeftest(MENL_bic, vcov. = vcovHAC(MENL_bic, adjust = TRUE, sandwich = TRUE))

summary(MENL_bic)
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MENL_bic <- predict(MENL_bic, newdata = data_CR_enl, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MENL_bic`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación con efectos no lineales - Modelo BIC sin efectos fijos de tiempo")
pronostico <- c(pred_int_MENL_bic)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MENL_bic, modelo, pred_int_MENL_bic, pronostico, MENL_bic_corr_hc, MENL_bic_corr_hac, var, modelo_completo_enl, formula_con_cont_enl)
```

#### MENL_aicc_et \*\*

##### Estimación

::: {style="text-align: justify"}
Se proporciona el código completo para su ejecución en R. No obstante, debido al tiempo excesivo que toma el procesamiento automático, se han incorporado manualmente las variables que optimizan el modelo. Estas variables fueron identificadas previamente mediante un proceso de selección realizado en Stata, con el fin de reducir significativamente el tiempo de ejecución en R.
:::

```{r, include=TRUE}
# Ajustar el modelo (automático R)
#modelo_completo_et_enl <- lm(formula_con_cont_t_enl, data = datos_int_con_cont_enl)

# Ajustar el modelo (manual, selección previa de variables en Stata)
modelo_completo_et_enl <- lm(Intercambio ~ Regulacion_intercambio + 
                                        Crimen + Cuenta + Credito + Debito + 
                                        Densidad_Pob + PIB_pc_PPA + ICH + Pago_digital +                                          Europa + Cajero_Autom_sqrt + Crimen_sqrt + 
                                        Cuenta_sqrt + Credito_sqrt + PIB_pc_PPA_sqrt +
                                        ICH_sqrt + Pago_digital_sqrt + Año_2013 + Año_2014                                         + Año_2015 + Año_2016 + Año_2017 + Año_2018 +
                                        Año_2019 + Año_2020 + Año_2021 + Año_2022 +   
                                        Año_2023 + Año_2024, 
                                        data = datos_int_con_cont_enl) 

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MENL_aicc_et <- dredge(modelo_completo_et_enl, rank = "AICc", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MENL_aicc_et <- get.models(MENL_aicc_et, 1)[[1]]
MENL_aicc_et

# Errores robustos
MENL_aicc_et_corr_hc <- coeftest(MENL_aicc_et, vcov. = vcovHC(MENL_aicc_et, "HC1", adjust = TRUE, sandwich = TRUE))

MENL_aicc_et_corr_hac <- coeftest(MENL_aicc_et, vcov. = vcovHAC(MENL_aicc_et, adjust = TRUE, sandwich = TRUE))

summary(MENL_aicc_et)
```

##### Predicción

```{r}
#Pronóstico para CR
pred_int_MENL_aicc_et <- predict(MENL_aicc_et, newdata = data_CR_enl, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MENL_aicc_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación con efectos no lineales - Modelo AICc con efectos fijos de tiempo")
pronostico <- c(pred_int_MENL_aicc_et)
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MENL_aicc_et, modelo, pred_int_MENL_aicc_et, pronostico, MENL_aicc_et_corr_hc, MENL_aicc_et_corr_hac, modelo_completo_et_enl)
```

#### MENL_bic_et

##### Estimación

::: {style="text-align: justify"}
Se proporciona el código completo para su ejecución en R. No obstante, debido al tiempo excesivo que toma el procesamiento automático, se han incorporado manualmente las variables que optimizan el modelo. Estas variables fueron identificadas previamente mediante un proceso de selección realizado en Stata, con el fin de reducir significativamente el tiempo de ejecución en R.
:::

```{r, include=TRUE}
# Ajustar el modelo (automático R)
#modelo_completo_et_enl <- lm(formula_con_cont_t_enl, data = datos_int_con_cont_enl)

# Ajustar el modelo (manual, selección previa de variables en Stata)
modelo_completo_et_enl <- lm(Intercambio ~ Regulacion_intercambio + 
                                       Cuenta + Credito + Densidad_Pob + Europa + 
                                       Credito_sqrt + Año_2013 + Año_2014 + Año_2015 + 
                                       Año_2016 + Año_2017 + Año_2018 + Año_2019 + 
                                       Año_2020 + Año_2021 + Año_2022 + Año_2023 +  
                                       Año_2024, data = datos_int_con_cont_enl) 

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MENL_bic_et <- dredge(modelo_completo_et_enl, rank = "BIC", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MENL_bic_et <- get.models(MENL_bic_et, 1)[[1]]
MENL_bic_et

# Errores robustos
MENL_bic_et_corr_hc <- coeftest(MENL_bic_et, vcov. = vcovHC(MENL_bic_et, "HC1", adjust = TRUE, sandwich = TRUE))

MENL_bic_et_corr_hac <- coeftest(MENL_bic_et, vcov. = vcovHAC(MENL_bic_et, adjust = TRUE, sandwich = TRUE))

summary(MENL_bic_et)
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MENL_bic_et <- predict(MENL_bic_et, newdata = data_CR_enl, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MENL_bic_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación con efectos no lineales - Modelo BIC con efectos fijos de tiempo")
pronostico <- c(pred_int_MENL_bic_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MENL_bic_et, modelo, pred_int_MENL_bic_et, pronostico, MENL_bic_et_corr_hc, MENL_bic_et_corr_hac, datos_int_con_cont_enl, data_CR_enl, modelo_completo_et_enl, controles_numericos, controles_sqrt, formula_con_cont_t_enl)
```

### Método Bootstrap

::: {style="text-align: justify"}
Una manera alternativa de obtener los intervalos de confianza de las predicciones de las comisiones es mediante el método bootstrap. Este método de remuestreo permite obtener una estimación no paramétrica de la distribución empírica de las comisiones de intercambio y adquirencia.

Para desarrollar este ejercicio, se realizan 1.000 remuestreos, y para cada muestra obtenida se sigue el mismo tipo de análisis presentado en la estimación principal. Esto da como resultado una comisión de adquirencia o intercambio promedio, según sea el caso, para cada iteración.

A continuación, se presenta el código completo en R. No obstante, debido al tiempo estimado de ejecución (aproximadamente 16 horas), los chunks se encuentran deshabilitados.
:::

```{r, include= TRUE, eval=FALSE}
start_time <- Sys.time()
```

```{r, include= TRUE, eval=FALSE}
# Función segura para convertir factores con un solo nivel
convert_factors <- function(df) {
  df %>%
    mutate(across(where(is.factor), ~ if (nlevels(.) < 2) as.character(.) else .))
}

# Configurar plan de ejecución paralela (puedes ajustar el número de workers)


bootstrap_iter <- function(i) {
  cat("Iteración:", i, "\n")
  
  resultado <- list(
    Intercambio_AICC = NA_real_,
    Intercambio_BIC = NA_real_,
    Intercambio_AICC_EF = NA_real_,
    Intercambio_BIC_EF = NA_real_
  )
  
  datos <- readxl::read_excel("D:/Base_datos_adquirencia_intercambio.xlsx")
  
  datos <- datos %>%
    mutate(across(all_of(c("Regulacion_intercambio", "Regulacion_adquirencia",
                           "Europa", "Asia_Este", "Am_lat", "OCDE", "EEE",  
                           "Dato_mas_reciente",
                           paste0("Año_", 2010:2025))), as.factor))
  
  datos_int_con_cont <- datos %>%
    dplyr::select(-Adquirencia) %>%
    na.omit()
  
  data_CR <- datos[which(datos$Codigo_Pais == "CRI"), which(names(datos) != "Adquirencia")]
  
  controles <- c("Cajero_Autom", "Crimen", "Cuenta", "Credito", "Debito", "Densidad_Pob", "Densidad_Neg", "PIB_pc_PPA", "ICH", "Pago_digital", "Valor_prom", "Europa", "Asia_Este", "Am_lat")
  ef_tiempo <- paste0("Año_", 2013:2025)
  variables_fijas <- c("Regulacion_intercambio", paste0("Año_", 2010:2025))
  
  muestra <- tryCatch({
    temp <- datos_int_con_cont[sample(nrow(datos_int_con_cont), replace = TRUE), ]
    temp <- convert_factors(temp)
    temp <- droplevels(temp)
    
    # Detectar factores con un solo nivel
    factores_un_nivel <- names(Filter(function(x) is.factor(x) && nlevels(x) < 2, temp))
    
    # Ajustar variables de control y efectos fijos
    controles_validos <- setdiff(controles, factores_un_nivel)
    ef_tiempo_validos <- setdiff(ef_tiempo, factores_un_nivel)
    variables_fijas_validas <- setdiff(variables_fijas, factores_un_nivel)
    
    # Crear fórmulas dinámicas
    formula_con_cont <- as.formula(
      paste("Intercambio ~ Regulacion_intercambio +", paste(controles_validos, collapse = " + "))
    )
    
    formula_con_cont_t <- as.formula(
      paste("Intercambio ~ Regulacion_intercambio +", paste(c(controles_validos, ef_tiempo_validos), collapse = " + "))
    )
    
    list(
      muestra = temp,
      formula_con_cont = formula_con_cont,
      formula_con_cont_t = formula_con_cont_t,
      variables_fijas_validas = variables_fijas_validas
    )
    
  }, error = function(e) {
    cat("Error al generar muestra en iteración", i, ":", conditionMessage(e), "\n")
    return(NULL)
  })
  
  if (is.null(muestra)) return(resultado)
  
  temp <- muestra$muestra
  formula_con_cont <- muestra$formula_con_cont
  formula_con_cont_t <- muestra$formula_con_cont_t
  variables_fijas_validas <- muestra$variables_fijas_validas
  
  # AICC sin EF
  tryCatch({
    modelo <- lm(formula_con_cont, data = temp)
    options(na.action = "na.fail")
    modelos <- dredge(modelo, fixed = "Regulacion_intercambio", rank = "AICc")
    modelo <- get.models(modelos, 1)[[1]]
    datos_pred <- rbind(temp, data_CR)
    datos_pred$fitted <- predict(modelo, newdata = datos_pred)
    resultado$Intercambio_AICC <- mean(datos_pred$fitted[datos_pred$Codigo_Pais == "CRI"], na.rm = TRUE)
  }, error = function(e) {
    cat("Error en AICC sin EF, iteración", i, ":", conditionMessage(e), "\n")
  })
  
  # BIC sin EF
  tryCatch({
    modelo <- lm(formula_con_cont, data = temp)
    options(na.action = "na.fail")
    modelos <- dredge(modelo, fixed = "Regulacion_intercambio", rank = "BIC")
    modelo <- get.models(modelos, 1)[[1]]
    datos_pred <- rbind(temp, data_CR)
    datos_pred$fitted <- predict(modelo, newdata = datos_pred)
    resultado$Intercambio_BIC <- mean(datos_pred$fitted[datos_pred$Codigo_Pais == "CRI"], na.rm = TRUE)
  }, error = function(e) {
    cat("Error en BIC sin EF, iteración", i, ":", conditionMessage(e), "\n")
  })
  
  # AICC con EF
  tryCatch({
    modelo <- lm(formula_con_cont_t, data = temp)
    options(na.action = "na.fail")
    modelos <- dredge(modelo, fixed = variables_fijas_validas, rank = "AICc")
    modelo <- get.models(modelos, 1)[[1]]
    datos_pred <- rbind(temp, data_CR)
    datos_pred$fitted <- predict(modelo, newdata = datos_pred)
    resultado$Intercambio_AICC_EF <- mean(datos_pred$fitted[datos_pred$Codigo_Pais == "CRI"], na.rm = TRUE)
  }, error = function(e) {
    cat("Error en AICC con EF, iteración", i, ":", conditionMessage(e), "\n")
  })
  
  # BIC con EF
  tryCatch({
    modelo <- lm(formula_con_cont_t, data = temp)
    options(na.action = "na.fail")
    modelos <- dredge(modelo, fixed = variables_fijas_validas, rank = "BIC")
    modelo <- get.models(modelos, 1)[[1]]
    datos_pred <- rbind(temp, data_CR)
    datos_pred$fitted <- predict(modelo, newdata = datos_pred)
    resultado$Intercambio_BIC_EF <- mean(datos_pred$fitted[datos_pred$Codigo_Pais == "CRI"], na.rm = TRUE)
  }, error = function(e) {
    cat("Error en BIC con EF, iteración", i, ":", conditionMessage(e), "\n")
  })
  
  return(resultado)
}

plan(multisession, workers = parallel::detectCores() - 1)

# Número de iteraciones
n_iter <- 1000

# Ejecutar en paralelo con semilla reproducible
resultados_lista <- future_map(
  1:n_iter,
  bootstrap_iter,
  .progress = TRUE,
  .options = furrr_options(seed = 472195)
)

# Convertir a data.table
resultados <- rbindlist(resultados_lista, fill = TRUE)

```

```{r, include= TRUE, eval=FALSE}
## Post Bootstrap
# Calcular el promedio
resultados$Intercambio <- rowMeans(resultados[, c("Intercambio_AICC",
                                                  "Intercambio_BIC",
                                                  "Intercambio_AICC_EF",
                                                  "Intercambio_BIC_EF")])


# Calcular estadísticas descriptivas
media <- mean(resultados$Intercambio, na.rm = TRUE)
mediana <- median(resultados$Intercambio, na.rm = TRUE)


# Formatear con dos decimales
media_fmt <- sprintf("%.2f", media)
mediana_fmt <- sprintf("%.2f", mediana)

# Mostrar resultados
cat("Media:", media_fmt, "\n")
cat("Mediana:", mediana_fmt, "\n")


#Obtener intervalo de confianza del 95% para Intercambio
c1 <- round(quantile(resultados$Intercambio, 0.025), 2)
c2 <- round(quantile(resultados$Intercambio, 0.975), 2)
ci <- quantile(resultados$Intercambio, probs = c(0.025, 0.975))

# Mostrar los valores
cat("IC 95%:", c1, "-", c2, "\n")

#Generar gráfico para Intercambio

p <- ggplot(resultados, aes(x = Intercambio)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "lightblue", color = "black") +
  theme_minimal(base_size = 14) +
  labs(
    x = "Comisión de intercambio (%)",
    y = "Densidad",
    title = "Distribución de la Comisión de Intercambio"
  ) +
  annotate("text", x = media, y = 0.5, label = paste0("Media: ", round(media, 2), "%"), size = 5) +
  annotate("text", x = mediana, y = 0.4, label = paste0("Mediana: ", round(mediana, 2), "%"), size = 5) +
  annotate("text", x = mean(ci), y = 0.3,
           label = paste0("IC 95%: [", round(ci[1], 2), "%, ", round(ci[2], 2), "%]"),
           size = 5)
p

```

```{r, include= TRUE, eval=FALSE}
end_time <- Sys.time()
duration <- end_time - start_time
cat("Tiempo total de ejecución:", duration, "\n")
```

### Regresión Robusta

::: {style="text-align: justify"}
El método de regresión robusta permite analizar cambios en los resultados obtenidos producto de observaciones atípicas. Esto se consigue mediante la asignación de ponderadores a cada observación, que varían en sentido inverso a su nivel de influencia (medida por el tamaño de sus residuos), que cambian de manera iterativa hasta alcanzar un punto de convergencia.
:::

#### MRR_aicc

##### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRR_aicc <- dredge(modelo_completo, rank = "AICc", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MRR_aicc <- get.models(MRR_aicc, 1)[[1]]

# Extraer variables seleccionadas
vars_seleccionadas <- names(coef(MRR_aicc))[-1]  # sin intercepto

# Ajustar regresión robusta con las variables seleccionadas

# Paso 1: Ajuste robusto con Huber
formula_robusta <- as.formula(paste("Intercambio ~", paste(vars_seleccionadas, collapse = " + ")))
modelo_huber <- rlm(formula_robusta, data = datos_int_con_cont,  psi = psi.huber)


# Paso 2: Ajuste robusto con Bisquare usando coef iniciales de Huber
modelo_bisquare <- rlm(formula_robusta, data = datos_int_con_cont, psi = psi.bisquare, init = coef(modelo_huber))

MRR_aicc <- modelo_bisquare

#Visualización
stargazer(MRR_aicc, type = "html", 
          title = "Regresión robusta - Modelo AICc sin efectos fijos de tiempo", 
          column.labels = c("MCO"), model.names = F, 
          dep.var.labels = c("Tasa intercambio"), dep.var.caption = "Variable dependiente")

```

##### Predicción

```{r, include=TRUE}
pred_int_MRR_aicc=predict(MRR_aicc, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MRR_aicc`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Regresión robusta - Modelo AICc sin efectos fijos de tiempo")
pronostico <- c(pred_int_MRR_aicc)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRR_aicc, modelo, pred_int_MRR_aicc, pronostico, modelo_bisquare, modelo_huber, vars_seleccionadas, formula_robusta)
```

#### MRR_bic

##### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRR_bic <- dredge(modelo_completo, rank = "BIC", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MRR_bic <- get.models(MRR_bic, 1)[[1]]

# Extraer variables seleccionadas
vars_seleccionadas <- names(coef(MRR_bic))[-1]  # sin intercepto

# Ajustar regresión robusta con las variables seleccionadas

# Paso 1: Ajuste robusto con Huber
formula_robusta <- as.formula(paste("Intercambio ~", paste(vars_seleccionadas, collapse = " + ")))
modelo_huber <- rlm(formula_robusta, data = datos_int_con_cont,  psi = psi.huber)


# Paso 2: Ajuste robusto con Bisquare usando coef iniciales de Huber
modelo_bisquare <- rlm(formula_robusta, data = datos_int_con_cont, psi = psi.bisquare, init = coef(modelo_huber))

MRR_bic <- modelo_bisquare

#Visualización
stargazer(MRR_bic, type = "html", 
          title = "Regresión robusta - Modelo BIC sin efectos fijos de tiempo", 
          column.labels = c("MCO"), model.names = F, 
          dep.var.labels = c("Tasa intercambio"), dep.var.caption = "Variable dependiente")

```

##### Predicción

```{r, include=TRUE}

pred_int_MRR_bic=predict(MRR_bic, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MRR_bic`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Regresión robusta - Modelo BIC sin efectos fijos de tiempo")
pronostico <- c(pred_int_MRR_bic)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRR_bic, modelo, pred_int_MRR_bic, pronostico, modelo_bisquare, modelo_huber, vars_seleccionadas, formula_robusta)
```

#### MRR_aicc_et

##### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRR_aicc_et <- dredge(modelo_completo_et, rank = "AICc", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MRR_aicc_et <- get.models(MRR_aicc_et, 1)[[1]]

# Obtener coeficientes
coeficientes <- coef(MRR_aicc_et)

#Filtrar nombres de coeficientes que no sean NA
coeficientes_validos <- names(coeficientes[!is.na(coeficientes)])

# Extraer variables seleccionadas
vars_seleccionadas <- names(coef(MRR_aicc_et))[-1]  # sin intercepto
vars_seleccionadas <- intersect(coeficientes_validos, vars_seleccionadas)

# Ajustar regresión robusta con las variables seleccionadas
# Paso 1: Ajuste robusto con Huber
formula_robusta <- as.formula(paste("Intercambio ~", paste(vars_seleccionadas, collapse = " + ")))

modelo_huber <- rlm(formula_robusta, data = datos_int_con_cont,  psi = psi.huber)

# Paso 2: Ajuste robusto con Bisquare usando coef iniciales de Huber
modelo_bisquare <- rlm(formula_robusta, data = datos_int_con_cont, psi = psi.bisquare, init = coef(modelo_huber))

MRR_aicc_et <- modelo_bisquare

#Visualización
stargazer(MRR_aicc_et, type = "html", 
          title = "Regresión robusta - Modelo AICc con efectos fijos de tiempo", 
          column.labels = c("MCO"), model.names = F, 
          dep.var.labels = c("Tasa intercambio"), dep.var.caption = "Variable dependiente")

```

##### Predicción

```{r, include=TRUE}
pred_int_MRR_aicc_et=predict(MRR_aicc_et, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MRR_aicc_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Regresión robusta - Modelo AICc con efectos fijos de tiempo")
pronostico <- c(pred_int_MRR_aicc_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRR_aicc_et, modelo, pred_int_MRR_aicc_et, pronostico, modelo_bisquare, modelo_huber, vars_seleccionadas, formula_robusta, coeficientes, coeficientes_validos)
```

#### MRR_bic_et

##### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRR_bic_et <- dredge(modelo_completo_et, rank = "BIC", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MRR_bic_et <- get.models(MRR_bic_et, 1)[[1]]

# Obtener coeficientes
coeficientes <- coef(MRR_bic_et)

#Filtrar nombres de coeficientes que no sean NA
coeficientes_validos <- names(coeficientes[!is.na(coeficientes)])

# Extraer variables seleccionadas
vars_seleccionadas <- names(coef(MRR_bic_et))[-1]  # sin intercepto
vars_seleccionadas <- intersect(coeficientes_validos, vars_seleccionadas)

# Ajustar regresión robusta con las variables seleccionadas
# Paso 1: Ajuste robusto con Huber
formula_robusta <- as.formula(paste("Intercambio ~", paste(vars_seleccionadas, collapse = " + ")))

modelo_huber <- rlm(formula_robusta, data = datos_int_con_cont,  psi = psi.huber)

# Paso 2: Ajuste robusto con Bisquare usando coef iniciales de Huber
modelo_bisquare <- rlm(formula_robusta, data = datos_int_con_cont, psi = psi.bisquare, init = coef(modelo_huber))

MRR_bic_et <- modelo_bisquare

#Visualización
stargazer(MRR_bic_et, type = "html", 
          title = "Regresión robusta - Modelo BIC con efectos fijos de tiempo", 
          column.labels = c("MCO"), model.names = F, 
          dep.var.labels = c("Tasa intercambio"), dep.var.caption = "Variable dependiente")

```

##### Predicción

```{r, include=TRUE}
pred_int_MRR_bic_et <- predict(MRR_bic_et, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MRR_bic_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Regresión robusta - Modelo BIC con efectos fijos de tiempo")
pronostico <- c(pred_int_MRR_bic_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRR_bic_et, modelo, pred_int_MRR_bic_et, pronostico, modelo_bisquare, modelo_huber, vars_seleccionadas, formula_robusta, coeficientes, coeficientes_validos)
```

### Omitir al Espacio Económico Europeo (EEE) en la estimación

::: {style="text-align: justify"}
Los países miembros del Espacio Económico (EEE) Europeo sobresalen por contar con las comisiones de intercambio máximas más bajas de la muestra. Si bien existe heterogeneidad en las características sociales, económicas y financieras de sus miembros, es pertinente analizar hasta qué punto las estimaciones obtenidas para Costa Rica pueden ser influenciadas por la regulación de este bloque.

Para esto, se procede a estimar el modelo principal, pero sin considerar los 30 países que conforman el EEE. Además, se enfoca el análisis para el caso de la comisión de intercambio, que es la comisión regulada en el EEE.
:::

```{r, include=TRUE}
#Datos

#Se filtran los datos para excluir países del EEE
datos_sin_EEE <- datos[datos$EEE != 1,]
datos_sin_EEE <- datos_sin_EEE %>%
  select(all_of(variables_interes)) %>%
  na.omit()
```

#### MOEEE_aicc \*\*

##### Estimación

```{r, include=TRUE, results='asis'}
# Ajustar el modelo
modelo_completo_sin_eee <- lm(formula_con_cont, data = datos_sin_EEE)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MOEEE_aicc <- dredge(modelo_completo_sin_eee, rank = "AICc", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MOEEE_aicc <- get.models(MOEEE_aicc, 1)[[1]]
MOEEE_aicc

# Errores robustos
MOEEE_aicc_corr_hc <- coeftest(MOEEE_aicc, vcov. = vcovHC(MOEEE_aicc, "HC1", adjust = TRUE, sandwich = TRUE))

MOEEE_aicc_corr_hac <- coeftest(MOEEE_aicc, vcov. = vcovHAC(MOEEE_aicc, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MOEEE_aicc, MOEEE_aicc_corr_hc,
          type = "html", 
          title = "Estimación omitiendo EEE - Modelo AICc sin efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MOEEE_aicc <- predict(MOEEE_aicc, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MOEEE_aicc`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación omitiendo EEE - Modelo AICc sin efectos fijos de tiempo")
pronostico <- c(pred_int_MOEEE_aicc)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MOEEE_aicc, modelo, pred_int_MOEEE_aicc, pronostico, MOEEE_aicc_corr_hc, MOEEE_aicc_corr_hac)
```

#### MOEEE_bic \*\*

##### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MOEEE_bic <- dredge(modelo_completo_sin_eee, rank = "BIC", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MOEEE_bic <- get.models(MOEEE_bic, 1)[[1]]
MOEEE_bic

# Errores robustos
MOEEE_bic_corr_hc <- coeftest(MOEEE_bic, vcov. = vcovHC(MOEEE_bic, "HC1", adjust = TRUE, sandwich = TRUE))

MOEEE_bic_corr_hac <- coeftest(MOEEE_bic, vcov. = vcovHAC(MOEEE_bic, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MOEEE_bic, MOEEE_bic_corr_hc, type = "html",
          title = "Estimación omitiendo EEE - Modelo BIC sin efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F, 
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")

```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MOEEE_bic <- predict(MOEEE_bic, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MOEEE_bic`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación omitiendo EEE - Modelo BIC sin efectos fijos de tiempo")
pronostico <- c(pred_int_MOEEE_bic)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MOEEE_bic, modelo, pred_int_MOEEE_bic, pronostico, MOEEE_bic_corr_hc, MOEEE_bic_corr_hac, modelo_completo_sin_eee)
```

#### MOEEE_aicc_et \*\*

##### Estimación

```{r, include=TRUE, results='asis'}
# Ajustar el modelo
modelo_completo_sin_eee_et <- lm(formula_con_cont_t, data = datos_sin_EEE)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MOEEE_aicc_et <- dredge(modelo_completo_sin_eee_et, rank = "AICc", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MOEEE_aicc_et <- get.models(MOEEE_aicc_et, 1)[[1]]
MOEEE_aicc_et

# Errores robustos
MOEEE_aicc_et_corr_hc <- coeftest(MOEEE_aicc_et, vcov. = vcovHC(MOEEE_aicc_et, "HC1", adjust = TRUE, sandwich = TRUE))

MOEEE_aicc_et_corr_hac <- coeftest(MOEEE_aicc_et, vcov. = vcovHAC(MOEEE_aicc_et, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MOEEE_aicc_et, MOEEE_aicc_et_corr_hc,
          type = "html", 
          title = "Estimación omitiendo EEE - Modelo AICc con efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MOEEE_aicc_et <- predict(MOEEE_aicc_et, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MOEEE_aicc_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación omitiendo EEE - Modelo AICc con efectos fijos de tiempo")
pronostico <- c(pred_int_MOEEE_aicc_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MOEEE_aicc_et, modelo, pred_int_MOEEE_aicc_et, pronostico, MOEEE_aicc_et_corr_hc, MOEEE_aicc_et_corr_hac)
```

#### MOEEE_bic_et \*\*

##### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MOEEE_bic_et <- dredge(modelo_completo_sin_eee_et, rank = "BIC", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MOEEE_bic_et <- get.models(MOEEE_bic_et, 1)[[1]]
MOEEE_bic_et

# Errores robustos
MOEEE_bic_et_corr_hc <- coeftest(MOEEE_bic_et, vcov. = vcovHC(MOEEE_bic_et, "HC1", adjust = TRUE, sandwich = TRUE))

MOEEE_bic_et_corr_hac <- coeftest(MOEEE_bic_et, vcov. = vcovHAC(MOEEE_bic_et, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MOEEE_bic_et, MOEEE_bic_et_corr_hc, type = "html",
          title = "Estimación omitiendo EEE - Modelo BIC con efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F, 
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")

```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MOEEE_bic_et <- predict(MOEEE_bic_et, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MOEEE_bic_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación omitiendo EEE - Modelo BIC con efectos fijos de tiempo")
pronostico <- c(pred_int_MOEEE_bic_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MOEEE_bic_et, modelo, pred_int_MOEEE_bic_et, pronostico, MOEEE_bic_et_corr_hc, MOEEE_bic_et_corr_hac, modelo_completo_sin_eee_et, datos_sin_EEE)
```

### Restringiendo la muestra a solo países regulados

::: {style="text-align: justify"}
Como una prueba adicional, se estima el modelo principal, pero luego de restringir la muestra a los países que al igual que Costa Rica, son regulados.
:::

```{r, include=FALSE}
#Datos

#Se filtran los datos para excluir países no regulados
datos_reg <- datos[datos$Regulacion_intercambio == 1,]
datos_reg <- datos_reg %>%
  select(all_of(variables_interes)) %>%
  na.omit()
```

#### MPR_aicc

##### Estimación

```{r, include=TRUE, results='asis'}
# Ajustar el modelo
modelo_completo_reg <- lm(formula_con_cont, data = datos_reg)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MPR_aicc <- dredge(modelo_completo_reg, rank = "AICc", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MPR_aicc <- get.models(MPR_aicc, 1)[[1]]
MPR_aicc

# Errores robustos
MPR_aicc_corr_hc <- coeftest(MPR_aicc, vcov. = vcovHC(MPR_aicc, "HC1", adjust = TRUE, sandwich = TRUE))

MPR_aicc_corr_hac <- coeftest(MPR_aicc, vcov. = vcovHAC(MPR_aicc, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MPR_aicc, MPR_aicc_corr_hc,
          type = "html", 
          title = "Estimación restringiendo muestra a países regulados - Modelo AICc sin efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MPR_aicc <- predict(MPR_aicc, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MPR_aicc`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación restringiendo muestra a países regulados - Modelo AICc sin efectos fijos de tiempo")
pronostico <- c(pred_int_MPR_aicc)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MPR_aicc, modelo, pred_int_MPR_aicc, pronostico, MPR_aicc_corr_hc, MPR_aicc_corr_hac)
```

#### MPR_bic

##### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MPR_bic <- dredge(modelo_completo_reg, rank = "BIC", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MPR_bic <- get.models(MPR_bic, 1)[[1]]
MPR_bic

# Errores robustos
MPR_bic_corr_hc <- coeftest(MPR_bic, vcov. = vcovHC(MPR_bic, "HC1", adjust = TRUE, sandwich = TRUE))

MPR_bic_corr_hac <- coeftest(MPR_bic, vcov. = vcovHAC(MPR_bic, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MPR_bic, MPR_bic_corr_hc, type = "html",
          title = "Estimación restringiendo muestra a países regulados - Modelo BIC sin efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F, 
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")

```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MPR_bic <- predict(MPR_bic, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MPR_bic`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación restringiendo muestra a países regulados - Modelo BIC sin efectos fijos de tiempo")
pronostico <- c(pred_int_MPR_bic)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MPR_bic, modelo, pred_int_MPR_bic, pronostico, MPR_bic_corr_hc, MPR_bic_corr_hac, modelo_completo_reg)
```

#### MPR_aicc_et

##### Estimación

```{r, include=TRUE, results='asis'}
# Ajustar el modelo
modelo_completo_et_reg <- lm(formula_con_cont_t, data = datos_reg)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MPR_aicc_et <- dredge(modelo_completo_et_reg, rank = "AICc", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MPR_aicc_et <- get.models(MPR_aicc_et, 1)[[1]]
MPR_aicc_et

# Errores robustos
MPR_aicc_et_corr_hc <- coeftest(MPR_aicc_et, vcov. = vcovHC(MPR_aicc_et, "HC1", adjust = TRUE, sandwich = TRUE))

MPR_aicc_et_corr_hac <- coeftest(MPR_aicc_et, vcov. = vcovHAC(MPR_aicc_et, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MPR_aicc_et, MPR_aicc_et_corr_hc,
          type = "html", 
          title = "Estimación restringiendo muestra a países regulados - Modelo AICc con efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MPR_aicc_et <- predict(MPR_aicc_et, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MPR_aicc_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación restringiendo muestra a países regulados - Modelo AICc con efectos fijos de tiempo")
pronostico <- c(pred_int_MPR_aicc_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MPR_aicc_et, modelo, pred_int_MPR_aicc_et, pronostico, MPR_aicc_et_corr_hc, MPR_aicc_et_corr_hac)
```

#### MPR_bic_et

##### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MPR_bic_et <- dredge(modelo_completo_et_reg, rank = "BIC", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MPR_bic_et <- get.models(MPR_bic_et, 1)[[1]]
MPR_bic_et

# Errores robustos
MPR_bic_et_corr_hc <- coeftest(MPR_bic_et, vcov. = vcovHC(MPR_bic_et, "HC1", adjust = TRUE, sandwich = TRUE))

MPR_bic_et_corr_hac <- coeftest(MPR_bic_et, vcov. = vcovHAC(MPR_bic_et, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MPR_bic_et, MPR_bic_et_corr_hc,
          type = "html", 
          title = "Estimación restringiendo muestra a países regulados - Modelo BIC con efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MPR_bic_et <- predict(MPR_bic_et, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MPR_bic_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación restringiendo muestra a países regulados - Modelo BIC sin efectos fijos de tiempo")
pronostico <- c(pred_int_MPR_bic_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MPR_bic_et, modelo, pred_int_MPR_bic_et, pronostico, MPR_bic_et_corr_hc, MPR_bic_et_corr_hac, datos_reg, modelo_completo_et_reg)
```

### Ampliando la definición de regulación de intercambio o adquirencia

::: {style="text-align: justify"}
En el análisis principal, la variable regulación hace referencia a regulación en las comisiones de intercambio. Esto se debe a que son relativamente pocos los países que regulan adquirencia. Sin embargo, es posible generar una variable dicotómica que haga referencia a que el país en cuestión tiene regulación sobre intercambio o adquirencia.
:::

```{r, include=TRUE}
# Arreglo de datos
datos_reg_amp <- datos

# Crear una nueva variable lógica que sea TRUE si se cumple alguna de las dos condiciones

datos_reg_amp$Regulacion_inter_adq <- ifelse(datos_reg_amp$Regulacion_intercambio == 1 | datos_reg_amp$Regulacion_adquirencia == 1, 1, 0)

#Se filtran los datos
datos_reg_amp <- datos_reg_amp %>%
  select(all_of(variables_interes), Regulacion_inter_adq) %>% 
  select(-Regulacion_intercambio) %>%
  na.omit()

# Crear una nueva variable lógica que sea TRUE si se cumple alguna de las dos condiciones
data_CR_reg_amp <- data_CR
data_CR_reg_amp$Regulacion_inter_adq <- ifelse(data_CR_reg_amp$Regulacion_intercambio == 1 | data_CR_reg_amp$Regulacion_adquirencia == 1, 1, 0)

# Incluir las variables de control sin controles de tiempo
formula_con_cont_amp <- as.formula(
  paste("Intercambio ~ Regulacion_inter_adq +", paste(controles, collapse = " + "))
)

# Incluir las variables de control con controles de tiempo
formula_con_cont_t_amp <- as.formula(
  paste("Intercambio ~ Regulacion_inter_adq +", paste(controles, collapse = " + ")," + ", paste(ef_tiempo, collapse = " + "))
)

# Variables fijas
variables_fijas_amp <- c("Regulacion_inter_adq", paste0("Año_", 2010:2025))
```

#### MADR_aicc

##### Estimación

```{r, include=TRUE, results='asis'}
# Ajustar el modelo
modelo_completo_reg_amp <- lm(formula_con_cont_amp, data = datos_reg_amp)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MADR_aicc <- dredge(modelo_completo_reg_amp, rank = "AICc", fixed = "Regulacion_inter_adq")

# Extraer el mejor modelo según AICc
MADR_aicc <- get.models(MADR_aicc, 1)[[1]]
MADR_aicc

# Errores robustos
MADR_aicc_corr_hc <- coeftest(MADR_aicc, vcov. = vcovHC(MADR_aicc, "HC1", adjust = TRUE, sandwich = TRUE))

MADR_aicc_corr_hac <- coeftest(MADR_aicc, vcov. = vcovHAC(MADR_aicc, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MADR_aicc, MADR_aicc_corr_hc,
          type = "html", 
          title = "Estimación ampliando definición de regulación - Modelo AICc sin efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MADR_aicc <- predict(MADR_aicc, newdata = data_CR_reg_amp, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MADR_aicc`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación ampliando definición de regulación - Modelo AICc sin efectos fijos de tiempo")
pronostico <- c(pred_int_MADR_aicc)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MADR_aicc, modelo, pred_int_MADR_aicc, pronostico, MADR_aicc_corr_hc, MADR_aicc_corr_hac)
```

#### MADR_bic

##### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MADR_bic <- dredge(modelo_completo_reg_amp, rank = "BIC", fixed = "Regulacion_inter_adq")

# Extraer el mejor modelo según AICc
MADR_bic <- get.models(MADR_bic, 1)[[1]]
MADR_bic

# Errores robustos
MADR_bic_corr_hc <- coeftest(MADR_bic, vcov. = vcovHC(MADR_bic, "HC1", adjust = TRUE, sandwich = TRUE))

MADR_bic_corr_hac <- coeftest(MADR_bic, vcov. = vcovHAC(MADR_bic, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MADR_bic, MADR_bic_corr_hc, type = "html",
          title = "Estimación ampliando definición de regulación - Modelo BIC sin efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F, 
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")

```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MADR_bic <- predict(MADR_bic, newdata = data_CR_reg_amp, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MADR_bic`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación ampliando definición de regulación - Modelo BIC sin efectos fijos de tiempo")
pronostico <- c(pred_int_MADR_bic)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MADR_bic, modelo, pred_int_MADR_bic, pronostico, MADR_bic_corr_hc, MADR_bic_corr_hac, modelo_completo_reg_amp, formula_con_cont_amp)
```

#### MADR_aicc_et \*\*

##### Estimación

```{r, include=TRUE, results='asis'}
# Ajustar el modelo
modelo_completo_t_reg_amp <- lm(formula_con_cont_t_amp, data = datos_reg_amp)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MADR_aicc_et <- dredge(modelo_completo_t_reg_amp, rank = "AICc", fixed = variables_fijas_amp)

# Extraer el mejor modelo según AICc
MADR_aicc_et <- get.models(MADR_aicc_et, 1)[[1]]
MADR_aicc_et

# Errores robustos
MADR_aicc_et_corr_hc <- coeftest(MADR_aicc_et, vcov. = vcovHC(MADR_aicc_et, "HC1", adjust = TRUE, sandwich = TRUE))

MADR_aicc_et_corr_hac <- coeftest(MADR_aicc_et, vcov. = vcovHAC(MADR_aicc_et, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MADR_aicc_et, MADR_aicc_et_corr_hc,
          type = "html", 
          title = "Estimación ampliando definición de regulación - Modelo AICc con efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MADR_aicc_et <- predict(MADR_aicc_et, newdata = data_CR_reg_amp, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MADR_aicc_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación ampliando definición de regulación - Modelo AICc con efectos fijos de tiempo")
pronostico <- c(pred_int_MADR_aicc_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MADR_aicc_et, modelo, pred_int_MADR_aicc_et, pronostico, MADR_aicc_et_corr_hc, MADR_aicc_et_corr_hac)
```

#### MADR_bic_et

##### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MADR_bic_et <- dredge(modelo_completo_t_reg_amp, rank = "BIC", fixed = variables_fijas_amp)

# Extraer el mejor modelo según AICc
MADR_bic_et <- get.models(MADR_bic_et, 1)[[1]]
MADR_bic_et

# Errores robustos
MADR_bic_et_corr_hc <- coeftest(MADR_bic_et, vcov. = vcovHC(MADR_bic_et, "HC1", adjust = TRUE, sandwich = TRUE))

MADR_bic_et_corr_hac <- coeftest(MADR_bic_et, vcov. = vcovHAC(MADR_bic_et, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MADR_bic_et, MADR_bic_et_corr_hc,
          type = "html", 
          title = "Estimación ampliando definición de regulación - Modelo BIC con efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MADR_bic_et <- predict(MADR_bic_et, newdata = data_CR_reg_amp, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MADR_bic_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación ampliando definición de regulación - Modelo BIC con efectos fijos de tiempo")
pronostico <- c(pred_int_MADR_bic_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MADR_bic_et, modelo, pred_int_MADR_bic_et, pronostico, MADR_bic_et_corr_hc, MADR_bic_et_corr_hac, data_CR_reg_amp, datos_reg_amp, modelo_completo_t_reg_amp, formula_con_cont_t_amp, variables_fijas_amp)
```

### Regresión que utilice las características observadas más recientes para cada país

::: {style="text-align: justify"}
La estimación principal utiliza las características observadas para los países al momento de la fijación de determinado porcentaje de comisión. Una forma alternativa de estimar el modelo es considerar para cada país las comisiones de intercambio o de adquirencia vigentes y las características observables más recientes, que en la mayoría de los países corresponde a información para el 2024.
:::

```{r, include=TRUE}
#Datos
datos_2024 <- read_excel("D:/Base_datos_adquirencia_intercambio_2024.xlsx")

datos_2024 <- datos_2024 %>% 
  dplyr :: select(-Adquirencia, -Regulacion_adquirencia) %>% 
  na.omit
```

#### MCR_aicc

##### Estimación

```{r, include=TRUE, results='asis'}
# Ajustar el modelo
modelo_completo_2024 <- lm(formula_con_cont, data = datos_2024)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MCR_aicc <- dredge(modelo_completo_2024, rank = "AICc", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MCR_aicc <- get.models(MCR_aicc, 1)[[1]]
MCR_aicc

# Errores robustos
MCR_aicc_corr_hc <- coeftest(MCR_aicc, vcov. = vcovHC(MCR_aicc, "HC1", adjust = TRUE, sandwich = TRUE))

MCR_aicc_corr_hac <- coeftest(MCR_aicc, vcov. = vcovHAC(MCR_aicc, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MCR_aicc, MCR_aicc_corr_hc,
          type = "html", 
          title = "Estimación usando características más recientes - Modelo AICc sin efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MCR_aicc <- predict(MCR_aicc, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MCR_aicc`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación usando características más recientes - Modelo AICc sin efectos fijos de tiempo")
pronostico <- c(pred_int_MCR_aicc)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MCR_aicc, modelo, pred_int_MCR_aicc, pronostico, MCR_aicc_corr_hc, MCR_aicc_corr_hac)
```

#### MCR_bic

##### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MCR_bic <- dredge(modelo_completo_2024, rank = "BIC", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MCR_bic <- get.models(MCR_bic, 1)[[1]]
MCR_bic

# Errores robustos
MCR_bic_corr_hc <- coeftest(MCR_bic, vcov. = vcovHC(MCR_bic, "HC1", adjust = TRUE, sandwich = TRUE))

MCR_bic_corr_hac <- coeftest(MCR_bic, vcov. = vcovHAC(MCR_bic, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MCR_bic, MCR_bic_corr_hc,
          type = "html", 
          title = "Estimación usando características más recientes - Modelo BIC sin efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MCR_bic <- predict(MCR_bic, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MCR_bic`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación usando características más recientes - Modelo BIC sin efectos fijos de tiempo")
pronostico <- c(pred_int_MCR_bic)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MCR_bic, modelo, pred_int_MCR_bic, pronostico, MCR_bic_corr_hc, MCR_bic_corr_hac, datos_2024, modelo_completo_2024)
```

### Regresión con variables propuestas por el Fondo Monetario Internacional (FMI)

::: {style="text-align: justify"}
El Fondo Monetario Internacional (FMI) publicó en diciembre de 2023 el estudio “Determinants of Card Payments Fees in A Cross-Country Analysis”, donde se replica y amplía el análisis de comparación internacional que realiza el BCCR para establecer los topes máximos de comisiones. La forma de ampliar el estudio fue incorporar las variables margen de interés neto y el grado de concentración de activos bancarios.
:::

```{r, include=TRUE}
# Datos
datos_fmi <- datos %>%
  select(all_of(variables_interes), c("Margen_neto", "Concentracion_bancaria")) %>%
  na.omit()

#Definir el conjunto de posibles variables de control

controles_fmi <- c("Cajero_Autom", "Crimen", "Cuenta", "Credito", "Debito", "Densidad_Pob", "Densidad_Neg", "PIB_pc_PPA", "ICH", "Pago_digital", "Valor_prom", "Europa", "Asia_Este", "Am_lat", "Margen_neto", "Concentracion_bancaria")

# Incluir las variables de control
formula_con_cont_fmi <- as.formula(
  paste("Intercambio ~ Regulacion_intercambio +", paste(controles_fmi, collapse = " + "))
)

# Incluir las variables de control
formula_con_cont_t_fmi <- as.formula(
  paste("Intercambio ~ Regulacion_intercambio +", paste(controles_fmi, collapse = " + ")," + ", paste(ef_tiempo, collapse = " + "))
)

```

#### MFMI_aicc

##### Estimación

```{r, include=TRUE, results='asis'}
# Ajustar el modelo
modelo_completo_fmi <- lm(formula_con_cont_fmi, data = datos_fmi)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MFMI_aicc <- dredge(modelo_completo_fmi, rank = "AICc", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MFMI_aicc <- get.models(MFMI_aicc, 1)[[1]]
MFMI_aicc

# Errores robustos
MFMI_aicc_corr_hc <- coeftest(MFMI_aicc, vcov. = vcovHC(MFMI_aicc, "HC1", adjust = TRUE, sandwich = TRUE))

MFMI_aicc_corr_hac <- coeftest(MFMI_aicc, vcov. = vcovHAC(MFMI_aicc, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MFMI_aicc, MFMI_aicc_corr_hc,
          type = "html", 
          title = "Réplica análisis FMI - Modelo AICc sin efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MFMI_aicc <- predict(MFMI_aicc, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MFMI_aicc`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Réplica análisis FMI - Modelo AICc sin efectos fijos de tiempo")
pronostico <- c(pred_int_MFMI_aicc)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MFMI_aicc, modelo, pred_int_MFMI_aicc, pronostico, MFMI_aicc_corr_hac, MFMI_aicc_corr_hc)
```

#### MFMI_bic

##### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MFMI_bic <- dredge(modelo_completo_fmi, rank = "BIC", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MFMI_bic <- get.models(MFMI_bic, 1)[[1]]
MFMI_bic

# Errores robustos
MFMI_bic_corr_hc <- coeftest(MFMI_bic, vcov. = vcovHC(MFMI_bic, "HC1", adjust = TRUE, sandwich = TRUE))

MFMI_bic_corr_hac <- coeftest(MFMI_bic, vcov. = vcovHAC(MFMI_bic, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MFMI_bic, MFMI_bic_corr_hc,
          type = "html", 
          title = "Réplica análisis FMI - Modelo BIC sin efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MFMI_bic <- predict(MFMI_bic, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MFMI_bic`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Réplica análisis FMI - Modelo BIC sin efectos fijos de tiempo")
pronostico <- c(pred_int_MFMI_bic)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MFMI_bic, modelo, pred_int_MFMI_bic, pronostico, MFMI_bic_corr_hac, MFMI_bic_corr_hc, modelo_completo_fmi, formula_con_cont_fmi)
```

#### MFMI_aicc_et

##### Estimación

```{r , include=TRUE, results='asis'}
# Ajustar el modelo con controles
modelo_completo_fmi_et <- lm(formula_con_cont_t_fmi, data = datos_fmi)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MFMI_aicc_et <- dredge(modelo_completo_fmi_et, rank = "AICc", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MFMI_aicc_et <- get.models(MFMI_aicc_et, 1)[[1]]
MFMI_aicc_et

# Errores robustos
MFMI_aicc_et_corr_hc <- coeftest(MFMI_aicc_et, vcov. = vcovHC(MFMI_aicc_et, "HC1", adjust = TRUE, sandwich = TRUE))

MFMI_aicc_et_corr_hac <- coeftest(MFMI_aicc_et, vcov. = vcovHAC(MFMI_aicc_et, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MFMI_aicc_et, MFMI_aicc_et_corr_hc,
          type = "html", 
          title = "Réplica análisis FMI - Modelo AICc con efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MFMI_aicc_et <- predict(MFMI_aicc_et, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MFMI_aicc_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Réplica análisis FMI - Modelo AICc con efectos fijos de tiempo")
pronostico <- c(pred_int_MFMI_aicc_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MFMI_aicc_et, modelo, pred_int_MFMI_aicc_et, pronostico, MFMI_aicc_et_corr_hac, MFMI_aicc_et_corr_hc)
```

#### MFMI_bic_et \*\*

##### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MFMI_bic_et <- dredge(modelo_completo_fmi_et, rank = "BIC", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MFMI_bic_et <- get.models(MFMI_bic_et, 1)[[1]]
MFMI_bic_et

# Errores robustos
MFMI_bic_et_corr_hc <- coeftest(MFMI_bic_et, vcov. = vcovHC(MFMI_bic_et, "HC1", adjust = TRUE, sandwich = TRUE))

MFMI_bic_et_corr_hac <- coeftest(MFMI_bic_et, vcov. = vcovHAC(MFMI_bic_et, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MFMI_bic_et, MFMI_bic_et_corr_hc,
          type = "html", 
          title = "Réplica análisis FMI - Modelo BIC con efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MFMI_bic_et <- predict(MFMI_bic_et, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MFMI_bic_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Réplica análisis FMI - Modelo BIC con efectos fijos de tiempo")
pronostico <- c(pred_int_MFMI_bic_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MFMI_bic_et, modelo, pred_int_MFMI_bic_et, pronostico, MFMI_bic_et_corr_hac, MFMI_bic_et_corr_hc, datos_fmi, modelo_completo_fmi_et, controles_fmi, formula_con_cont_t_fmi)
```

### Regresión con variables instrumentales

::: {style="text-align: justify"}
Una posible observación es que las comisiones de intercambio y adquirencia afectan la adopción de tarjetas de pago, y estas a su vez afectan las comisiones, lo que introduciría un sesgo de endogeneidad en las estimaciones. Para atender esta observación se recurre a dos enfoques.

El primer enfoque consiste en estimar la especificación principal, pero sin considerar como variables explicativas la información relativa directamente al uso de tarjetas.

Para el segundo enfoque, se plantea un modelo de variables instrumentales. En este caso, se instrumenta el valor promedio de la transacción realizada con tarjetas de pago con el Índice de Capital Humano.
:::

#### Enfoque1: regresión sin variables relacionadas al uso de tarjetas

::: {style="text-align: justify"}
En esta regresión se estima la especificación principal, pero sin considerar como variables explicativas la información relativa directamente al uso de tarjetas. De esta manera, se excluye el porcentaje de adultos con tarjeta de crédito o débito, el porcentaje de adultos con cuenta en una institución financiera, el porcentaje de adultos que realizaron o recibieron un pago digital, el número de cajeros automáticos, y el valor de la transacción promedio realizada con tarjetas de pago.
:::

```{r, include=FALSE}
#Definir el conjunto de posibles variables de control

controles_red <- c("Crimen", "Densidad_Pob", "Densidad_Neg", "PIB_pc_PPA", "ICH", "Europa", "Asia_Este", "Am_lat")

# Incluir las variables de control
formula_con_cont_red <- as.formula(
  paste("Intercambio ~ Regulacion_intercambio +", paste(controles_red, collapse = " + "))
)

# Incluir las variables de control
formula_con_cont_t_red <- as.formula(
  paste("Intercambio ~ Regulacion_intercambio +", paste(controles_red, collapse = " + ")," + ", paste(ef_tiempo, collapse = " + "))
)

#Datos sin variables relacionadas al uso de tarjetas
variables_interes_red <- c(variables_constantes, controles_red, ef_tiempo)

datos_red <- datos %>%
  select(all_of(variables_interes_red)) %>%
  na.omit()
```

##### MSVR_aicc

###### Estimación

```{r, include=TRUE, results='asis'}
# Ajustar el modelo
modelo_completo_red <- lm(formula_con_cont_red, data = datos_red)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MSVR_aicc <- dredge(modelo_completo_red, rank = "AICc", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MSVR_aicc <- get.models(MSVR_aicc, 1)[[1]]
MSVR_aicc

# Errores robustos
MSVR_aicc_corr_hc <- coeftest(MSVR_aicc, vcov. = vcovHC(MSVR_aicc, "HC1", adjust = TRUE, sandwich = TRUE))

MSVR_aicc_corr_hac <- coeftest(MSVR_aicc, vcov. = vcovHAC(MSVR_aicc, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MSVR_aicc, MSVR_aicc_corr_hc,
          type = "html", 
          title = "Regresión sin variables relacionadas al uso de tarjetas - Modelo AICc sin efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")
```

###### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MSVR_aicc <- predict(MSVR_aicc, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MSVR_aicc`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Regresión sin variables relacionadas al uso de tarjetas - Modelo AICc sin efectos fijos de tiempo")
pronostico <- c(pred_int_MSVR_aicc)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MSVR_aicc, modelo, pred_int_MSVR_aicc, pronostico, MSVR_aicc_corr_hc, MSVR_aicc_corr_hac)
```

##### MSVR_bic

###### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MSVR_bic <- dredge(modelo_completo_red, rank = "BIC", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MSVR_bic <- get.models(MSVR_bic, 1)[[1]]
MSVR_bic

# Errores robustos
MSVR_bic_corr_hc<-coeftest(MSVR_bic, vcov. = vcovHC(MSVR_bic, "HC1", adjust = TRUE, sandwich = TRUE))

MSVR_bic_corr_hac<-coeftest(MSVR_bic, vcov. = vcovHAC(MSVR_bic, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MSVR_bic, MSVR_bic_corr_hc,
          type = "html", 
          title = "Regresión sin variables relacionadas al uso de tarjetas - Modelo BIC  sin efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")
```

###### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MSVR_bic <- predict(MSVR_bic, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MSVR_bic`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Regresión sin variables relacionadas al uso de tarjetas - Modelo BIC sin efectos fijos de tiempo")
pronostico <- c(pred_int_MSVR_bic)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MSVR_bic, modelo, pred_int_MSVR_bic, pronostico, MSVR_bic_corr_hc, MSVR_bic_corr_hac, modelo_completo_red, formula_con_cont_red)
```

##### MSVR_aicc_et \*\*

###### Estimación

```{r, include=TRUE, results='asis'}
# Ajustar el modelo
modelo_completo_et_red <- lm(formula_con_cont_t_red, data = datos_red)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MSVR_aicc_et <- dredge(modelo_completo_et_red, rank = "AICc", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MSVR_aicc_et <- get.models(MSVR_aicc_et, 1)[[1]]
MSVR_aicc_et

# Errores robustos
MSVR_aicc_et_corr_hc <- coeftest(MSVR_aicc_et, vcov. = vcovHC(MSVR_aicc_et, "HC1", adjust = TRUE, sandwich = TRUE))

MSVR_aicc_et_corr_hac <- coeftest(MSVR_aicc_et, vcov. = vcovHAC(MSVR_aicc_et, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MSVR_aicc_et, MSVR_aicc_et_corr_hc,
          type = "html", 
          title = "Regresión sin variables relacionadas al uso de tarjetas - Modelo AICc con efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")
```

###### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MSVR_aicc_et <- predict(MSVR_aicc_et, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MSVR_aicc_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Regresión sin variables relacionadas al uso de tarjetas - Modelo AICc con efectos fijos de tiempo")
pronostico <- c(pred_int_MSVR_aicc_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MSVR_aicc_et, modelo, pred_int_MSVR_aicc_et, pronostico, MSVR_aicc_et_corr_hc, MSVR_aicc_et_corr_hac)
```

##### MSVR_bic_et \*\*

###### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MSVR_bic_et <- dredge(modelo_completo_et_red, rank = "BIC", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MSVR_bic_et <- get.models(MSVR_bic_et, 1)[[1]]
MSVR_bic_et

# Errores robustos
MSVR_bic_et_corr_hc <- coeftest(MSVR_bic_et, vcov. = vcovHC(MSVR_bic_et, "HC1", adjust = TRUE, sandwich = TRUE))

MSVR_bic_et_corr_hac <- coeftest(MSVR_bic_et, vcov. = vcovHAC(MSVR_bic_et, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MSVR_bic_et, MSVR_bic_et_corr_hc,
          type = "html", 
          title = "Regresión con variables relacionadas al uso de tarjetas - Modelo BIC con efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")
```

###### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_int_MSVR_bic_et <- predict(MSVR_bic_et, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MSVR_bic_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Regresión sin variables relacionadas al uso de tarjetas - Modelo BIC con efectos fijos de tiempo")
pronostico <- c(pred_int_MSVR_bic_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MSVR_bic_et, modelo, pred_int_MSVR_bic_et, pronostico, MSVR_bic_et_corr_hc, MSVR_bic_et_corr_hac, datos_red, modelo_completo_et_red, controles_red, formula_con_cont_t_red, variables_interes_red)
```

#### Enfoque2: modelo con variables instrumentales

::: {style="text-align: justify"}
Para el segundo enfoque, se plantea un modelo de variables instrumentales. En este caso, se instrumenta el valor promedio de la transacción realizada con tarjetas de pago con el Índice de Capital Humano. Este índice incorpora una medición del nivel educativo de un país, que es una característica del cliente que puede estar relacionada con su disposición a adoptar y utilizar tecnologías de pago distintas al efectivo.

Dado que esta característica impacta el lado de la demanda, es posible utilizar el índice como un instrumento para trazar la oferta.
:::

```{r, include=FALSE}
#Datos sin variables relacionadas al uso de tarjetas
v_1 <- c("Pais", "Codigo_Pais", "Año", "Intercambio", "Regulacion_intercambio")
v_2 <- c("Crimen", "Densidad_Pob", "Densidad_Neg", "PIB_pc_PPA", "Europa", "Asia_Este", "Am_lat", "Valor_prom")
v_3 <- c("Crimen", "Densidad_Pob", "Densidad_Neg", "PIB_pc_PPA", "Europa", "Asia_Este", "Am_lat", "ICH")
ef_tiempo_vi <- c("Año_2012", ef_tiempo)
vi <- c(v_1, v_2, "ICH", ef_tiempo_vi)

datos_vi <- datos %>%
  select(all_of(vi)) %>%
  na.omit()

# Fórmula sin efectos de tiempo
formula_vi <- as.formula(
  paste("Intercambio ~ Regulacion_intercambio +", paste(v_2, collapse = " + "), "| Regulacion_intercambio + ", paste(v_3, collapse = " + ") )
)

# Fórmula con efectos de tiempo
formula_vi_et <- as.formula(
  paste("Intercambio ~ Regulacion_intercambio + ", paste(c(v_2,ef_tiempo_vi), collapse = " + "), "| Regulacion_intercambio + ", paste(c(v_3, ef_tiempo_vi), collapse = " + ") )
)

```

##### Modelo sin efectos de tiempo

###### Estimación

```{r, include=TRUE, results='asis'}
#Regresión
modelo_vi <- ivreg(formula_vi,data = datos_vi)
modelo_vi

# Errores robustos
modelo_vi_corr_hc <- coeftest(modelo_vi, vcov. = vcovHC(modelo_vi, "HC1", adjust = TRUE, sandwich = TRUE))

modelo_vi_corr_hac <- coeftest(modelo_vi, vcov. = vcovHAC(modelo_vi, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(modelo_vi, modelo_vi_corr_hc, type = "html", 
          title = "Regresión con variables instrumentales - Sin efectos fijos de tiempo", 
          column.labels = c("Instrumental","Instrumental HC"), model.names = F, 
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")
```

###### Predicción

```{r, include=TRUE}
#Pronósticos para Costa Rica
pred_int_vi <- predict(modelo_vi, newdata = data_CR,
                                   interval = "prediction", level = 0.95,
                                   type = "response")[which(data_CR$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_vi`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Regresión con variables instrumentales - Sin efectos fijos de tiempo")
pronostico <- c(pred_int_vi)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(modelo_vi, modelo, pred_int_vi, pronostico, formula_vi, modelo_vi_corr_hc, modelo_vi_corr_hac)
```

##### Modelo con efectos de tiempo

###### Estimación

```{r, include=TRUE, results='asis'}
#Regresión
modelo_vi_et <- ivreg(formula_vi_et,data = datos_vi)
modelo_vi_et

# Errores robustos
modelo_vi_et_corr_hc <- coeftest(modelo_vi_et, vcov. = vcovHC(modelo_vi_et, "HC1", adjust = TRUE, sandwich = TRUE))

modelo_vi_et_corr_hac <- coeftest(modelo_vi_et, vcov. = vcovHAC(modelo_vi_et, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(modelo_vi_et, modelo_vi_et_corr_hc, type = "html", 
          title = "Regresión con variables instrumentales - Con efectos fijos de tiempo", 
          column.labels = c("Instrumental","Instrumental HC"), model.names = F, 
          dep.var.labels = c("Tasa intercambio", "Tasa intercambio"), 
          dep.var.caption = "Variable dependiente")
```

###### Predicción

```{r, include=TRUE}
#Pronósticos para Costa Rica
pred_int_vi_et <- predict(modelo_vi_et, newdata = data_CR,
                                   interval = "prediction", level = 0.95,
                                   type = "response")[which(data_CR$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_vi_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Regresión con variables instrumentales - Con efectos fijos de tiempo")
pronostico <- c(pred_int_vi_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(modelo_vi_et, modelo, pred_int_vi_et, pronostico, formula_vi_et, modelo_vi_et_corr_hc, modelo_vi_et_corr_hac, datos_vi, ef_tiempo_vi, v_1, v_2, v_3, vi)
```

### Modelo de panel de datos con efectos aleatorios

::: {style="text-align: justify"}
La base de datos recopilada para este estudio contiene para algunos países el seguimiento del comportamiento de las comisiones de adquirencia e intercambio a lo largo del tiempo. Este seguimiento, sin embargo, es limitado. Para el 74% de los países se disponen como máximo de dos observaciones. No obstante, es posible considerar en el presente análisis de sensibilidad un modelo de datos de panel que permita explorar el papel de la heterogeneidad no observable en los resultados.

Para este ejercicio, se utiliza un modelo de efectos aleatorios. Esto ya que Costa Rica no se incluye en la estimación del modelo, sino hasta que se realiza el pronóstico, lo que hace poco práctico un modelo de efectos fijos si el objetivo final del modelo es obtener valores esperados de las comisiones para Costa Rica.
:::

```{r, include=FALSE}
# Crear el identificador de panel
datos_int_con_cont_panel <- datos_int_con_cont %>%
  mutate(id = as.numeric(factor(Codigo_Pais)))

#Declarar estructura de datos de panel
pdata <- pdata.frame(datos_int_con_cont_panel, index = c("id", "Año"))

```

#### MPDEA_aicc

##### Estimación

```{r, include=TRUE}
# Modelo completo
modelo_completo_panel <- lm(formula_con_cont, data = datos_int_con_cont_panel)

# Selección hacia atrás con AICc, fijando Regulacion_intercambio
MPDEA_aicc <- dredge(modelo_completo_panel, fixed = "Regulacion_intercambio", rank = "AICc")

# Extraer el mejor modelo
MPDEA_aicc <- get.models(MPDEA_aicc, 1)[[1]]

# Ver variables seleccionadas
var_panel <- names(coef(MPDEA_aicc))[-1]  # excluye intercepto
var_panel <- var_panel[var_panel != "Regulacion_intercambio"]

# Fórmula para lmer
formula_lmer <- as.formula(paste("Intercambio ~ Regulacion_intercambio +", paste(var_panel, collapse = " + "), " + (1 | Codigo_Pais)" ))

#Regresión con efectos aleatorios y errores agrupados
MPDEA_aicc <- lmer(formula_lmer, data = pdata)

resultado_panel_aicc <- summary(MPDEA_aicc)

summary(MPDEA_aicc)
```

##### Predicción

```{r, include=TRUE}
#Nombre de variables
var_finales <- attr(terms(MPDEA_aicc), "term.labels")

#Pronóstico
Caso_CR <- datos %>% filter(Codigo_Pais=="CRI") %>% dplyr::select(var_finales)
Caso_CR <- cbind(1,Caso_CR) # agregar un 1 para que se multiple con el intercepto general

pred_int_MPDEA_aicc <- sum(resultado_panel_aicc$coefficients*Caso_CR)
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MPDEA_aicc`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación con efectos aleatorios - Modelo AICc sin efectos fijos de tiempo")
pronostico <- c(pred_int_MPDEA_aicc)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MPDEA_aicc, modelo, pred_int_MPDEA_aicc, pronostico, resultado_panel_aicc, formula_lmer, var_finales)
```

#### MPDEA_bic

##### Estimación

```{r, include=TRUE}
# Selección hacia atrás con AICc, fijando Regulacion_intercambio
MPDEA_bic <- dredge(modelo_completo_panel, fixed = "Regulacion_intercambio", rank = "BIC")

# Extraer el mejor modelo
MPDEA_bic <- get.models(MPDEA_bic, 1)[[1]]

# Ver variables seleccionadas
var_panel <- names(coef(MPDEA_bic))[-1]  # excluye intercepto
var_panel <- var_panel[var_panel != "Regulacion_intercambio"]

# Fórmula para lmer
formula_lmer <- as.formula(paste("Intercambio ~ Regulacion_intercambio +", paste(var_panel, collapse = " + "), " + (1 | Codigo_Pais)" ))

#Regresión con efectos aleatorios y errores agrupados
# Modelo de efectos aleatorios
MPDEA_bic <- lmer(formula_lmer, data = pdata)

# Resumen con errores agrupados por id
resultado_panel_bic <- summary(MPDEA_bic)

summary(MPDEA_bic)
```

##### Predicción

```{r, include=TRUE}
#Nombre de variables
var_finales <- attr(terms(MPDEA_bic), "term.labels")

#Pronóstico
Caso_CR <- datos %>% filter(Codigo_Pais=="CRI") %>% dplyr::select(var_finales)
Caso_CR <- cbind(1,Caso_CR) # agregar un 1 para que se multiple con el intercepto general

pred_int_MPDEA_bic <- sum(resultado_panel_bic$coefficients*Caso_CR)
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MPDEA_bic`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación con efectos aleatorios - Modelo BIC sin efectos fijos de tiempo")
pronostico <- c(pred_int_MPDEA_bic)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MPDEA_bic, modelo, pred_int_MPDEA_bic, pronostico, resultado_panel_bic, formula_lmer, var_finales, modelo_completo_panel)
```

#### MPDEA_aicc_et

##### Estimación

```{r, include=TRUE}
# Modelo completo
modelo_completo_panel_et <- lm(formula_con_cont_t, data = datos_int_con_cont_panel)

# Selección hacia atrás con AICc, fijando Regulacion_intercambio
MPDEA_aicc_et <- dredge(modelo_completo_panel_et, fixed = variables_fijas, rank = "AICc")

# Extraer el mejor modelo
MPDEA_aicc_et <- get.models(MPDEA_aicc_et, 1)[[1]]

# Ver variables seleccionadas
coef_original <- coef(MPDEA_aicc_et)
coef_sin_na <- coef_original[!is.na(coef_original)]
var_panel <- names(coef_sin_na)
var_panel <- var_panel[var_panel != "(Intercept)"] # excluye intercepto
var_panel <- var_panel[var_panel != "Regulacion_intercambio"] 

# Fórmula para lmer
formula_lmer <- as.formula(paste("Intercambio ~ Regulacion_intercambio +", paste(var_panel, collapse = " + "), " + (1 | Codigo_Pais)" ))

#Regresión con efectos aleatorios y errores agrupados
MPDEA_aicc_et <- lmer(formula_lmer, data = pdata)

# Resumen con errores agrupados por id
resultado_panel_aicc_et <- summary(MPDEA_aicc_et)

summary(MPDEA_aicc_et)
```

##### Predicción

```{r}
#Nombre de variables
var_finales <- attr(terms(MPDEA_aicc_et), "term.labels")

#Pronóstico
Caso_CR <- datos %>% filter(Codigo_Pais=="CRI") %>% dplyr::select(var_finales)
Caso_CR <- cbind(1,Caso_CR) # agregar un 1 para que se multiple con el intercepto general

pred_int_MPDEA_aicc_et <- sum(resultado_panel_aicc_et$coefficients*Caso_CR)
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MPDEA_aicc_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación con efectos aleatorios - Modelo AICc con efectos fijos de tiempo")
pronostico <- c(pred_int_MPDEA_aicc_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MPDEA_aicc_et, modelo, pred_int_MPDEA_aicc_et, pronostico, resultado_panel_aicc_et, formula_lmer, var_finales, coef_original, coef_sin_na, var_panel)
```

#### MPDEA_bic_et

##### Estimación

```{r, include=TRUE}
# Selección hacia atrás con AICc, fijando Regulacion_intercambio
MPDEA_bic_et <- dredge(modelo_completo_panel_et, fixed = variables_fijas, rank = "BIC")

# Extraer el mejor modelo
MPDEA_bic_et <- get.models(MPDEA_bic_et, 1)[[1]]

# Ver variables seleccionadas
coef_original <- coef(MPDEA_bic_et)
coef_sin_na <- coef_original[!is.na(coef_original)]
var_panel <- names(coef_sin_na)
var_panel <- var_panel[var_panel != "(Intercept)"] # excluye intercepto
var_panel <- var_panel[var_panel != "Regulacion_intercambio"] 

# Fórmula para lmer
formula_lmer <- as.formula(paste("Intercambio ~ Regulacion_intercambio +", paste(var_panel, collapse = " + "), " + (1 | Codigo_Pais)" ))

#Regresión con efectos aleatorios y errores agrupados
MPDEA_bic_et <- lmer(formula_lmer, data = pdata)

# Resumen con errores agrupados por id
resultado_panel_bic_et <- summary(MPDEA_bic_et)

summary(MPDEA_bic_et)
```

##### Predicción

```{r, include=TRUE}
#Nombre de variables
var_finales <- attr(terms(MPDEA_bic_et), "term.labels")

#Pronóstico
Caso_CR <- datos %>% filter(Codigo_Pais=="CRI") %>% dplyr::select(var_finales)
Caso_CR <- cbind(1,Caso_CR) # agregar un 1 para que se multiple con el intercepto general

pred_int_MPDEA_bic_et <- sum(resultado_panel_bic_et$coefficients*Caso_CR)
```

::: {style="text-align: justify"}
Predicción tasa de intercambio para Costa Rica: **`r pred_int_MPDEA_bic_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación con efectos aleatorios - Modelo BIC con efectos fijos de tiempo")
pronostico <- c(pred_int_MPDEA_bic_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MPDEA_bic_et, modelo, pred_int_MPDEA_bic_et, pronostico, resultado_panel_bic_et, formula_lmer, var_finales, coef_original, coef_sin_na, var_panel, datos_int_con_cont_panel, modelo_completo_panel_et, pdata, Caso_CR)
```

## Estimación de las comisiones de intercambio

::: {style="text-align: justify"}
Posterior a la etapa de modelación, se realiza un pronóstico para conocer la tasa de intercambio esperada para Costa Rica. A continuación, se calcula un promedio para los valores esperados de cada modelo estimado.
:::

```{r, include=TRUE}
Pronosticos$Tipo_modelo <- Pronosticos$Modelo

Pronosticos <- Pronosticos %>%
  separate(Tipo_modelo, into = c("Tipo_modelo", "Especificacion"), sep = "-")


Pronosticos_resumen <- Pronosticos[-1, ] %>%
  group_by(Tipo_modelo) %>%
  summarise(promedio = mean(Pronostico, na.rm = TRUE))

promedio_intercambio_CR <- Pronosticos_resumen %>%
                            filter(Tipo_modelo == "Estimación principal ") %>%
                            pull(promedio)

print(Pronosticos_resumen)

Pronosticos_resumen_intercambio <- Pronosticos_resumen

```

::: {style="text-align: justify"}
En promedio, la tasa de intercambio estimada para Costa Rica es: **`r promedio_intercambio_CR`**.
:::

# Comisión de adquirencia

## Estimación principal

::: {style="text-align: justify"}
Se estima una regresión lineal múltiple (MRL). A partir de los coeficientes estimados y las características observables para Costa Rica, se calcula la comisión de adquirencia esperada para el país. Se calculan 5 modelos:

* MRL_sc: estima una regresión lineal sin variables de control; únicamente se incluye la variable dicotómica que indica si el país cuenta con regulación.
* MRL_aicc: corresponde a los resultados obtenidos con el criterio de información de Akaike (AICc) con una corrección para tamaños de muestra pequeñas.
* MRL_bic: incluye las estimaciones al controlar por las variables seleccionadas mediante el criterio de información bayesiano (BIC).
* MRL_aicc_et: incluyen efectos fijos de tiempo en las regresiones, para el modelo MRL_aicc.
* MRL_bic_et: incluyen efectos fijos de tiempo en las regresiones, para el modelo MRL_bic.
:::

```{r, include=FALSE}
#Países con información sobre intercambio en la muestra (distintos de C.R.)

n_paises_adq <- datos[!is.na(datos$Adquirencia), ] %>%
  filter(Codigo_Pais != "CRI") %>%
  distinct(Codigo_Pais) %>%
  nrow()

```

::: {style="text-align: justify"}
El total de países con información sobre adquirencia es **`r n_paises_adq`**.
:::

### Definiciones

``` {r, include=TRUE}
# Variables de control
controles <- c("Cajero_Autom", "Crimen", "Cuenta", "Credito", "Debito", "Densidad_Pob", "Densidad_Neg", "PIB_pc_PPA", "ICH", "Pago_digital", "Valor_prom", "Asia_Este", "Europa", "Am_lat")

# Efectos fijos de tiempo
ef_tiempo <- c("Año_2013", "Año_2014", "Año_2015", "Año_2016", "Año_2017", "Año_2018", "Año_2019", "Año_2020", "Año_2021", "Año_2022", "Año_2023", "Año_2024", "Año_2025")

# Fórmula que incluye las variables de control
formula_con_cont <- as.formula(
  paste("Adquirencia ~ Regulacion_intercambio +", paste(controles, collapse = " + "))
)

# Fórmula que incluye las variables de control y los efectos fijos de tiempo
formula_con_cont_t <- as.formula(
  paste("Adquirencia ~ Regulacion_intercambio +", paste(controles, collapse = " + ")," + ", paste(ef_tiempo, collapse = " + "))
)

# Variables fijas para los modelos con efectos de tiempo
variables_fijas <- c("Regulacion_intercambio", paste0("Año_", 2010:2025))

# Definir el conjuntos de datos para los modelos con variables de control
variables_constantes <- c("Pais", "Codigo_Pais", "Año", "Adquirencia", "Regulacion_intercambio")
variables_interes <- c(variables_constantes, controles, ef_tiempo)

datos_adq_con_cont <- datos %>%
  select(all_of(variables_interes)) %>%
  na.omit()

```

::: {style="text-align: justify"}
Variables de control: `r paste(capture.output(controles), collapse = "<br>")`.

Efectos fijos de tiempo: `r paste(capture.output(ef_tiempo), collapse = "<br>")`.
:::

### MRL_sc

::: {style="text-align: justify"}
Regresión lineal sin variables de control.
:::

#### Estimación

```{r, include=TRUE, results='asis'}
#Remover NA's de la columna "Intercambio" y dejar los demás.
datos_adq_sin_cont <- datos %>% 
  filter(!is.na(Adquirencia))

MRL_sc <- lm(Adquirencia ~ Regulacion_intercambio, datos_adq_sin_cont)

# Errores robustos
MRL_sc_corr_hc <- coeftest(MRL_sc, vcov. = vcovHC(MRL_sc, "HC1", adjust = TRUE, sandwich = TRUE))

MRL_sc_corr_hac <- coeftest(MRL_sc, vcov. = vcovHAC(MRL_sc, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MRL_sc, MRL_sc_corr_hc, type = "html", 
          title = "Estimación principal - Modelo sin variables de control", 
          column.labels = c("MCO","MCO HC"), model.names = F, 
          dep.var.labels = c("Tasa adquirencia", "Tasa adquirencia"), 
          dep.var.caption = "Variable dependiente")
```

#### Pruebas estadísticas

```{r, include=TRUE}
#Prueba de normalidad
residuos_MRL_sc <- resid(MRL_sc) #Se guardan los residuos
jarque.test(residuos_MRL_sc) #Se computa la prueba con la función Nula: distribución normal residuos
car::qqPlot(MRL_sc)

#Test de White sobre el modelo original. Hipótesis nula: homocedasticidad
white(MRL_sc, interactions = TRUE)

```

#### Predicción

```{r, include=TRUE}
#Pronóstico para CR
data_CR <- datos
data_CR[ which(datos$Codigo_Pais=="CRI"), which(names(datos)=="Regulacion_intercambio") ]<-1
pred_adq_MRL_sc <- predict(MRL_sc, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MRL_sc`**.
:::

```{r, include=FALSE}
#Crea el data frame para guardar la predicción
modelo <- c("Estimación principal - Modelo sin variables de control")
pronostico <- c(pred_adq_MRL_sc)  
Pronosticos <- data.frame(Modelo = modelo, Pronostico = pronostico)
```


```{r, include=FALSE}
rm(datos_adq_sin_cont, MRL_sc, MRL_sc_corr_hac, MRL_sc_corr_hc, residuos_MRL_sc, pred_adq_MRL_sc, modelo, pronostico)
```

###  MRL_aicc

::: {style="text-align: justify"}
Regresión lineal con variables de control, sin efectos fijos de tiempo. Utiliza el criterio de información de Akaike (AICc) con una corrección para tamaños de muestra pequeñas.
:::

#### Estimación

```{r, include=TRUE, results='asis'}
# Ajustar el modelo
modelo_completo <- lm(formula_con_cont, data = datos_adq_con_cont)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRL_aicc <- dredge(modelo_completo, rank = "AICc", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MRL_aicc <- get.models(MRL_aicc, 3)[[1]] #Ranking=3
MRL_aicc

# Errores robustos
MRL_aicc_corr_hc <- coeftest(MRL_aicc, vcov. = vcovHC(MRL_aicc, "HC1", adjust = TRUE, sandwich = TRUE))

MRL_aicc_corr_hac <- coeftest(MRL_aicc, vcov. = vcovHAC(MRL_aicc, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MRL_aicc, MRL_aicc_corr_hc,
          type = "html", 
          title = "Estimación principal - Modelo AICc sin efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa adquirencia", "Tasa adquirencia"), 
          dep.var.caption = "Variable dependiente")
```

#### Pruebas estadísticas

```{r, include=TRUE}
#Prueba de normalidad

residuos_MRL_aicc<-resid(MRL_aicc) #Se guardan los residuos
jarque.test(residuos_MRL_aicc) #Se computa la prueba con la función Nula: distribución normal residuos
car::qqPlot(MRL_aicc)

#Se obtiene el resultado del índice de condición de multicolinealidad
#create vector of VIF values
print("El VIF del modelo es:")
vif(MRL_aicc)

#Se ejecuta el test de White sobre el modelo original Nula: homocedasticidad
#white_lm(modelo_int_con_cont_aicc, interactions = TRUE)
white(MRL_aicc, interactions = TRUE)
plot(MRL_aicc$fit,MRL_aicc$res,xlab="Valores ajustados",ylab="Residuales")

#Prueba de especificación lineal RESET Nula: especificación correcta lineal
resettest(MRL_aicc, power = 2:3, type = c("fitted"), data = datos_adq_con_cont, vcov = vcovHC) #Corregido por heterocedasticidad)

```

#### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MRL_aicc <- predict(MRL_aicc, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MRL_aicc`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación principal - Modelo AICc sin efectos fijos de tiempo")
pronostico <- c(pred_adq_MRL_aicc)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRL_aicc, modelo, MRL_aicc_corr_hac, MRL_aicc_corr_hc, pred_adq_MRL_aicc, pronostico, residuos_MRL_aicc)
```

### MRL_bic

::: {style="text-align: justify"}
Regresión lineal con variables de control, sin efectos fijos de tiempo. Utiliza el criterio de información bayesiana (BIC).
:::

#### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRL_bic <- dredge(modelo_completo, rank = "BIC", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MRL_bic <- get.models(MRL_bic, 20)[[1]] #Ranking=20
MRL_bic

# Errores robustos
MRL_bic_corr_hc <- coeftest(MRL_bic, vcov. = vcovHC(MRL_bic, "HC1", adjust = TRUE, sandwich = TRUE))

MRL_bic_corr_hac <- coeftest(MRL_bic, vcov. = vcovHAC(MRL_bic, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MRL_bic, MRL_bic_corr_hc, type = "html",
          title = "Estimación principal - Modelo BIC sin efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F, 
          dep.var.labels = c("Tasa adquirencia", "Tasa adquirencia"), 
          dep.var.caption = "Variable dependiente")


```

#### Pruebas estadísticas

```{r, include=TRUE}
#Prueba de normalidad
#Se guardan los residuos
residuos_MRL_bic<-resid(MRL_bic)
#Se computa la prueba con la función Nula: distribución normal residuos
jarque.test(residuos_MRL_bic)

#Se obtiene el resultado del índice de condición de multicolinealidad
#create vector of VIF values
print("El VIF del modelo es:")
vif(MRL_bic)

#Se ejecuta el test de White sobre el modelo original Nula: homocedasticidad
#white_lm(modelo_int_con_cont_bic, interactions = TRUE)
white(MRL_bic, interactions = TRUE)
#Prueba de especificación lineal RESET Nula: especificación correcta lineal
resettest(MRL_bic, power = 2:3, type = c("fitted"), data = datos_adq_con_cont, vcov = vcovHC)#Corregido por heterocedasticidad)

```

#### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MRL_bic<-predict(MRL_bic, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MRL_bic`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación principal - Modelo BIC sin efectos fijos de tiempo")
pronostico <- c(pred_adq_MRL_bic)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRL_bic, modelo, MRL_bic_corr_hac, MRL_bic_corr_hc, pred_adq_MRL_bic, pronostico, residuos_MRL_bic)
```

###  MRL_aicc_et

::: {style="text-align: justify"}
Regresión lineal con variables de control, con efectos fijos de tiempo. Utiliza el criterio de información de Akaike (AICc) con una corrección para tamaños de muestra pequeñas.
:::

#### Estimación

```{r, include=TRUE, results='asis'}
# Ajustar el modelo
modelo_completo_et <- lm(formula_con_cont_t, data = datos_adq_con_cont)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRL_aicc_et <- dredge(modelo_completo_et, rank = "AICc", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MRL_aicc_et <- get.models(MRL_aicc_et, 11)[[1]] #Ranking = 11
MRL_aicc_et

# Errores robustos
MRL_aicc_et_corr_hc <- coeftest(MRL_aicc_et, vcov. = vcovHC(MRL_aicc_et, "HC1", adjust = TRUE, sandwich = TRUE))

MRL_aicc_et_corr_hac <- coeftest(MRL_aicc_et, vcov. = vcovHAC(MRL_aicc_et, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MRL_aicc_et, MRL_aicc_et_corr_hc,
          type = "html", 
          title = "Estimación principal - Modelo AICc con efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa adquirencia", "Tasa adquirencia"), 
          dep.var.caption = "Variable dependiente")
```

#### Pruebas estadísticas

```{r, include=TRUE}

#Prueba de normalidad
#Se guardan los residuos
residuos_MRL_aicc_et <- resid(MRL_aicc_et)
#Se computa la prueba con la función Nula: distribución normal residuos
jarque.test(residuos_MRL_aicc_et)
car::qqPlot(residuos_MRL_aicc_et)

#Se ejecuta el test de White sobre el modelo original Nula: homocedasticidad
#white_lm(modelo_int_con_cont_ef_aicc, interactions = TRUE)
white(MRL_aicc_et, interactions = TRUE)

```

#### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MRL_aicc_et<-predict(MRL_aicc_et, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MRL_aicc_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación principal - Modelo AICc con efectos fijos de tiempo")
pronostico <- c(pred_adq_MRL_aicc_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRL_aicc_et, modelo, MRL_aicc_et_corr_hac, MRL_aicc_et_corr_hc, pred_adq_MRL_aicc_et, pronostico, residuos_MRL_aicc_et)
```

### MRL_bic_et

::: {style="text-align: justify"}
Regresión lineal con variables de control, con efectos fijos de tiempo. Utiliza el criterio de información bayesiana (BIC).
:::

#### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRL_bic_et <- dredge(modelo_completo_et, rank = "BIC", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MRL_bic_et <- get.models(MRL_bic_et, 10)[[1]] #Ranking=10
MRL_bic_et

# Errores robustos
MRL_bic_et_corr_hc <- coeftest(MRL_bic_et, vcov. = vcovHC(MRL_bic_et, "HC1", adjust = TRUE, sandwich = TRUE))

MRL_bic_et_corr_hac <- coeftest(MRL_bic_et, vcov. = vcovHAC(MRL_bic_et, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MRL_bic_et, MRL_bic_et_corr_hc, type = "html",
          title = "Estimación principal - Modelo BIC con efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F, 
          dep.var.labels = c("Tasa adquirencia", "Tasa adquirencia"), 
          dep.var.caption = "Variable dependiente")

```

#### Pruebas estadísticas

```{r, include=TRUE}
#Prueba de normalidad
#Se guardan los residuos
residuos_MRL_bic_et <- resid(MRL_bic_et)
#Se computa la prueba con la función Nula: distribución normal residuos
jarque.test(residuos_MRL_bic_et)
white(MRL_bic_et, interactions = TRUE)
```

#### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MRL_bic_et <- predict(MRL_bic_et, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MRL_bic_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación principal - Modelo BIC con efectos fijos de tiempo")
pronostico <- c(pred_adq_MRL_bic_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRL_bic_et, modelo, MRL_bic_et_corr_hac, MRL_bic_et_corr_hc, pred_adq_MRL_bic_et, pronostico, residuos_MRL_bic_et)
```

## Pruebas de robustez

::: {style="text-align: justify"}
A continuación, se evalúa cómo cambian los resultados anteriores al estimar las comisiones con modelos alternativos.
:::

###	Modelos de regresión truncados

::: {style="text-align: justify"}
Mientras que el modelo de regresión lineal no impone ninguna restricción a las comisiones, es de esperar que las comisiones sean necesariamente mayores a cero como mínimo, o mayores a los costos asociados en particular. Para tratar estos casos, se pueden utilizar modelos de variable dependiente truncada. 

Se estiman estos modelos bajo dos supuestos sobre el punto de truncamiento. Un primer punto de partida es que el punto de truncamiento es en cero. El segundo supuesto es que el punto de truncamiento sucede en el mínimo de la comisión de intercambio observada en la muestra. 

Además, se estima cada modelo con las variables seleccionadas por los criterios de información de AICc y BIC. 
:::

#### MRT_aicc_0

##### Estimación

```{r, include=TRUE}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRT_aicc_0 <- dredge(modelo_completo, rank = "AICc", fixed = "Regulación_intercambio")

# Extraer el mejor modelo según AICc
MRT_aicc_0 <- names(coef(get.models(MRT_aicc_0, 3)[[1]])) #Ranking=3
MRT_aicc_0 <- MRT_aicc_0[MRT_aicc_0 != "(Intercept)"]  # Excluir intercepto

# Ajustar modelo truncado con límite inferior 0
formula_trunc <- as.formula(paste("Adquirencia ~", paste(MRT_aicc_0, collapse = " + ")))
MRT_aicc_0 <- truncreg(formula_trunc, data = datos_adq_con_cont, point = 0, direction = "left")

summary(MRT_aicc_0)
```

##### Predicción

```{r, include=TRUE, warning=FALSE}
#Pronóstico para CR
pred_adq_MRT_aicc_0 <- predict(MRT_aicc_0, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MRT_aicc_0`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación truncada en 0 - Modelo AICc sin efectos fijos de tiempo")
pronostico <- c(pred_adq_MRT_aicc_0)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRT_aicc_0, modelo, pred_adq_MRT_aicc_0, pronostico, lastFuncGrad, lastFuncParam)
```

#### MRT_bic_0 **

##### Estimación

```{r, include=TRUE}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRT_bic_0 <- dredge(modelo_completo, rank = "BIC", fixed = "Regulación_intercambio")

# Extraer el mejor modelo según AICc
MRT_bic_0 <- names(coef(get.models(MRT_bic_0, 20)[[1]])) #Ranking=20
MRT_bic_0 <- MRT_bic_0[MRT_bic_0 != "(Intercept)"]  # Excluir intercepto

# Ajustar modelo truncado con límite inferior 0
formula_trunc <- as.formula(paste("Adquirencia ~", paste(MRT_bic_0, collapse = " + ")))
MRT_bic_0 <- truncreg(formula_trunc, data = datos_adq_con_cont, point = 0, direction = "left")

summary(MRT_bic_0)
```

##### Predicción

```{r, include=TRUE, warning=FALSE}
#Pronóstico para CR
pred_adq_MRT_bic_0<-predict(MRT_bic_0, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MRT_bic_0`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación truncada en 0 - Modelo BIC sin efectos fijos de tiempo")
pronostico <- c(pred_adq_MRT_bic_0)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRT_bic_0, modelo, pred_adq_MRT_bic_0, pronostico, lastFuncGrad, lastFuncParam)
```

#### MRT_aicc_et_0

##### Estimación

```{r, include=TRUE, warning=FALSE}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRT_aicc_et_0 <- dredge(modelo_completo_et, rank = "AICc", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MRT_aicc_et_0 <- get.models(MRT_aicc_et_0, 11)[[1]] #Rankig=11

# Obtener coeficientes
coeficientes <- coef(MRT_aicc_et_0)

# Filtrar nombres de coeficientes que no sean NA
coeficientes_validos <- names(coeficientes[!is.na(coeficientes)])

# Excluir el intercepto
variables_truncadas <- coeficientes_validos[coeficientes_validos != "(Intercept)"]

# Construir fórmula para modelo truncado
formula_trunc <- as.formula(paste("Adquirencia ~", paste(variables_truncadas, collapse = " + ")))

# Ajustar modelo truncado con límite inferior 0
MRT_aicc_et_0 <- truncreg(formula_trunc, data = datos_adq_con_cont, point = 0, direction = "left")

summary(MRT_aicc_et_0)
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MRT_aicc_et_0<-predict(MRT_aicc_et_0, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MRT_aicc_et_0`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación truncada en 0 - Modelo AICc con efectos fijos de tiempo")
pronostico <- c(pred_adq_MRT_aicc_et_0)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRT_aicc_et_0, modelo, pred_adq_MRT_aicc_et_0, pronostico, lastFuncGrad, lastFuncParam, coeficientes, coeficientes_validos)
```

#### MRT_bic_et_0

##### Estimación

```{r, include=TRUE, warning=FALSE}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRT_bic_et_0 <- dredge(modelo_completo_et, rank = "BIC", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MRT_bic_et_0 <- get.models(MRT_bic_et_0, 10)[[1]] #Ranking=10

# Obtener coeficientes
coeficientes <- coef(MRT_bic_et_0)

# Filtrar nombres de coeficientes que no sean NA
coeficientes_validos <- names(coeficientes[!is.na(coeficientes)])

# Excluir el intercepto
variables_truncadas <- coeficientes_validos[coeficientes_validos != "(Intercept)"]

# Construir fórmula para modelo truncado
formula_trunc <- as.formula(paste("Adquirencia ~", paste(variables_truncadas, collapse = " + ")))

# Ajustar modelo truncado con límite inferior 0
MRT_bic_et_0 <- truncreg(formula_trunc, data = datos_adq_con_cont, point = 0, direction = "left")

summary(MRT_bic_et_0)
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MRT_bic_et_0<-predict(MRT_bic_et_0, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MRT_bic_et_0`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación truncada en 0 - Modelo BIC con efectos fijos de tiempo")
pronostico <- c(pred_adq_MRT_bic_et_0)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRT_bic_et_0, modelo, pred_adq_MRT_bic_et_0, pronostico, lastFuncGrad, lastFuncParam, coeficientes, coeficientes_validos)
```

#### MRT_aicc_033

##### Estimación

```{r, include=TRUE}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRT_aicc_033 <- dredge(modelo_completo, rank = "AICc", fixed = "Regulación_intercambio")

# Extraer el mejor modelo según AICc
MRT_aicc_033 <- names(coef(get.models(MRT_aicc_033, 3)[[1]])) #Ranking=3
MRT_aicc_033 <- MRT_aicc_033[MRT_aicc_033 != "(Intercept)"]  # Excluir intercepto

# Ajustar modelo truncado con límite inferior 0.33
formula_trunc <- as.formula(paste("Adquirencia ~", paste(MRT_aicc_033, collapse = " + ")))
MRT_aicc_033 <- truncreg(formula_trunc, data = datos_adq_con_cont, point = 0.33, direction = "left")

summary(MRT_aicc_033)
```

##### Predicción

```{r, include=TRUE, warning=FALSE}
#Pronóstico para CR
pred_adq_MRT_aicc_033 <- predict(MRT_aicc_033, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MRT_aicc_033`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación truncada en 0.33 - Modelo AICc sin efectos fijos de tiempo")
pronostico <- c(pred_adq_MRT_aicc_033)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRT_aicc_033, modelo, pred_adq_MRT_aicc_033, pronostico, lastFuncGrad, lastFuncParam)
```

#### MRT_bic_033

##### Estimación

```{r, include=TRUE}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRT_bic_033 <- dredge(modelo_completo, rank = "BIC", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MRT_bic_033 <- names(coef(get.models(MRT_bic_033, 20)[[1]])) #Ranking=20
MRT_bic_033 <- MRT_bic_033[MRT_bic_033 !="(Intercept)"]  # Excluir intercepto

# Ajustar modelo truncado con límite inferior 0.33
formula_trunc <- as.formula(paste("Adquirencia ~", paste(MRT_bic_033, collapse = " + ")))
MRT_bic_033 <- truncreg(formula_trunc, data = datos_adq_con_cont, point = 0.33, direction = "left")

summary(MRT_bic_033)
```

##### Predicción

```{r, include=TRUE, warning=FALSE}
#Pronóstico para CR
pred_adq_MRT_bic_033 <- predict(MRT_bic_033, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MRT_bic_033`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación truncada en 0.33 - Modelo BIC sin efectos fijos de tiempo")
pronostico <- c(pred_adq_MRT_bic_033)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRT_bic_033, modelo, pred_adq_MRT_bic_033, pronostico, lastFuncGrad, lastFuncParam)
```

#### MRT_aicc_et_033

##### Estimación

```{r, include=TRUE, warning=FALSE}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRT_aicc_et_033 <- dredge(modelo_completo_et, rank = "AICc", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MRT_aicc_et_033 <- get.models(MRT_aicc_et_033, 11)[[1]] #Ranking=11

# Obtener coeficientes
coeficientes <- coef(MRT_aicc_et_033)

# Filtrar nombres de coeficientes que no sean NA
coeficientes_validos <- names(coeficientes[!is.na(coeficientes)])

# Excluir el intercepto
variables_truncadas <- coeficientes_validos[coeficientes_validos != "(Intercept)"]

# Construir fórmula para modelo truncado
formula_trunc <- as.formula(paste("Adquirencia ~", paste(variables_truncadas, collapse = " + ")))

# Ajustar modelo truncado con límite inferior 0.33
MRT_aicc_et_033 <- truncreg(formula_trunc, data = datos_adq_con_cont, point = 0.33, direction = "left")

summary(MRT_aicc_et_033)
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MRT_aicc_et_033 <- predict(MRT_aicc_et_033, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MRT_aicc_et_033`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación truncada en 0.33 - Modelo AICc con efectos fijos de tiempo")
pronostico <- c(pred_adq_MRT_aicc_et_033)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRT_aicc_et_033, modelo, pred_adq_MRT_aicc_et_033, pronostico, lastFuncGrad, lastFuncParam, coeficientes, coeficientes_validos)
```

#### MRT_bic_et_033

##### Estimación

```{r, include=TRUE, warning=FALSE}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRT_bic_et_033 <- dredge(modelo_completo_et, rank = "BIC", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MRT_bic_et_033 <- get.models(MRT_bic_et_033, 10)[[1]] #Ranking=10

# Obtener coeficientes
coeficientes <- coef(MRT_bic_et_033)

# Filtrar nombres de coeficientes que no sean NA
coeficientes_validos <- names(coeficientes[!is.na(coeficientes)])

# Excluir el intercepto
variables_truncadas <- coeficientes_validos[coeficientes_validos != "(Intercept)"]

# Construir fórmula para modelo truncado
formula_trunc <- as.formula(paste("Adquirencia ~", paste(variables_truncadas, collapse = " + ")))

# Ajustar modelo truncado con límite inferior 0.33
MRT_bic_et_033 <- truncreg(formula_trunc, data = datos_adq_con_cont, point = 0.33, direction = "left")

summary(MRT_bic_et_033)
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MRT_bic_et_033 <- predict(MRT_bic_et_033, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MRT_bic_et_033`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación truncada en 0.33 - Modelo BIC con efectos fijos de tiempo")
pronostico <- c(pred_adq_MRT_bic_et_033)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRT_bic_et_033, modelo, pred_adq_MRT_bic_et_033, pronostico, lastFuncGrad, lastFuncParam, coeficientes, coeficientes_validos, variables_truncadas, formula_trunc)
```

###	Modelos de regresión con transformación logarítmica de la variable dependiente

::: {style="text-align: justify"}
Al no existir ningún tipo de restricción en la estimación de los coeficientes, cabe la posibilidad de que un modelo de regresión lineal genere predicciones negativas. Sin embargo, ninguna de las predicciones obtenidas para el intercambio o adquirencia son negativas.
 
A pesar de no significar un potencial problema, una manera de garantizar que las predicciones sean positivas es recurrir a una transformación logarítmica de la variable dependiente.
:::

```{r, include=FALSE}
## Transformación a logaritmos de las comisiones
datos_adq_con_cont_log <- datos_adq_con_cont
datos_adq_con_cont_log$Adquirencia <- log(datos_adq_con_cont_log$Adquirencia)
```

#### MTL_aicc

##### Estimación

```{r, include=TRUE}
# Ajustar el modelo
modelo_completo_log <- lm(formula_con_cont, data = datos_adq_con_cont_log)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MTL_aicc <- dredge(modelo_completo_log, rank = "AICc", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MTL_aicc <- get.models(MTL_aicc, 1)[[1]]
MTL_aicc

# Errores robustos
MTL_aicc_corr_hc <- coeftest(MTL_aicc, vcov. = vcovHC(MTL_aicc, "HC1", adjust = TRUE, sandwich = TRUE))

MTL_aicc_corr_hac <- coeftest(MTL_aicc, vcov. = vcovHAC(MTL_aicc, adjust = TRUE, sandwich = TRUE))

summary(MTL_aicc)
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MTL_aicc <- predict(MTL_aicc, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MTL_aicc`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación con transformación logarítmica - Modelo AICc sin efectos fijos de tiempo")
pronostico <- c(pred_adq_MTL_aicc)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MTL_aicc, modelo, pred_adq_MTL_aicc, pronostico, MTL_aicc_corr_hc, MTL_aicc_corr_hac)
```

#### MTL_bic

##### Estimación

```{r, include=TRUE}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MTL_bic <- dredge(modelo_completo_log, rank = "BIC", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MTL_bic <- get.models(MTL_bic, 3)[[1]] #Ranking=3
MTL_bic

# Errores robustos
MTL_bic_corr_hc <- coeftest(MTL_bic, vcov. = vcovHC(MTL_bic, "HC1", adjust = TRUE, sandwich = TRUE))

MTL_bic_corr_hac <- coeftest(MTL_bic, vcov. = vcovHAC(MTL_bic, adjust = TRUE, sandwich = TRUE))

summary(MTL_bic)
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MTL_bic <- predict(MTL_bic, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MTL_bic`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación con transformación logarítmica - Modelo BIC sin efectos fijos de tiempo")
pronostico <- c(pred_adq_MTL_bic)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MTL_bic, modelo, pred_adq_MTL_bic, pronostico, MTL_bic_corr_hc, MTL_bic_corr_hac, modelo_completo_log)
```

#### MTL_aicc_et **

##### Estimación

```{r, include=TRUE}
# Ajustar el modelo
#modelo_completo_log_et <- lm(formula_con_cont_t, data = datos_adq_con_cont_log)
modelo_completo_log_et <- lm(Adquirencia ~ Regulacion_intercambio + Pago_digital +                                     Valor_prom + Asia_Este + Europa + Año_2014 + Año_2015 +  
                             Año_2016 + Año_2017 + Año_2018 + Año_2019 + Año_2020 + 
                             Año_2021 + Año_2022 + Año_2023 + Año_2024 + Año_2025, 
                             data = datos_adq_con_cont_log)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MTL_aicc_et <- dredge(modelo_completo_log_et, rank = "AICc", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MTL_aicc_et <- get.models(MTL_aicc_et, 1)[[1]]
MTL_aicc_et

# Errores robustos
MTL_aicc_et_corr_hc <- coeftest(MTL_aicc_et, vcov. = vcovHC(MTL_aicc_et, "HC1", adjust = TRUE, sandwich = TRUE))

MTL_aicc_et_corr_hac <- coeftest(MTL_aicc_et, vcov. = vcovHAC(MTL_aicc_et, adjust = TRUE, sandwich = TRUE))

summary(MTL_aicc_et)
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MTL_aicc_et <- predict(MTL_aicc_et, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MTL_aicc_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación con transformación logarítmica - Modelo AICc con efectos fijos de tiempo")
pronostico <- c(pred_adq_MTL_aicc_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MTL_aicc_et, modelo, pred_adq_MTL_aicc_et, pronostico, MTL_aicc_et_corr_hc, MTL_aicc_et_corr_hac)
```

#### MTL_bic_et **

##### Estimación

```{r, include=TRUE}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MTL_bic_et <- dredge(modelo_completo_log_et, rank = "BIC", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MTL_bic_et <- get.models(MTL_bic_et, 1)[[1]]
MTL_bic_et

# Errores robustos
MTL_bic_et_corr_hc <- coeftest(MTL_bic_et, vcov. = vcovHC(MTL_bic_et, "HC1", adjust = TRUE, sandwich = TRUE))

MTL_bic_et_corr_hac <- coeftest(MTL_bic_et, vcov. = vcovHAC(MTL_bic_et, adjust = TRUE, sandwich = TRUE))

summary(MTL_bic_et)
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MTL_bic_et <- predict(MTL_bic_et, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MTL_bic_et**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación con transformación logarítmica - Modelo BIC con efectos fijos de tiempo")
pronostico <- c(pred_adq_MTL_bic_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MTL_bic_et, modelo, pred_adq_MTL_bic_et, pronostico, MTL_bic_et_corr_hc, MTL_bic_et_corr_hac, datos_adq_con_cont_log, modelo_completo_log_et)
```

### Modelo de regresión con efectos no lineales

::: {style="text-align: justify"}
Para permitir la presencia de efectos no lineales en la regresión, se estima una especificación donde se incluyen como potenciales candidatos en la selección del modelo los términos cuadráticos de las variables de control.
:::

```{r, include=FALSE}
## Arreglo de los datos
data_CR_enl <- data_CR
datos_adq_con_cont_enl <- datos_adq_con_cont

controles_numericos <- setdiff(controles, c("Europa", "Asia_Este", "Am_lat"))

for (var in controles_numericos) {
  datos_adq_con_cont_enl[[paste0(var, "_sqrt")]] <- datos_adq_con_cont_enl[[var]]^2
}

for (var in controles_numericos) {
  data_CR_enl[[paste0(var, "_sqrt")]] <- data_CR_enl[[var]]^2
}

controles_sqrt <- c(controles, paste0(controles_numericos, "_sqrt"))


# Fórmula para incluir las variables de control
formula_con_cont_enl <- as.formula(
  paste("Adquirencia ~ Regulacion_intercambio +", paste(controles_sqrt, collapse = " + "))
)

# Fórmula para incluir las variables de control con efectos fijos de tiempo
formula_con_cont_t_enl <- as.formula(
  paste("Adquirencia ~ Regulacion_intercambio +", paste(controles_sqrt, collapse = " + ")," + ", paste(ef_tiempo, collapse = " + "))
)

```

#### MENL_aicc

##### Estimación

::: {style="text-align: justify"}
Se proporciona el código completo para su ejecución en R. No obstante, debido al tiempo excesivo que toma el procesamiento automático, se han incorporado manualmente las variables que optimizan el modelo. Estas variables fueron identificadas previamente mediante un proceso de selección realizado en Stata, con el fin de reducir significativamente el tiempo de ejecución en R.
:::

```{r, include=TRUE}
# Ajustar el modelo (automático R)
#modelo_completo_enl <- lm(formula_con_cont_enl, data = datos_adq_con_cont_enl)

# Ajustar el modelo (manual, selección previa de variables en Stata)
modelo_completo_enl <- lm(Adquirencia ~ Regulacion_intercambio + 
                                     Crimen + Cuenta + Credito + PIB_pc_PPA + 
                                     ICH + Pago_digital + Am_lat + Credito_sqrt +
                                     Debito_sqrt + Densidad_Neg_sqrt + PIB_pc_PPA_sqrt + 
                                     ICH_sqrt, data = datos_adq_con_cont_enl) 

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MENL_aicc <- dredge(modelo_completo_enl, rank = "AICc", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MENL_aicc <- get.models(MENL_aicc, 5)[[1]] #Ranking=5
MENL_aicc

# Errores robustos
MENL_aicc_corr_hc <- coeftest(MENL_aicc, vcov. = vcovHC(MENL_aicc, "HC1", adjust = TRUE, sandwich = TRUE))

MENL_aicc_corr_hac <- coeftest(MENL_aicc, vcov. = vcovHAC(MENL_aicc, adjust = TRUE, sandwich = TRUE))

summary(MENL_aicc)
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MENL_aicc <- predict(MENL_aicc, newdata = data_CR_enl, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MENL_aicc`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación con efectos no lineales - Modelo AICc sin efectos fijos de tiempo")
pronostico <- c(pred_adq_MENL_aicc)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MENL_aicc, modelo, pred_adq_MENL_aicc, pronostico, MENL_aicc_corr_hc, MENL_aicc_corr_hac, var, modelo_completo_enl)
```

#### MENL_bic

##### Estimación

::: {style="text-align: justify"}
Se proporciona el código completo para su ejecución en R. No obstante, debido al tiempo excesivo que toma el procesamiento automático, se han incorporado manualmente las variables que optimizan el modelo. Estas variables fueron identificadas previamente mediante un proceso de selección realizado en Stata, con el fin de reducir significativamente el tiempo de ejecución en R.
:::

```{r, include=TRUE}
# Ajustar el modelo (automático R)
#modelo_completo_enl <- lm(formula_con_cont_enl, data = datos_adq_con_cont_enl)

# Ajustar el modelo (manual, selección previa de variables en Stata)
modelo_completo_enl <- lm(Adquirencia ~ Regulacion_intercambio + Cuenta + Credito +
                                     ICH + Pago_digital + Am_lat + Credito_sqrt + 
                                     Densidad_Neg_sqrt + ICH_sqrt, 
                                     data = datos_adq_con_cont_enl) 

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MENL_bic <- dredge(modelo_completo_enl, rank = "BIC", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MENL_bic <- get.models(MENL_bic, 3)[[1]] #Ranking=3
MENL_bic

# Errores robustos
MENL_bic_corr_hc <- coeftest(MENL_bic, vcov. = vcovHC(MENL_bic, "HC1", adjust = TRUE, sandwich = TRUE))

MENL_bic_corr_hac <- coeftest(MENL_bic, vcov. = vcovHAC(MENL_bic, adjust = TRUE, sandwich = TRUE))

summary(MENL_bic)
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MENL_bic <- predict(MENL_bic, newdata = data_CR_enl, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MENL_bic`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación con efectos no lineales - Modelo BIC sin efectos fijos de tiempo")
pronostico <- c(pred_adq_MENL_bic)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MENL_bic, modelo, pred_adq_MENL_bic, pronostico, MENL_bic_corr_hc, MENL_bic_corr_hac, var, modelo_completo_enl, formula_con_cont_enl)
```

#### MENL_aicc_et

##### Estimación

::: {style="text-align: justify"}
Se proporciona el código completo para su ejecución en R. No obstante, debido al tiempo excesivo que toma el procesamiento automático, se han incorporado manualmente las variables que optimizan el modelo. Estas variables fueron identificadas previamente mediante un proceso de selección realizado en Stata, con el fin de reducir significativamente el tiempo de ejecución en R.
:::

```{r, include= TRUE}
# Ajustar el modelo (automático R)
#modelo_completo_et_enl <- lm(formula_con_cont_t_enl, data = datos_adq_con_cont_enl)

# Ajustar el modelo (manual, selección previa de variables en Stata)
modelo_completo_et_enl <- lm(Adquirencia ~ Regulacion_intercambio + 
                                        Credito + ICH + Valor_prom + Am_lat + Credito_sqrt                                         + Densidad_Neg_sqrt + ICH_sqrt + Valor_prom_sqrt +                                         + Año_2014 + Año_2015 + Año_2016 + Año_2017 +  
                                        Año_2018 + Año_2019 + Año_2020 + Año_2021 + 
                                        Año_2022 + Año_2023 + Año_2024 + Año_2025, 
                                        data = datos_adq_con_cont_enl) 

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MENL_aicc_et <- dredge(modelo_completo_et_enl, rank = "AICc", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MENL_aicc_et <- get.models(MENL_aicc_et, 2)[[1]] #Ranking=2
MENL_aicc_et

# Errores robustos
MENL_aicc_et_corr_hc <- coeftest(MENL_aicc_et, vcov. = vcovHC(MENL_aicc_et, "HC1", adjust = TRUE, sandwich = TRUE))

MENL_aicc_et_corr_hac <- coeftest(MENL_aicc_et, vcov. = vcovHAC(MENL_aicc_et, adjust = TRUE, sandwich = TRUE))

summary(MENL_aicc_et)
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MENL_aicc_et <- predict(MENL_aicc_et, newdata = data_CR_enl, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MENL_aicc_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación con efectos no lineales - Modelo AICc con efectos fijos de tiempo")
pronostico <- c(pred_adq_MENL_aicc_et)
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MENL_aicc_et, modelo, pred_adq_MENL_aicc_et, pronostico, MENL_aicc_et_corr_hc, MENL_aicc_et_corr_hac, modelo_completo_et_enl)
```

#### MENL_bic_et

##### Estimación

::: {style="text-align: justify"}
Se proporciona el código completo para su ejecución en R. No obstante, debido al tiempo excesivo que toma el procesamiento automático, se han incorporado manualmente las variables que optimizan el modelo. Estas variables fueron identificadas previamente mediante un proceso de selección realizado en Stata, con el fin de reducir significativamente el tiempo de ejecución en R.
:::

```{r, include=FALSE}
# Ajustar el modelo (automático R)
#modelo_completo_et_enl <- lm(formula_con_cont_t_enl, data = datos_adq_con_cont_enl)

# Ajustar el modelo (manual, selección previa de variables en Stata)
modelo_completo_et_enl <- lm(Adquirencia ~ Regulacion_intercambio + Credito + ICH +
                                       Valor_prom + Am_lat + Credito_sqrt + 
                                       Densidad_Neg_sqrt + ICH_sqrt + Valor_prom_sqrt + 
                                       Año_2014 + Año_2015 + Año_2016 + Año_2017 +  
                                       Año_2018 + Año_2019 + Año_2020 + Año_2021 +  
                                       Año_2022 + Año_2023 + Año_2024 + Año_2025, 
                                       data = datos_adq_con_cont_enl) 

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MENL_bic_et <- dredge(modelo_completo_et_enl, rank = "BIC", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MENL_bic_et <- get.models(MENL_bic_et, 3)[[1]] #Ranking=3
MENL_bic_et

# Errores robustos
MENL_bic_et_corr_hc<-coeftest(MENL_bic_et, vcov. = vcovHC(MENL_bic_et, "HC1", adjust = TRUE, sandwich = TRUE))

MENL_bic_et_corr_hac<-coeftest(MENL_bic_et, vcov. = vcovHAC(MENL_bic_et, adjust = TRUE, sandwich = TRUE))
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MENL_bic_et <- predict(MENL_bic_et, newdata = data_CR_enl, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MENL_bic_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación con efectos no lineales - Modelo BIC con efectos fijos de tiempo")
pronostico <- c(pred_adq_MENL_bic_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MENL_bic_et, modelo, pred_adq_MENL_bic_et, pronostico, MENL_bic_et_corr_hc, MENL_bic_et_corr_hac, datos_adq_con_cont_enl, data_CR_enl, modelo_completo_et_enl, controles_numericos, controles_sqrt, formula_con_cont_t_enl)
```

### Método Bootstrap

::: {style="text-align: justify"}
Una manera alternativa de obtener los intervalos de confianza de las predicciones de las comisiones es mediante el método bootstrap. Este método de remuestreo permite obtener una estimación no paramétrica de la distribución empírica de las comisiones de intercambio y adquirencia. 

Para desarrollar este ejercicio, se realizan 1.000 remuestreos, y para cada muestra obtenida se sigue el mismo tipo de análisis presentado en la estimación principal. Esto da como resultado una comisión de adquirencia o intercambio promedio, según sea el caso, para cada iteración. 

A continuación, se presenta el código completo en R. No obstante, debido al tiempo estimado de ejecución (aproximadamente 16 horas), los chunks se encuentran deshabilitados.
:::

```{r, include= TRUE, eval=FALSE}
start_time <- Sys.time()
```

```{r, include= TRUE, eval=FALSE}
# Función segura para convertir factores con un solo nivel
convert_factors <- function(df) {
  df %>%
    mutate(across(where(is.factor), ~ if (nlevels(.) < 2) as.character(.) else .))
}

# Configurar plan de ejecución paralela (puedes ajustar el número de workers)


bootstrap_iter <- function(i) {
  cat("Iteración:", i, "\n")
  
  resultado <- list(
    Adquirencia_AICC = NA_real_,
    Adquirencia_BIC = NA_real_,
    Adquirencia_AICC_EF = NA_real_,
    Adquirencia_BIC_EF = NA_real_
  )
  
  datos <- readxl::read_excel("D:/Base_datos_adquirencia_intercambio.xlsx")
  
  datos <- datos %>%
    mutate(across(all_of(c("Regulacion_intercambio", "Regulacion_adquirencia",
                           "Europa", "Asia_Este", "Am_lat", "OCDE", "EEE",  
                           "Dato_mas_reciente",
                           paste0("Año_", 2010:2025))), as.factor))
  
  datos_int_con_cont <- datos %>%
    dplyr::select(-Intercambio) %>%
    na.omit()
  
  data_CR <- datos[which(datos$Codigo_Pais == "CRI"), which(names(datos) != "Intercambio")]
  
  controles <- c("Cajero_Autom", "Crimen", "Cuenta", "Credito", "Debito", "Densidad_Pob", "Densidad_Neg", "PIB_pc_PPA", "ICH", "Pago_digital", "Valor_prom", "Europa", "Asia_Este", "Am_lat")
  ef_tiempo <- paste0("Año_", 2013:2025)
  variables_fijas <- c("Regulacion_intercambio", paste0("Año_", 2010:2025))
  
  muestra <- tryCatch({
    temp <- datos_int_con_cont[sample(nrow(datos_int_con_cont), replace = TRUE), ]
    temp <- convert_factors(temp)
    temp <- droplevels(temp)
    
    # Detectar factores con un solo nivel
    factores_un_nivel <- names(Filter(function(x) is.factor(x) && nlevels(x) < 2, temp))
    
    # Ajustar variables de control y efectos fijos
    controles_validos <- setdiff(controles, factores_un_nivel)
    ef_tiempo_validos <- setdiff(ef_tiempo, factores_un_nivel)
    variables_fijas_validas <- setdiff(variables_fijas, factores_un_nivel)
    
    # Crear fórmulas dinámicas
    formula_con_cont <- as.formula(
      paste("Adquirencia ~ Regulacion_intercambio +", paste(controles_validos, collapse = " + "))
    )
    
    formula_con_cont_t <- as.formula(
      paste("Adquirencia ~ Regulacion_intercambio +", paste(c(controles_validos, ef_tiempo_validos), collapse = " + "))
    )
    
    list(
      muestra = temp,
      formula_con_cont = formula_con_cont,
      formula_con_cont_t = formula_con_cont_t,
      variables_fijas_validas = variables_fijas_validas
    )
    
  }, error = function(e) {
    cat("Error al generar muestra en iteración", i, ":", conditionMessage(e), "\n")
    return(NULL)
  })
  
  if (is.null(muestra)) return(resultado)
  
  temp <- muestra$muestra
  formula_con_cont <- muestra$formula_con_cont
  formula_con_cont_t <- muestra$formula_con_cont_t
  variables_fijas_validas <- muestra$variables_fijas_validas
  
  # AICC sin EF
  tryCatch({
    modelo <- lm(formula_con_cont, data = temp)
    options(na.action = "na.fail")
    modelos <- dredge(modelo, fixed = "Regulacion_intercambio", rank = "AICc")
    modelo <- get.models(modelos, 1)[[1]]
    datos_pred <- rbind(temp, data_CR)
    datos_pred$fitted <- predict(modelo, newdata = datos_pred)
    resultado$Adquirencia_AICC <- mean(datos_pred$fitted[datos_pred$Codigo_Pais == "CRI"], na.rm = TRUE)
  }, error = function(e) {
    cat("Error en AICC sin EF, iteración", i, ":", conditionMessage(e), "\n")
  })
  
  # BIC sin EF
  tryCatch({
    modelo <- lm(formula_con_cont, data = temp)
    options(na.action = "na.fail")
    modelos <- dredge(modelo, fixed = "Regulacion_intercambio", rank = "BIC")
    modelo <- get.models(modelos, 1)[[1]]
    datos_pred <- rbind(temp, data_CR)
    datos_pred$fitted <- predict(modelo, newdata = datos_pred)
    resultado$Adquirencia_BIC <- mean(datos_pred$fitted[datos_pred$Codigo_Pais == "CRI"], na.rm = TRUE)
  }, error = function(e) {
    cat("Error en BIC sin EF, iteración", i, ":", conditionMessage(e), "\n")
  })
  
  # AICC con EF
  tryCatch({
    modelo <- lm(formula_con_cont_t, data = temp)
    options(na.action = "na.fail")
    modelos <- dredge(modelo, fixed = variables_fijas_validas, rank = "AICc")
    modelo <- get.models(modelos, 1)[[1]]
    datos_pred <- rbind(temp, data_CR)
    datos_pred$fitted <- predict(modelo, newdata = datos_pred)
    resultado$Adquirencia_AICC_EF <- mean(datos_pred$fitted[datos_pred$Codigo_Pais == "CRI"], na.rm = TRUE)
  }, error = function(e) {
    cat("Error en AICC con EF, iteración", i, ":", conditionMessage(e), "\n")
  })
  
  # BIC con EF
  tryCatch({
    modelo <- lm(formula_con_cont_t, data = temp)
    options(na.action = "na.fail")
    modelos <- dredge(modelo, fixed = variables_fijas_validas, rank = "BIC")
    modelo <- get.models(modelos, 1)[[1]]
    datos_pred <- rbind(temp, data_CR)
    datos_pred$fitted <- predict(modelo, newdata = datos_pred)
    resultado$Adquirencia_BIC_EF <- mean(datos_pred$fitted[datos_pred$Codigo_Pais == "CRI"], na.rm = TRUE)
  }, error = function(e) {
    cat("Error en BIC con EF, iteración", i, ":", conditionMessage(e), "\n")
  })
  
  return(resultado)
}

plan(multisession, workers = parallel::detectCores() - 1)

# Número de iteraciones
n_iter <- 1000

# Ejecutar en paralelo con semilla reproducible
resultados_lista <- future_map(
  1:n_iter,
  bootstrap_iter,
  .progress = TRUE,
  .options = furrr_options(seed = 472195)
)

# Convertir a data.table
resultados <- rbindlist(resultados_lista, fill = TRUE)

```

```{r, include= TRUE, eval=FALSE}
## Post Bootstrap
# Calcular el promedio
resultados$Adquirencia <- rowMeans(resultados[, c("Adquirencia_AICC",
                                                  "Adquirencia_BIC",
                                                  "Adquirencia_AICC_EF",
                                                  "Adquirencia_BIC_EF")])


# Calcular estadísticas descriptivas
media <- mean(resultados$Adquirencia, na.rm = TRUE)
mediana <- median(resultados$Adquirencia, na.rm = TRUE)


# Formatear con dos decimales
media_fmt <- sprintf("%.2f", media)
mediana_fmt <- sprintf("%.2f", mediana)

# Mostrar resultados
cat("Media:", media_fmt, "\n")
cat("Mediana:", mediana_fmt, "\n")


#Obtener intervalo de confianza del 95% para Intercambio
c1 <- round(quantile(resultados$Adquirencia, 0.025), 2)
c2 <- round(quantile(resultados$Adquirencia, 0.975), 2)
ci <- quantile(resultados$Adquirencia, probs = c(0.025, 0.975))

# Mostrar los valores
cat("IC 95%:", c1, "-", c2, "\n")

#Generar gráfico para Intercambio

p <- ggplot(resultados, aes(x = Adquirencia)) +
  geom_histogram(aes(y = after_stat(density)), bins = 30, fill = "lightblue", color = "black") +
  theme_minimal(base_size = 14) +
  labs(
    x = "Comisión de adquirencia (%)",
    y = "Densidad",
    title = "Distribución de la Comisión de adquirencia"
  ) +
  annotate("text", x = media, y = 0.5, label = paste0("Media: ", round(media, 2), "%"), size = 5) +
  annotate("text", x = mediana, y = 0.4, label = paste0("Mediana: ", round(mediana, 2), "%"), size = 5) +
  annotate("text", x = mean(ci), y = 0.3,
           label = paste0("IC 95%: [", round(ci[1], 2), "%, ", round(ci[2], 2), "%]"),
           size = 5)
p

```

```{r, include= TRUE, eval=FALSE}
end_time <- Sys.time()
duration <- end_time - start_time
cat("Tiempo total de ejecución:", duration, "\n")
```

### Regresión Robusta

::: {style="text-align: justify"}
El método de regresión robusta permite analizar cambios en los resultados obtenidos producto de observaciones atípicas. Esto se consigue mediante la asignación de ponderadores a cada observación, que varían en sentido inverso a su nivel de influencia (medida por el tamaño de sus residuos), que cambian de manera iterativa hasta alcanzar un punto de convergencia.
:::

#### MRR_aicc **

##### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRR_aicc <- dredge(modelo_completo, rank = "AICc", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MRR_aicc <- get.models(MRR_aicc, 3)[[1]] #Ranking=3

# Extraer variables seleccionadas
vars_seleccionadas <- names(coef(MRR_aicc))[-1]  # sin intercepto

# Ajustar regresión robusta con las variables seleccionadas

# Paso 1: Ajuste robusto con Huber
formula_robusta <- as.formula(paste("Adquirencia ~", paste(vars_seleccionadas, collapse = " + ")))
modelo_huber <- rlm(formula_robusta, data = datos_adq_con_cont,  psi = psi.huber)


# Paso 2: Ajuste robusto con Bisquare usando coef iniciales de Huber
modelo_bisquare <- rlm(formula_robusta, data = datos_adq_con_cont, psi = psi.bisquare, init = coef(modelo_huber))

MRR_aicc <- modelo_bisquare

#Visualización
stargazer(MRR_aicc, type = "html", 
          title = "Regresión robusta - Modelo AICc sin efectos fijos de tiempo", 
          column.labels = c("MCO"), model.names = F, 
          dep.var.labels = c("Tasa adquirencia"), dep.var.caption = "Variable dependiente")

```

##### Predicción

```{r, include=TRUE}
pred_adq_MRR_aicc=predict(MRR_aicc, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MRR_aicc`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Regresión robusta - Modelo AICc sin efectos fijos de tiempo")
pronostico <- c(pred_adq_MRR_aicc)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRR_aicc, modelo, pred_adq_MRR_aicc, pronostico, modelo_bisquare, modelo_huber, vars_seleccionadas, formula_robusta)
```

#### MRR_bic **
##### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRR_bic <- dredge(modelo_completo, rank = "BIC", fixed = "Regulacion_intercambio") 

# Extraer el mejor modelo según AICc
MRR_bic <- get.models(MRR_bic, 20)[[1]] #Ranking=20

# Extraer variables seleccionadas
vars_seleccionadas <- names(coef(MRR_bic))[-1]  # sin intercepto

# Ajustar regresión robusta con las variables seleccionadas

# Paso 1: Ajuste robusto con Huber
formula_robusta <- as.formula(paste("Adquirencia ~", paste(vars_seleccionadas, collapse = " + ")))
modelo_huber <- rlm(formula_robusta, data = datos_adq_con_cont,  psi = psi.huber)


# Paso 2: Ajuste robusto con Bisquare usando coef iniciales de Huber
modelo_bisquare <- rlm(formula_robusta, data = datos_adq_con_cont, psi = psi.bisquare, init = coef(modelo_huber))

MRR_bic <- modelo_bisquare

#Visualización
stargazer(MRR_bic, type = "html", 
          title = "Regresión robusta - Modelo BIC sin efectos fijos de tiempo", 
          column.labels = c("MCO"), model.names = F, 
          dep.var.labels = c("Tasa adquirencia"), dep.var.caption = "Variable dependiente")

```

##### Predicción

```{r, include=TRUE}

pred_adq_MRR_bic=predict(MRR_bic, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MRR_bic`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Regresión robusta - Modelo BIC sin efectos fijos de tiempo")
pronostico <- c(pred_adq_MRR_bic)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRR_bic, modelo, pred_adq_MRR_bic, pronostico, modelo_bisquare, modelo_huber, vars_seleccionadas, formula_robusta)
```

#### MRR_aicc_et **

##### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRR_aicc_et <- dredge(modelo_completo_et, rank = "AICc", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MRR_aicc_et <- get.models(MRR_aicc_et, 11)[[1]] #Ranking=11

# Obtener coeficientes
coeficientes <- coef(MRR_aicc_et)

#Filtrar nombres de coeficientes que no sean NA
coeficientes_validos <- names(coeficientes[!is.na(coeficientes)])

# Extraer variables seleccionadas
vars_seleccionadas <- names(coef(MRR_aicc_et))[-1]  # sin intercepto
vars_seleccionadas <- intersect(coeficientes_validos, vars_seleccionadas)

# Ajustar regresión robusta con las variables seleccionadas
# Paso 1: Ajuste robusto con Huber
formula_robusta <- as.formula(paste("Adquirencia ~", paste(vars_seleccionadas, collapse = " + ")))

modelo_huber <- rlm(formula_robusta, data = datos_adq_con_cont,  psi = psi.huber)

# Paso 2: Ajuste robusto con Bisquare usando coef iniciales de Huber
modelo_bisquare <- rlm(formula_robusta, data = datos_adq_con_cont, psi = psi.bisquare, init = coef(modelo_huber))

MRR_aicc_et <- modelo_bisquare

#Visualización
stargazer(MRR_aicc_et, type = "html", 
          title = "Regresión robusta - Modelo AICc con efectos fijos de tiempo", 
          column.labels = c("MCO"), model.names = F, 
          dep.var.labels = c("Tasa adquirencia"), dep.var.caption = "Variable dependiente")

```

##### Predicción

```{r, include=TRUE}
pred_adq_MRR_aicc_et=predict(MRR_aicc_et, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MRR_aicc_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Regresión robusta - Modelo AICc con efectos fijos de tiempo")
pronostico <- c(pred_adq_MRR_aicc_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRR_aicc_et, modelo, pred_adq_MRR_aicc_et, pronostico, modelo_bisquare, modelo_huber, vars_seleccionadas, formula_robusta, coeficientes, coeficientes_validos)
```

#### MRR_bic_et **

##### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MRR_bic_et <- dredge(modelo_completo_et, rank = "BIC", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MRR_bic_et <- get.models(MRR_bic_et, 10)[[1]] #Ranking=10

# Obtener coeficientes
coeficientes <- coef(MRR_bic_et)

#Filtrar nombres de coeficientes que no sean NA
coeficientes_validos <- names(coeficientes[!is.na(coeficientes)])

# Extraer variables seleccionadas
vars_seleccionadas <- names(coef(MRR_bic_et))[-1]  # sin intercepto
vars_seleccionadas <- intersect(coeficientes_validos, vars_seleccionadas)

# Ajustar regresión robusta con las variables seleccionadas
# Paso 1: Ajuste robusto con Huber
formula_robusta <- as.formula(paste("Adquirencia ~", paste(vars_seleccionadas, collapse = " + ")))

modelo_huber <- rlm(formula_robusta, data = datos_adq_con_cont,  psi = psi.huber)

# Paso 2: Ajuste robusto con Bisquare usando coef iniciales de Huber
modelo_bisquare <- rlm(formula_robusta, data = datos_adq_con_cont, psi = psi.bisquare, init = coef(modelo_huber))

MRR_bic_et <- modelo_bisquare

#Visualización
stargazer(MRR_bic_et, type = "html", 
          title = "Regresión robusta - Modelo BIC con efectos fijos de tiempo", 
          column.labels = c("MCO"), model.names = F, 
          dep.var.labels = c("Tasa adquirencia"), dep.var.caption = "Variable dependiente")

```

##### Predicción

```{r, include=TRUE}
pred_adq_MRR_bic_et <- predict(MRR_bic_et, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MRR_bic_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Regresión robusta - Modelo BIC con efectos fijos de tiempo")
pronostico <- c(pred_adq_MRR_bic_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MRR_bic_et, modelo, pred_adq_MRR_bic_et, pronostico, modelo_bisquare, modelo_huber, vars_seleccionadas, formula_robusta, coeficientes, coeficientes_validos)
```

### Restringiendo la muestra a solo países regulados

::: {style="text-align: justify"}
Como una prueba adicional, se estima el modelo principal, pero luego de restringir la muestra a los países que al igual que Costa Rica, son regulados.
:::

```{r, include=FALSE}
#Datos

#Se filtran los datos para excluir países no regulados
datos_reg <- datos[datos$Regulacion_intercambio == 1,]
datos_reg <- datos_reg %>%
  select(all_of(variables_interes)) %>%
  na.omit()
```

#### MPR_aicc

##### Estimación

```{r, include=TRUE, results='asis'}
# Ajustar el modelo
modelo_completo_reg <- lm(formula_con_cont, data = datos_reg)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MPR_aicc <- dredge(modelo_completo_reg, rank = "AICc", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MPR_aicc <- get.models(MPR_aicc, 3)[[1]] #Ranking=3
MPR_aicc

# Errores robustos
MPR_aicc_corr_hc <- coeftest(MPR_aicc, vcov. = vcovHC(MPR_aicc, "HC1", adjust = TRUE, sandwich = TRUE))

MPR_aicc_corr_hac <- coeftest(MPR_aicc, vcov. = vcovHAC(MPR_aicc, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MPR_aicc, MPR_aicc_corr_hc,
          type = "html", 
          title = "Estimación restringiendo muestra a países regulados - Modelo AICc sin efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa adquirencia", "Tasa adquirencia"), 
          dep.var.caption = "Variable dependiente")
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MPR_aicc <- predict(MPR_aicc, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MPR_aicc`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación restringiendo muestra a países regulados - Modelo AICc sin efectos fijos de tiempo")
pronostico <- c(pred_adq_MPR_aicc)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MPR_aicc, modelo, pred_adq_MPR_aicc, pronostico, MPR_aicc_corr_hc, MPR_aicc_corr_hac)
```

#### MPR_bic

##### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MPR_bic <- dredge(modelo_completo_reg, rank = "BIC", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MPR_bic <- get.models(MPR_bic, 2)[[1]] #Ranking=2
MPR_bic

# Errores robustos
MPR_bic_corr_hc <- coeftest(MPR_bic, vcov. = vcovHC(MPR_bic, "HC1", adjust = TRUE, sandwich = TRUE)) 

MPR_bic_corr_hac <- coeftest(MPR_bic, vcov. = vcovHAC(MPR_bic, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MPR_bic, MPR_bic_corr_hc, type = "html",
          title = "Estimación restringiendo muestra a países regulados - Modelo BIC sin efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F, 
          dep.var.labels = c("Tasa adquirencia", "Tasa adquirencia"), 
          dep.var.caption = "Variable dependiente")

```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MPR_bic <- predict(MPR_bic, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MPR_bic`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación restringiendo muestra a países regulados - Modelo BIC sin efectos fijos de tiempo")
pronostico <- c(pred_adq_MPR_bic)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MPR_bic, modelo, pred_adq_MPR_bic, pronostico, MPR_bic_corr_hc, MPR_bic_corr_hac, modelo_completo_reg)
```

#### MPR_aicc_et

##### Estimación

```{r, include=TRUE, results='asis'}
# Ajustar el modelo
#modelo_completo_et_reg <- lm(formula_con_cont_t, data = datos_reg)
modelo_completo_et_reg <- lm(Adquirencia ~ Regulacion_intercambio + Cajero_Autom + Cuenta 
                             + Debito + Densidad_Pob + Asia_Este + Europa + Año_2015 + 
                             Año_2016 + Año_2017 + Año_2018 + Año_2019 + Año_2020 + 
                             Año_2021 + Año_2022 + Año_2023 + Año_2024, data = datos_reg)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MPR_aicc_et <- dredge(modelo_completo_et_reg, rank = "AICc", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MPR_aicc_et <- get.models(MPR_aicc_et, 1)[[1]] 
MPR_aicc_et

# Errores robustos
MPR_aicc_et_corr_hc <- coeftest(MPR_aicc_et, vcov. = vcovHC(MPR_aicc_et, "HC1", adjust = TRUE, sandwich = TRUE))

MPR_aicc_et_corr_hac <- coeftest(MPR_aicc_et, vcov. = vcovHAC(MPR_aicc_et, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MPR_aicc_et, MPR_aicc_et_corr_hc,
          type = "html", 
          title = "Estimación restringiendo muestra a países regulados - Modelo AICc con efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa adquirencia", "Tasa adquirencia"), 
          dep.var.caption = "Variable dependiente")
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MPR_aicc_et <- predict(MPR_aicc_et, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MPR_aicc_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación restringiendo muestra a países regulados - Modelo AICc con efectos fijos de tiempo")
pronostico <- c(pred_adq_MPR_aicc_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MPR_aicc_et, modelo, pred_adq_MPR_aicc_et, pronostico, MPR_aicc_et_corr_hc, MPR_aicc_et_corr_hac)
```

#### MPR_bic_et

##### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MPR_bic_et <- dredge(modelo_completo_et_reg, rank = "BIC", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MPR_bic_et <- get.models(MPR_bic_et, 1)[[1]] 
MPR_bic_et

# Errores robustos
MPR_bic_et_corr_hc <- coeftest(MPR_bic_et, vcov. = vcovHC(MPR_bic_et, "HC1", adjust = TRUE, sandwich = TRUE))

MPR_bic_et_corr_hac <- coeftest(MPR_bic_et, vcov. = vcovHAC(MPR_bic_et, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MPR_bic_et, MPR_bic_et_corr_hc,
          type = "html", 
          title = "Estimación restringiendo muestra a países regulados - Modelo BIC con efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa adquirencia", "Tasa adquirencia"), 
          dep.var.caption = "Variable dependiente")
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MPR_bic_et <- predict(MPR_bic_et, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MPR_bic_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación restringiendo muestra a países regulados - Modelo BIC sin efectos fijos de tiempo")
pronostico <- c(pred_adq_MPR_bic_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MPR_bic_et, modelo, pred_adq_MPR_bic_et, pronostico, MPR_bic_et_corr_hc, MPR_bic_et_corr_hac, datos_reg, modelo_completo_et_reg)
```

### Ampliando la definición de regulación de intercambio o adquirencia

::: {style="text-align: justify"}
En el análisis principal, la variable regulación hace referencia a regulación en las comisiones de intercambio. Esto se debe a que son relativamente pocos los países que regulan adquirencia. Sin embargo, es posible generar una variable dicotómica que haga referencia a que el país en cuestión tiene regulación sobre intercambio o adquirencia. 
:::

```{r, include=TRUE}
# Arreglo de datos
datos_reg_amp <- datos

# Crear una nueva variable lógica que sea TRUE si se cumple alguna de las dos condiciones

datos_reg_amp$Regulacion_inter_adq <- ifelse(datos_reg_amp$Regulacion_intercambio == 1 | datos_reg_amp$Regulacion_adquirencia == 1, 1, 0)

#Se filtran los datos
datos_reg_amp <- datos_reg_amp %>%
  select(all_of(variables_interes), Regulacion_inter_adq) %>% 
  select(-Regulacion_intercambio) %>%
  na.omit()

# Crear una nueva variable lógica que sea TRUE si se cumple alguna de las dos condiciones
data_CR_reg_amp <- data_CR
data_CR_reg_amp$Regulacion_inter_adq <- ifelse(data_CR_reg_amp$Regulacion_intercambio == 1 | data_CR_reg_amp$Regulacion_adquirencia == 1, 1, 0)

# Incluir las variables de control sin controles de tiempo
formula_con_cont_amp <- as.formula(
  paste("Adquirencia ~ Regulacion_inter_adq +", paste(controles, collapse = " + "))
)

# Incluir las variables de control con controles de tiempo
formula_con_cont_t_amp <- as.formula(
  paste("Adquirencia ~ Regulacion_inter_adq +", paste(controles, collapse = " + ")," + ", paste(ef_tiempo, collapse = " + "))
)

# Variables fijas
variables_fijas_amp <- c("Regulacion_inter_adq", paste0("Año_", 2010:2025))
```

#### MADR_aicc

##### Estimación

```{r, include=TRUE, results='asis'}
# Ajustar el modelo
modelo_completo_reg_amp <- lm(formula_con_cont_amp, data = datos_reg_amp)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MADR_aicc <- dredge(modelo_completo_reg_amp, rank = "AICc", fixed = "Regulacion_inter_adq")

# Extraer el mejor modelo según AICc
MADR_aicc <- get.models(MADR_aicc, 2)[[1]] #Ranking=2
MADR_aicc

# Errores robustos
MADR_aicc_corr_hc <- coeftest(MADR_aicc, vcov. = vcovHC(MADR_aicc, "HC1", adjust = TRUE, sandwich = TRUE))

MADR_aicc_corr_hac <- coeftest(MADR_aicc, vcov. = vcovHAC(MADR_aicc, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MADR_aicc, MADR_aicc_corr_hc,
          type = "html", 
          title = "Estimación ampliando definición de regulación - Modelo AICc sin efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa adquirencia", "Tasa adquirencia"), 
          dep.var.caption = "Variable dependiente")
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MADR_aicc <- predict(MADR_aicc, newdata = data_CR_reg_amp, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MADR_aicc`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación ampliando definición de regulación - Modelo AICc sin efectos fijos de tiempo")
pronostico <- c(pred_adq_MADR_aicc)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MADR_aicc, modelo, pred_adq_MADR_aicc, pronostico, MADR_aicc_corr_hc, MADR_aicc_corr_hac)
```

#### MADR_bic

##### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MADR_bic <- dredge(modelo_completo_reg_amp, rank = "BIC", fixed = "Regulacion_inter_adq")

# Extraer el mejor modelo según AICc
MADR_bic <- get.models(MADR_bic, 15)[[1]] #Ranking=15
MADR_bic

# Errores robustos
MADR_bic_corr_hc <- coeftest(MADR_bic, vcov. = vcovHC(MADR_bic, "HC1", adjust = TRUE, sandwich = TRUE))

MADR_bic_corr_hac <- coeftest(MADR_bic, vcov. = vcovHAC(MADR_bic, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MADR_bic, MADR_bic_corr_hc, type = "html",
          title = "Estimación ampliando definición de regulación - Modelo BIC sin efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F, 
          dep.var.labels = c("Tasa adquirencia", "Tasa adquirencia"), 
          dep.var.caption = "Variable dependiente")

```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MADR_bic <- predict(MADR_bic, newdata = data_CR_reg_amp, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MADR_bic`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación ampliando definición de regulación - Modelo BIC sin efectos fijos de tiempo")
pronostico <- c(pred_adq_MADR_bic)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MADR_bic, modelo, pred_adq_MADR_bic, pronostico, MADR_bic_corr_hc, MADR_bic_corr_hac, modelo_completo_reg_amp, formula_con_cont_amp)
```

#### MADR_aicc_et **

##### Estimación

```{r, include=TRUE, results='asis'}
# Ajustar el modelo
#modelo_completo_t_reg_amp <- lm(formula_con_cont_t_amp, data = datos_reg_amp)
modelo_completo_t_reg_amp <- lm(Adquirencia ~ Regulacion_inter_adq + Debito + ICH + 
                                  Valor_prom + Asia_Este + Am_lat + Año_2014 + Año_2015 + 
                                  Año_2016 + Año_2017 + Año_2018 + Año_2019 + Año_2020 + 
                                  Año_2021 + Año_2022 + Año_2023 + Año_2024 + Año_2025, 
                                  data = datos_reg_amp)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MADR_aicc_et <- dredge(modelo_completo_t_reg_amp, rank = "AICc", fixed = variables_fijas_amp)

# Extraer el mejor modelo según AICc
MADR_aicc_et <- get.models(MADR_aicc_et, 1)[[1]]
MADR_aicc_et

# Errores robustos
MADR_aicc_et_corr_hc <- coeftest(MADR_aicc_et, vcov. = vcovHC(MADR_aicc_et, "HC1", adjust = TRUE, sandwich = TRUE))

MADR_aicc_et_corr_hac <- coeftest(MADR_aicc_et, vcov. = vcovHAC(MADR_aicc_et, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MADR_aicc_et, MADR_aicc_et_corr_hc,
          type = "html", 
          title = "Estimación ampliando definición de regulación - Modelo AICc con efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa adquirencia", "Tasa adquirencia"), 
          dep.var.caption = "Variable dependiente")
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MADR_aicc_et <- predict(MADR_aicc_et, newdata = data_CR_reg_amp, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MADR_aicc_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación ampliando definición de regulación - Modelo AICc con efectos fijos de tiempo")
pronostico <- c(pred_adq_MADR_aicc_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MADR_aicc_et, modelo, pred_adq_MADR_aicc_et, pronostico, MADR_aicc_et_corr_hc, MADR_aicc_et_corr_hac)
```

#### MADR_bic_et

##### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MADR_bic_et <- dredge(modelo_completo_t_reg_amp, rank = "BIC", fixed = variables_fijas_amp)

# Extraer el mejor modelo según AICc
MADR_bic_et <- get.models(MADR_bic_et, 1)[[1]]
MADR_bic_et

# Errores robustos
MADR_bic_et_corr_hc <- coeftest(MADR_bic_et, vcov. = vcovHC(MADR_bic_et, "HC1", adjust = TRUE, sandwich = TRUE))

MADR_bic_et_corr_hac <- coeftest(MADR_bic_et, vcov. = vcovHAC(MADR_bic_et, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MADR_bic_et, MADR_bic_et_corr_hc,
          type = "html", 
          title = "Estimación ampliando definición de regulación - Modelo BIC con efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa adquirencia", "Tasa adquirencia"), 
          dep.var.caption = "Variable dependiente")
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MADR_bic_et <- predict(MADR_bic_et, newdata = data_CR_reg_amp, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MADR_bic_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación ampliando definición de regulación - Modelo BIC con efectos fijos de tiempo")
pronostico <- c(pred_adq_MADR_bic_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MADR_bic_et, modelo, pred_adq_MADR_bic_et, pronostico, MADR_bic_et_corr_hc, MADR_bic_et_corr_hac, data_CR_reg_amp, datos_reg_amp, modelo_completo_t_reg_amp, formula_con_cont_t_amp, variables_fijas_amp)
```

### Regresión que utilice las características observadas más recientes para cada país

::: {style="text-align: justify"}
La estimación principal utiliza las características observadas para los países al momento de la fijación de determinado porcentaje de comisión. Una forma alternativa de estimar el modelo es considerar para cada país las comisiones de intercambio o de adquirencia vigentes y las características observables más recientes, que en la mayoría de los países corresponde a información para el 2024.
:::

```{r, include=TRUE}
#Datos
datos_2024 <- read_excel("D:/Base_datos_adquirencia_intercambio_2024.xlsx")

datos_2024 <- datos_2024 %>% 
  dplyr :: select(-Intercambio) %>% 
  na.omit
```

#### MCR_aicc

##### Estimación

```{r, include=TRUE, results='asis'}
# Ajustar el modelo
modelo_completo_2024 <- lm(formula_con_cont, data = datos_2024)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MCR_aicc <- dredge(modelo_completo_2024, rank = "AICc", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MCR_aicc <- get.models(MCR_aicc, 3)[[1]] #Ranking=3
MCR_aicc

# Errores robustos
MCR_aicc_corr_hc <- coeftest(MCR_aicc, vcov. = vcovHC(MCR_aicc, "HC1", adjust = TRUE, sandwich = TRUE))

MCR_aicc_corr_hac <- coeftest(MCR_aicc, vcov. = vcovHAC(MCR_aicc, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MCR_aicc, MCR_aicc_corr_hc,
          type = "html", 
          title = "Estimación usando características más recientes - Modelo AICc sin efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa adquirencia", "Tasa adquirencia"), 
          dep.var.caption = "Variable dependiente")
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MCR_aicc <- predict(MCR_aicc, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MCR_aicc`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación usando características más recientes - Modelo AICc sin efectos fijos de tiempo")
pronostico <- c(pred_adq_MCR_aicc)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MCR_aicc, modelo, pred_adq_MCR_aicc, pronostico, MCR_aicc_corr_hc, MCR_aicc_corr_hac)
```

#### MCR_bic

##### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MCR_bic <- dredge(modelo_completo_2024, rank = "BIC", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MCR_bic <- get.models(MCR_bic, 4)[[1]] #Ranking=4
MCR_bic

# Errores robustos
MCR_bic_corr_hc <- coeftest(MCR_bic, vcov. = vcovHC(MCR_bic, "HC1", adjust = TRUE, sandwich = TRUE))

MCR_bic_corr_hac <- coeftest(MCR_bic, vcov. = vcovHAC(MCR_bic, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MCR_bic, MCR_bic_corr_hc,
          type = "html", 
          title = "Estimación usando características más recientes - Modelo BIC sin efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa adquirencia", "Tasa adquirencia"), 
          dep.var.caption = "Variable dependiente")
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MCR_bic <- predict(MCR_bic, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MCR_bic`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación usando características más recientes - Modelo BIC sin efectos fijos de tiempo")
pronostico <- c(pred_adq_MCR_bic)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MCR_bic, modelo, pred_adq_MCR_bic, pronostico, MCR_bic_corr_hc, MCR_bic_corr_hac, datos_2024, modelo_completo_2024)
```

### Regresión con variables propuestas por el Fondo Monetario Internacional (FMI)

::: {style="text-align: justify"}
El Fondo Monetario Internacional (FMI) publicó en diciembre de 2023 el estudio “Determinants of Card Payments Fees in A Cross-Country Analysis”, donde se replica y amplía el análisis de comparación internacional que realiza el BCCR para establecer los topes máximos de comisiones. La forma de ampliar el estudio fue incorporar las variables margen de interés neto y el grado de concentración de activos bancarios.
:::

```{r, include=TRUE}
# Datos
datos_fmi <- datos %>%
  select(all_of(variables_interes), c("Margen_neto", "Concentracion_bancaria")) %>%
  na.omit()

#Definir el conjunto de posibles variables de control

controles_fmi <- c("Cajero_Autom", "Crimen", "Cuenta", "Credito", "Debito", "Densidad_Pob", "Densidad_Neg", "PIB_pc_PPA", "ICH", "Pago_digital", "Valor_prom", "Europa", "Asia_Este", "Am_lat", "Margen_neto", "Concentracion_bancaria")

# Incluir las variables de control
formula_con_cont_fmi <- as.formula(
  paste("Adquirencia ~ Regulacion_intercambio +", paste(controles_fmi, collapse = " + "))
)

# Incluir las variables de control
formula_con_cont_t_fmi <- as.formula(
  paste("Adquirencia ~ Regulacion_intercambio +", paste(controles_fmi, collapse = " + ")," + ", paste(ef_tiempo, collapse = " + "))
)

```

#### MFMI_aicc

##### Estimación

```{r, include=TRUE, results='asis'}
# Ajustar el modelo
modelo_completo_fmi <- lm(formula_con_cont_fmi, data = datos_fmi)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MFMI_aicc <- dredge(modelo_completo_fmi, rank = "AICc", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MFMI_aicc <- get.models(MFMI_aicc, 1)[[1]]
MFMI_aicc

# Errores robustos
MFMI_aicc_corr_hc <- coeftest(MFMI_aicc, vcov. = vcovHC(MFMI_aicc, "HC1", adjust = TRUE, sandwich = TRUE))

MFMI_aicc_corr_hac <- coeftest(MFMI_aicc, vcov. = vcovHAC(MFMI_aicc, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MFMI_aicc, MFMI_aicc_corr_hc,
          type = "html", 
          title = "Réplica análisis FMI - Modelo AICc sin efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa adquirencia", "Tasa adquirencia"), 
          dep.var.caption = "Variable dependiente")
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MFMI_aicc <- predict(MFMI_aicc, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MFMI_aicc`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Réplica análisis FMI - Modelo AICc sin efectos fijos de tiempo")
pronostico <- c(pred_adq_MFMI_aicc)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MFMI_aicc, modelo, pred_adq_MFMI_aicc, pronostico, MFMI_aicc_corr_hac, MFMI_aicc_corr_hc)
```

#### MFMI_bic

##### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MFMI_bic <- dredge(modelo_completo_fmi, rank = "BIC", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MFMI_bic <- get.models(MFMI_bic, 6)[[1]] #Ranking=6
MFMI_bic

# Errores robustos
MFMI_bic_corr_hc <- coeftest(MFMI_bic, vcov. = vcovHC(MFMI_bic, "HC1", adjust = TRUE, sandwich = TRUE))

MFMI_bic_corr_hac <- coeftest(MFMI_bic, vcov. = vcovHAC(MFMI_bic, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MFMI_bic, MFMI_bic_corr_hc,
          type = "html", 
          title = "Réplica análisis FMI - Modelo BIC sin efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa adquirencia", "Tasa adquirencia"), 
          dep.var.caption = "Variable dependiente")
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MFMI_bic <- predict(MFMI_bic, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MFMI_bic`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Réplica análisis FMI - Modelo BIC sin efectos fijos de tiempo")
pronostico <- c(pred_adq_MFMI_bic)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MFMI_bic, modelo, pred_adq_MFMI_bic, pronostico, MFMI_bic_corr_hac, MFMI_bic_corr_hc, modelo_completo_fmi, formula_con_cont_fmi)
```

#### MFMI_aicc_et

##### Estimación

```{r, include=TRUE, results='asis'}
# Ajustar el modelo
#modelo_completo_fmi_et <- lm(formula_con_cont_t_fmi, data = datos_fmi)
modelo_completo_fmi_et <- lm(Adquirencia ~ Regulacion_intercambio + Cajero_Autom + Cuenta 
                             + Valor_prom +Asia_Este + Margen_neto + 
                             Concentracion_bancaria + Año_2014 + Año_2015 + Año_2016 + 
                             Año_2017 + Año_2018 + Año_2019 + Año_2020 + Año_2021 + 
                             Año_2022 + Año_2023 + Año_2024 + Año_2025, data = datos_fmi)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MFMI_aicc_et <- dredge(modelo_completo_fmi_et, rank = "AICc", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MFMI_aicc_et <- get.models(MFMI_aicc_et, 1)[[1]]
MFMI_aicc_et

# Errores robustos
MFMI_aicc_et_corr_hc <- coeftest(MFMI_aicc_et, vcov. = vcovHC(MFMI_aicc_et, "HC1", adjust = TRUE, sandwich = TRUE))

MFMI_aicc_et_corr_hac <- coeftest(MFMI_aicc_et, vcov. = vcovHAC(MFMI_aicc_et, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MFMI_aicc_et, MFMI_aicc_et_corr_hc,
          type = "html", 
          title = "Réplica análisis FMI - Modelo AICc con efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa adquirencia", "Tasa adquirencia"), 
          dep.var.caption = "Variable dependiente")
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MFMI_aicc_et <- predict(MFMI_aicc_et, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MFMI_aicc_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Réplica análisis FMI - Modelo AICc con efectos fijos de tiempo")
pronostico <- c(pred_adq_MFMI_aicc_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MFMI_aicc_et, modelo, pred_adq_MFMI_aicc_et, pronostico, MFMI_aicc_et_corr_hac, MFMI_aicc_et_corr_hc)
```

#### MFMI_bic_et 

##### Estimación

```{r, include=TRUE, results='asis'}
# Ajustar el modelo
#modelo_completo_fmi_et <- lm(formula_con_cont_t_fmi, data = datos_fmi)
modelo_completo_fmi_et <- lm(Adquirencia ~ Regulacion_intercambio + Cajero_Autom + Cuenta 
                             + Valor_prom +Asia_Este + Margen_neto + 
                             Año_2014 + Año_2015 + Año_2016 + 
                             Año_2017 + Año_2018 + Año_2019 + Año_2020 + Año_2021 + 
                             Año_2022 + Año_2023 + Año_2024 + Año_2025, data = datos_fmi)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MFMI_bic_et <- dredge(modelo_completo_fmi_et, rank = "BIC", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MFMI_bic_et <- get.models(MFMI_bic_et, 1)[[1]]
MFMI_bic_et

# Errores robustos
MFMI_bic_et_corr_hc <- coeftest(MFMI_bic_et, vcov. = vcovHC(MFMI_bic_et, "HC1", adjust = TRUE, sandwich = TRUE))

MFMI_bic_et_corr_hac <- coeftest(MFMI_bic_et, vcov. = vcovHAC(MFMI_bic_et, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MFMI_bic_et, MFMI_bic_et_corr_hc,
          type = "html", 
          title = "Réplica análisis FMI - Modelo BIC con efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa adquirencia", "Tasa adquirencia"), 
          dep.var.caption = "Variable dependiente")
```

##### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MFMI_bic_et <- predict(MFMI_bic_et, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MFMI_bic_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Réplica análisis FMI - Modelo BIC con efectos fijos de tiempo")
pronostico <- c(pred_adq_MFMI_bic_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MFMI_bic_et, modelo, pred_adq_MFMI_bic_et, pronostico, MFMI_bic_et_corr_hac, MFMI_bic_et_corr_hc, datos_fmi, modelo_completo_fmi_et, controles_fmi, formula_con_cont_t_fmi)
```

### Regresión con variables instrumentales

::: {style="text-align: justify"}
Una posible observación es que las comisiones de intercambio y adquirencia afectan la adopción de tarjetas de pago, y estas a su vez afectan las comisiones, lo que introduciría un sesgo de endogeneidad en las estimaciones. Para atender esta observación se recurre a dos enfoques.

El primer enfoque consiste en estimar la especificación principal, pero sin considerar como variables explicativas la información relativa directamente al uso de tarjetas.

Para el segundo enfoque, se plantea un modelo de variables instrumentales. En este caso, se instrumenta el valor promedio de la transacción realizada con tarjetas de pago con el Índice de Capital Humano.
:::

#### Enfoque1: regresión sin variables relacionadas al uso de tarjetas

::: {style="text-align: justify"}
En esta regresión se estima la especificación principal, pero sin considerar como variables explicativas la información relativa directamente al uso de tarjetas. De esta manera, se excluye el porcentaje de adultos con tarjeta de crédito o débito, el porcentaje de adultos con cuenta en una institución financiera, el porcentaje de adultos que realizaron o recibieron un pago digital, el número de cajeros automáticos, y el valor de la transacción promedio realizada con tarjetas de pago.  
:::

```{r, include=FALSE}
#Definir el conjunto de posibles variables de control

controles_red <- c("Crimen", "Densidad_Pob", "Densidad_Neg", "PIB_pc_PPA", "ICH", "Europa", "Asia_Este", "Am_lat")

# Incluir las variables de control
formula_con_cont_red <- as.formula(
  paste("Adquirencia ~ Regulacion_intercambio +", paste(controles_red, collapse = " + "))
)

# Incluir las variables de control
formula_con_cont_t_red <- as.formula(
  paste("Adquirencia ~ Regulacion_intercambio +", paste(controles_red, collapse = " + ")," + ", paste(ef_tiempo, collapse = " + "))
)

#Datos sin variables relacionadas al uso de tarjetas
variables_interes_red <- c(variables_constantes, controles_red, ef_tiempo)

datos_red <- datos %>%
  select(all_of(variables_interes_red)) %>%
  na.omit()
```

##### MSVR_aicc

###### Estimación

```{r, include=TRUE, results='asis'}
# Ajustar el modelo
modelo_completo_red <- lm(formula_con_cont_red, data = datos_red)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MSVR_aicc <- dredge(modelo_completo_red, rank = "AICc", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MSVR_aicc <- get.models(MSVR_aicc, 1)[[1]]
MSVR_aicc

# Errores robustos
MSVR_aicc_corr_hc <- coeftest(MSVR_aicc, vcov. = vcovHC(MSVR_aicc, "HC1", adjust = TRUE, sandwich = TRUE))

MSVR_aicc_corr_hac <- coeftest(MSVR_aicc, vcov. = vcovHAC(MSVR_aicc, adjust = TRUE, sandwich = TRUE)) 

# Visualización
stargazer(MSVR_aicc, MSVR_aicc_corr_hc,
          type = "html", 
          title = "Regresión sin variables relacionadas al uso de tarjetas - Modelo AICc sin efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa adquirencia", "Tasa adquirencia"), 
          dep.var.caption = "Variable dependiente")
```

###### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MSVR_aicc <- predict(MSVR_aicc, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MSVR_aicc`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Regresión sin variables relacionadas al uso de tarjetas - Modelo AICc sin efectos fijos de tiempo")
pronostico <- c(pred_adq_MSVR_aicc)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MSVR_aicc, modelo, pred_adq_MSVR_aicc, pronostico, MSVR_aicc_corr_hc, MSVR_aicc_corr_hac)
```

##### MSVR_bic

###### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MSVR_bic <- dredge(modelo_completo_red, rank = "BIC", fixed = "Regulacion_intercambio")

# Extraer el mejor modelo según AICc
MSVR_bic <- get.models(MSVR_bic, 1)[[1]]
MSVR_bic

# Errores robustos
MSVR_bic_corr_hc <- coeftest(MSVR_bic, vcov. = vcovHC(MSVR_bic, "HC1", adjust = TRUE, sandwich = TRUE))

MSVR_bic_corr_hac <- coeftest(MSVR_bic, vcov. = vcovHAC(MSVR_bic, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MSVR_bic, MSVR_bic_corr_hc,
          type = "html", 
          title = "Regresión sin variables relacionadas al uso de tarjetas - Modelo BIC  sin efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa adquirencia", "Tasa adquirencia"), 
          dep.var.caption = "Variable dependiente")
```

###### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MSVR_bic <- predict(MSVR_bic, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MSVR_bic`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Regresión sin variables relacionadas al uso de tarjetas - Modelo BIC sin efectos fijos de tiempo")
pronostico <- c(pred_adq_MSVR_bic)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MSVR_bic, modelo, pred_adq_MSVR_bic, pronostico, MSVR_bic_corr_hc, MSVR_bic_corr_hac, modelo_completo_red, formula_con_cont_red)
```

##### MSVR_aicc_et **

###### Estimación

```{r, include=TRUE, results='asis'}
# Ajustar el modelo
modelo_completo_et_red <- lm(formula_con_cont_t_red, data = datos_red)

# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MSVR_aicc_et <- dredge(modelo_completo_et_red, rank = "AICc", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MSVR_aicc_et <- get.models(MSVR_aicc_et, 1)[[1]]
MSVR_aicc_et

# Errores robustos
MSVR_aicc_et_corr_hc<-coeftest(MSVR_aicc_et, vcov. = vcovHC(MSVR_aicc_et, "HC1", adjust = TRUE, sandwich = TRUE))

MSVR_aicc_et_corr_hac<-coeftest(MSVR_aicc_et, vcov. = vcovHAC(MSVR_aicc_et, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MSVR_aicc_et, MSVR_aicc_et_corr_hc,
          type = "html", 
          title = "Regresión sin variables relacionadas al uso de tarjetas - Modelo AICc con efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa adquirencia", "Tasa adquirencia"), 
          dep.var.caption = "Variable dependiente")
```

###### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MSVR_aicc_et <- predict(MSVR_aicc_et, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MSVR_aicc_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Regresión sin variables relacionadas al uso de tarjetas - Modelo AICc con efectos fijos de tiempo")
pronostico <- c(pred_adq_MSVR_aicc_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MSVR_aicc_et, modelo, pred_adq_MSVR_aicc_et, pronostico, MSVR_aicc_et_corr_hc, MSVR_aicc_et_corr_hac)
```

##### MSVR_bic_et **

###### Estimación

```{r, include=TRUE, results='asis'}
# Esto es necesario para que dredge funcione correctamente
options(na.action = "na.fail")

#Generar todos los subconjuntos de modelos posibles
MSVR_bic_et <- dredge(modelo_completo_et_red, rank = "BIC", fixed = variables_fijas)

# Extraer el mejor modelo según AICc
MSVR_bic_et <- get.models(MSVR_bic_et, 1)[[1]]
MSVR_bic_et

# Errores robustos
MSVR_bic_et_corr_hc <- coeftest(MSVR_bic_et, vcov. = vcovHC(MSVR_bic_et, "HC1", adjust = TRUE, sandwich = TRUE))

MSVR_bic_et_corr_hac <- coeftest(MSVR_bic_et, vcov. = vcovHAC(MSVR_bic_et, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(MSVR_bic_et, MSVR_bic_et_corr_hc,
          type = "html", 
          title = "Regresión con variables relacionadas al uso de tarjetas - Modelo BIC con efectos fijos de tiempo", 
          column.labels = c("MCO","MCO HC"), model.names = F,
          dep.var.labels = c("Tasa adquirencia", "Tasa adquirencia"), 
          dep.var.caption = "Variable dependiente")
```

###### Predicción

```{r, include=TRUE}
#Pronóstico para CR
pred_adq_MSVR_bic_et <- predict(MSVR_bic_et, newdata = data_CR, interval = "prediction", level = 0.95, type = "response")[which(datos$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MSVR_bic_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Regresión sin variables relacionadas al uso de tarjetas - Modelo BIC con efectos fijos de tiempo")
pronostico <- c(pred_adq_MSVR_bic_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MSVR_bic_et, modelo, pred_adq_MSVR_bic_et, pronostico, MSVR_bic_et_corr_hc, MSVR_bic_et_corr_hac, datos_red, modelo_completo_et_red, controles_red, formula_con_cont_t_red, variables_interes_red)
```

#### Enfoque2: modelo con variables instrumentales

::: {style="text-align: justify"}
Para el segundo enfoque, se plantea un modelo de variables instrumentales. En este caso, se instrumenta el valor promedio de la transacción realizada con tarjetas de pago con el Índice de Capital Humano. Este índice incorpora una medición del nivel educativo de un país, que es una característica del cliente que puede estar relacionada con su disposición a adoptar y utilizar tecnologías de pago distintas al efectivo. 

Dado que esta característica impacta el lado de la demanda, es posible utilizar el índice como un instrumento para trazar la oferta.
:::

```{r, include=FALSE}
#Datos sin variables relacionadas al uso de tarjetas
v_1 <- c("Pais", "Codigo_Pais", "Año", "Adquirencia", "Regulacion_intercambio")
v_2 <- c("Crimen", "Densidad_Pob", "Densidad_Neg", "PIB_pc_PPA", "Europa", "Asia_Este", "Am_lat", "Valor_prom")
v_3 <- c("Crimen", "Densidad_Pob", "Densidad_Neg", "PIB_pc_PPA", "Europa", "Asia_Este", "Am_lat", "ICH")
ef_tiempo_vi <- c("Año_2012", ef_tiempo)
vi <- c(v_1, v_2, "ICH", ef_tiempo_vi)

datos_vi <- datos %>%
  select(all_of(vi)) %>%
  na.omit()

# Fórmula sin efectos de tiempo
formula_vi <- as.formula(
  paste("Adquirencia ~ Regulacion_intercambio +", paste(v_2, collapse = " + "), "| Regulacion_intercambio + ", paste(v_3, collapse = " + ") )
)

# Fórmula con efectos de tiempo
formula_vi_et <- as.formula(
  paste("Adquirencia ~ Regulacion_intercambio + ", paste(c(v_2,ef_tiempo_vi), collapse = " + "), "| Regulacion_intercambio + ", paste(c(v_3, ef_tiempo_vi), collapse = " + ") )
)

```

##### Modelo sin efectos de tiempo

###### Estimación

```{r, include=TRUE, results='asis'}
#Regresión
modelo_vi <- ivreg(formula_vi,data = datos_vi)
modelo_vi

# Errores robustos
modelo_vi_corr_hc <- coeftest(modelo_vi, vcov. = vcovHC(modelo_vi, "HC1", adjust = TRUE, sandwich = TRUE))

modelo_vi_corr_hac <- coeftest(modelo_vi, vcov. = vcovHAC(modelo_vi, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(modelo_vi, modelo_vi_corr_hc, type = "html", 
          title = "Regresión con variables instrumentales - Sin efectos fijos de tiempo", 
          column.labels = c("Instrumental","Instrumental HC"), model.names = F, 
          dep.var.labels = c("Tasa adquirencia", "Tasa adquirencia"), 
          dep.var.caption = "Variable dependiente")
```  
 
###### Predicción 
 
```{r, include=TRUE}
#Pronósticos para Costa Rica
pred_adq_vi <- predict(modelo_vi, newdata = data_CR,
                                   interval = "prediction", level = 0.95,
                                   type = "response")[which(data_CR$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_vi`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Regresión con variables instrumentales - Sin efectos fijos de tiempo")
pronostico <- c(pred_adq_vi)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(modelo_vi, modelo, pred_adq_vi, pronostico, formula_vi, modelo_vi_corr_hc, modelo_vi_corr_hac)
```

##### Modelo con efectos de tiempo

###### Estimación

```{r, include=TRUE, results='asis'}
#Regresión
modelo_vi_et <- ivreg(formula_vi_et,data = datos_vi)
modelo_vi_et

# Errores robustos
modelo_vi_et_corr_hc <- coeftest(modelo_vi_et, vcov. = vcovHC(modelo_vi_et, "HC1", adjust = TRUE, sandwich = TRUE))

modelo_vi_et_corr_hac <- coeftest(modelo_vi_et, vcov. = vcovHAC(modelo_vi_et, adjust = TRUE, sandwich = TRUE))

# Visualización
stargazer(modelo_vi_et, modelo_vi_et_corr_hc, type = "html", 
          title = "Regresión con variables instrumentales - Con efectos fijos de tiempo", 
          column.labels = c("Instrumental","Instrumental HC"), model.names = F, 
          dep.var.labels = c("Tasa adquirencia", "Tasa adquirencia"), 
          dep.var.caption = "Variable dependiente")
```  
 
###### Predicción
 
```{r, include=TRUE}
#Pronósticos para Costa Rica
pred_adq_vi_et <- predict(modelo_vi_et, newdata = data_CR,
                                   interval = "prediction", level = 0.95,
                                   type = "response")[which(data_CR$Codigo_Pais=="CRI")]
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_vi_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Regresión con variables instrumentales - Con efectos fijos de tiempo")
pronostico <- c(pred_adq_vi_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(modelo_vi_et, modelo, pred_adq_vi_et, pronostico, formula_vi_et, modelo_vi_et_corr_hc, modelo_vi_et_corr_hac, datos_vi, ef_tiempo_vi, v_1, v_2, v_3, vi)
```

### Modelo de panel de datos con efectos aleatorios

::: {style="text-align: justify"}
La base de datos recopilada para este estudio contiene para algunos países el seguimiento del comportamiento de las comisiones de adquirencia e intercambio a lo largo del tiempo. Este seguimiento, sin embargo, es limitado. Para el 74% de los países se disponen como máximo de dos observaciones. No obstante, es posible considerar en el presente análisis de sensibilidad un modelo de datos de panel que permita explorar el papel de la heterogeneidad no observable en los resultados. 

Para este ejercicio, se utiliza un modelo de efectos aleatorios. Esto ya que Costa Rica no se incluye en la estimación del modelo, sino hasta que se realiza el pronóstico, lo que hace poco práctico un modelo de efectos fijos si el objetivo final del modelo es obtener valores esperados de las comisiones para Costa Rica. 
:::

```{r, include=FALSE}
# Crear el identificador de panel
datos_adq_con_cont_panel <- datos_adq_con_cont %>%
  mutate(id = as.numeric(factor(Codigo_Pais)))

#Declarar estructura de datos de panel
pdata <- pdata.frame(datos_adq_con_cont_panel, index = c("id", "Año"))

```

#### MPDEA_aicc

##### Estimación

```{r, include=TRUE}
# Modelo completo
modelo_completo_panel <- lm(formula_con_cont, data = datos_adq_con_cont_panel)

# Selección hacia atrás con AICc, fijando Regulacion_intercambio
MPDEA_aicc <- dredge(modelo_completo_panel, fixed = "Regulacion_intercambio", rank = "AICc")

# Extraer el mejor modelo
MPDEA_aicc <- get.models(MPDEA_aicc, 3)[[1]] #Ranking=3

# Ver variables seleccionadas
var_panel <- names(coef(MPDEA_aicc))[-1]  # excluye intercepto
var_panel <- var_panel[var_panel != "Regulacion_intercambio"]

# Fórmula para lmer
formula_lmer <- as.formula(paste("Adquirencia ~ Regulacion_intercambio +", paste(var_panel, collapse = " + "), " + (1 | Codigo_Pais)" ))

#Regresión con efectos aleatorios y errores agrupados
MPDEA_aicc <- lmer(formula_lmer, data = pdata)

resultado_panel_aicc <- summary(MPDEA_aicc)

summary(MPDEA_aicc)
```

##### Predicción

```{r, include=TRUE}
#Nombre de variables
var_finales <- attr(terms(MPDEA_aicc), "term.labels")

#Pronóstico
Caso_CR <- datos %>% filter(Codigo_Pais=="CRI") %>% dplyr::select(var_finales)
Caso_CR <- cbind(1,Caso_CR) # agregar un 1 para que se multiple con el intercepto general

pred_adq_MPDEA_aicc <- sum(resultado_panel_aicc$coefficients*Caso_CR)
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MPDEA_aicc`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación con efectos aleatorios - Modelo AICc sin efectos fijos de tiempo")
pronostico <- c(pred_adq_MPDEA_aicc)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MPDEA_aicc, modelo, pred_adq_MPDEA_aicc, pronostico, resultado_panel_aicc, formula_lmer, var_finales)
```

#### MPDEA_bic

##### Estimación

```{r, include=TRUE}
# Selección hacia atrás con AICc, fijando Regulacion_intercambio
MPDEA_bic <- dredge(modelo_completo_panel, fixed = "Regulacion_intercambio", rank = "BIC")

# Extraer el mejor modelo
MPDEA_bic <- get.models(MPDEA_bic, 20)[[1]] #Ranking=20

# Ver variables seleccionadas
var_panel <- names(coef(MPDEA_bic))[-1]  # excluye intercepto
var_panel <- var_panel[var_panel != "Regulacion_intercambio"]

# Fórmula para lmer
formula_lmer <- as.formula(paste("Adquirencia ~ Regulacion_intercambio +", paste(var_panel, collapse = " + "), " + (1 | Codigo_Pais)" ))

#Regresión con efectos aleatorios y errores agrupados
# Modelo de efectos aleatorios
MPDEA_bic <- lmer(formula_lmer, data = pdata)

# Resumen con errores agrupados por id
resultado_panel_bic <- summary(MPDEA_bic)

summary(MPDEA_bic)
```

##### Predicción

```{r, include=TRUE}
#Nombre de variables
var_finales <- attr(terms(MPDEA_bic), "term.labels")

#Pronóstico
Caso_CR <- datos %>% filter(Codigo_Pais=="CRI") %>% dplyr::select(var_finales)
Caso_CR <- cbind(1,Caso_CR) # agregar un 1 para que se multiple con el intercepto general

pred_adq_MPDEA_bic <- sum(resultado_panel_bic$coefficients*Caso_CR)
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MPDEA_bic`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación con efectos aleatorios - Modelo BIC sin efectos fijos de tiempo")
pronostico <- c(pred_adq_MPDEA_bic)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MPDEA_bic, modelo, pred_adq_MPDEA_bic, pronostico, resultado_panel_bic, formula_lmer, var_finales, modelo_completo_panel)
```

#### MPDEA_aicc_et **

##### Estimación

```{r, include=TRUE}
# Modelo completo
modelo_completo_panel_et <- lm(formula_con_cont_t, data = datos_adq_con_cont_panel)

# Selección hacia atrás con AICc, fijando Regulacion_intercambio
MPDEA_aicc_et <- dredge(modelo_completo_panel_et, fixed = variables_fijas, rank = "AICc")

# Extraer el mejor modelo
MPDEA_aicc_et <- get.models(MPDEA_aicc_et, 11)[[1]] #Ranking=11

# Ver variables seleccionadas
coef_original <- coef(MPDEA_aicc_et)
coef_sin_na <- coef_original[!is.na(coef_original)]
var_panel <- names(coef_sin_na)
var_panel <- var_panel[var_panel != "(Intercept)"] # excluye intercepto
var_panel <- var_panel[var_panel != "Regulacion_intercambio"] 

# Fórmula para lmer
formula_lmer <- as.formula(paste("Adquirencia ~ Regulacion_intercambio +", paste(var_panel, collapse = " + "), " + (1 | Codigo_Pais)" ))

#Regresión con efectos aleatorios y errores agrupados
MPDEA_aicc_et <- lmer(formula_lmer, data = pdata)

# Resumen con errores agrupados por id
resultado_panel_aicc_et <- summary(MPDEA_aicc_et)

summary(MPDEA_aicc_et)
```

##### Predicción

```{r, include=TRUE}
#Nombre de variables
var_finales <- attr(terms(MPDEA_aicc_et), "term.labels")

#Pronóstico
Caso_CR <- datos %>% filter(Codigo_Pais=="CRI") %>% dplyr::select(var_finales)
Caso_CR<-cbind(1,Caso_CR) # agregar un 1 para que se multiple con el intercepto general

pred_adq_MPDEA_aicc_et <- sum(resultado_panel_aicc_et$coefficients*Caso_CR)
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MPDEA_aicc_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación con efectos aleatorios - Modelo AICc con efectos fijos de tiempo")
pronostico <- c(pred_adq_MPDEA_aicc_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MPDEA_aicc_et, modelo, pred_adq_MPDEA_aicc_et, pronostico, resultado_panel_aicc_et, formula_lmer, var_finales, coef_original, coef_sin_na, var_panel)
```

#### MPDEA_bic_et **

##### Estimación

```{r, include=TRUE}
# Selección hacia atrás con AICc, fijando Regulacion_intercambio
MPDEA_bic_et <- dredge(modelo_completo_panel_et, fixed = variables_fijas, rank = "BIC")

# Extraer el mejor modelo
MPDEA_bic_et <- get.models(MPDEA_bic_et, 10)[[1]] #Ranking=10

# Ver variables seleccionadas
coef_original <- coef(MPDEA_bic_et)
coef_sin_na <- coef_original[!is.na(coef_original)]
var_panel <- names(coef_sin_na)
var_panel <- var_panel[var_panel != "(Intercept)"] # excluye intercepto
var_panel <- var_panel[var_panel != "Regulacion_intercambio"] 

# Fórmula para lmer
formula_lmer <- as.formula(paste("Adquirencia ~ Regulacion_intercambio +", paste(var_panel, collapse = " + "), " + (1 | Codigo_Pais)" ))

#Regresión con efectos aleatorios y errores agrupados
MPDEA_bic_et <- lmer(formula_lmer, data = pdata)

# Resumen con errores agrupados por id
resultado_panel_bic_et <- summary(MPDEA_bic_et)

summary(MPDEA_bic_et)
```

##### Predicción

```{r, include=TRUE}
#Nombre de variables
var_finales <- attr(terms(MPDEA_bic_et), "term.labels")

#Pronóstico
Caso_CR <- datos %>% filter(Codigo_Pais=="CRI") %>% dplyr::select(var_finales)
Caso_CR<-cbind(1,Caso_CR) # agregar un 1 para que se multiple con el intercepto general

pred_adq_MPDEA_bic_et <- sum(resultado_panel_bic_et$coefficients*Caso_CR)
```

::: {style="text-align: justify"}
Predicción tasa de adquirencia para Costa Rica: **`r pred_adq_MPDEA_bic_et`**.
:::

```{r, include=FALSE}
#Actualiza el data frame para guardar la predicción
modelo <- c("Estimación con efectos aleatorios - Modelo BIC con efectos fijos de tiempo")
pronostico <- c(pred_adq_MPDEA_bic_et)  
Pronosticos <- rbind(Pronosticos, data.frame(Modelo = modelo, Pronostico = pronostico))
```

```{r, include=FALSE}
rm(MPDEA_bic_et, modelo, pred_adq_MPDEA_bic_et, pronostico, resultado_panel_bic_et, formula_lmer, var_finales, coef_original, coef_sin_na, var_panel, datos_adq_con_cont_panel, modelo_completo_panel_et, pdata, Caso_CR)
```

## Estimación de las comisiones de adquirencia

::: {style="text-align: justify"}
Posterior a la etapa de modelación, se realiza un pronóstico para conocer la tasa de adquirencia esperada para Costa Rica. A continuación, se calcula un promedio para los valores esperados de cada modelo estimado.
:::

```{r, include=TRUE}
Pronosticos$Tipo_modelo <- Pronosticos$Modelo

Pronosticos <- Pronosticos %>%
  separate(Tipo_modelo, into = c("Tipo_modelo", "Especificacion"), sep = "-")

Pronosticos$Pronostico_ajustado <- (Pronosticos$Pronostico - (1*0.09) - (1.25*0.06) - (1.5*0.07) - (1.75*0.15))/(0.63)


Pronosticos_resumen <- Pronosticos[-1, ] %>%
  group_by(Tipo_modelo) %>%
  summarise(promedio = mean(Pronostico_ajustado, na.rm = TRUE))

print(Pronosticos_resumen)

promedio_adquirencia_CR <- Pronosticos_resumen %>%
                            filter(Tipo_modelo == "Estimación principal ") %>%
                            pull(promedio)

Pronosticos_resumen_adquirencia <- Pronosticos_resumen

```

::: {style="text-align: justify"}
En promedio, la tasa de adquirencia estimada para Costa Rica es: **`r promedio_adquirencia_CR`**.
:::

