SIMULACIÓN

UNIVERSIDAD DE EL SALVADOR
FACULTAD DE CIENCIAS ECONÓMICAS
ESCUELA DE ECONOMÍA
CICLO I - 2023
“SIMULACIÓN”
ASIGNATURA: ECONOMETRÍA
CATEDRÁTICO: MSF. CARLOS ADEMIR PERÉZ
GRUPO TEÓRICO: 03
PRESENTADO POR: AGUILAR ZACARIAS, CRISTALI DAYAMARI AZ20006

Parte 1

SIMULACIÓN DE LA PRESENTACIÓN

1. Importación de datos

options(scipen = 999999)
library(lmtest)
library(stargazer)
library(equatiomatic)
library(mlbench) # Libreria con data frame BostonHousing
data("BostonHousing")
modelo_boston<-lm(formula = medv~.,data=BostonHousing)
extract_eq(modelo_boston,wrap = TRUE)

\[ \begin{aligned} \operatorname{medv} &= \alpha + \beta_{1}(\operatorname{crim}) + \beta_{2}(\operatorname{zn}) + \beta_{3}(\operatorname{indus})\ + \\ &\quad \beta_{4}(\operatorname{chas}_{\operatorname{1}}) + \beta_{5}(\operatorname{nox}) + \beta_{6}(\operatorname{rm}) + \beta_{7}(\operatorname{age})\ + \\ &\quad \beta_{8}(\operatorname{dis}) + \beta_{9}(\operatorname{rad}) + \beta_{10}(\operatorname{tax}) + \beta_{11}(\operatorname{ptratio})\ + \\ &\quad \beta_{12}(\operatorname{b}) + \beta_{13}(\operatorname{lstat}) + \epsilon \end{aligned} \]

coeftest(modelo_boston)
## 
## t test of coefficients:
## 
##                 Estimate   Std. Error  t value              Pr(>|t|)    
## (Intercept)  36.45948839   5.10345881   7.1441    0.0000000000032834 ***
## crim         -0.10801136   0.03286499  -3.2865             0.0010868 ** 
## zn            0.04642046   0.01372746   3.3816             0.0007781 ***
## indus         0.02055863   0.06149569   0.3343             0.7382881    
## chas1         2.68673382   0.86157976   3.1184             0.0019250 ** 
## nox         -17.76661123   3.81974371  -4.6513    0.0000042456438076 ***
## rm            3.80986521   0.41792525   9.1161 < 0.00000000000000022 ***
## age           0.00069222   0.01320978   0.0524             0.9582293    
## dis          -1.47556685   0.19945473  -7.3980    0.0000000000006013 ***
## rad           0.30604948   0.06634644   4.6129    0.0000050705290227 ***
## tax          -0.01233459   0.00376054  -3.2800             0.0011116 ** 
## ptratio      -0.95274723   0.13082676  -7.2825    0.0000000000013088 ***
## b             0.00931168   0.00268596   3.4668             0.0005729 ***
## lstat        -0.52475838   0.05071528 -10.3471 < 0.00000000000000022 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2. SIMULACIÓN

2.1. Definición de funciones para el cálculo de Theil y su descomposición

#Bias Proportion
Um <-function(pronosticado,observado){
  library(DescTools)
  ((mean(pronosticado)-mean(observado))^2)/MSE(pronosticado,observado) 
}
#Variance Proportion
Us <-function(pronosticado,observado){
  library(DescTools)
  ((sd(pronosticado)-sd(observado))^2)/MSE(pronosticado,observado)
}
#Covariance Proportion
Uc <-function(pronosticado,observado){
  library(DescTools)
  (2*(1-cor(pronosticado,observado))*sd(pronosticado)*sd(observado))/MSE(pronosticado,observado)}

#Coeficiente U de Theil (también aparece en la librería "DescTools")
THEIL_U <-function(pronosticado,observado){
   library(DescTools)
  RMSE(pronosticado,observado)/(sqrt(mean(pronosticado^2))+sqrt(mean(observado^2)))
}

2.2. Script de simulación

options(scipen = 999999)
library(dplyr) 
library(caret) 
library(DescTools) 
library(stargazer)
set.seed(50) 
numero_de_muestras<-1000 

# Se crea la lista con las 1000 muestras (indica la posición de la fila en cada data frame)
muestras<- BostonHousing$medv %>%
  createDataPartition(p = 0.8,
                      times = numero_de_muestras,
                      list = TRUE)

# Listas vacias, que contendran los datos de entrenamiento, los pronosticos para los datos de prueba, y para las estadisticas de cada muestra
Modelos_Entrenamiento<-vector(mode = "list",
                              length = numero_de_muestras)
Pronostico_Prueba<-vector(mode = "list",
                              length = numero_de_muestras)
Resultados_Performance_data_entrenamiento<-vector(mode = "list",
                              length = numero_de_muestras)

Resultados_Performance<-vector(mode = "list",
                              length = numero_de_muestras)

#Estimación de los modelos lineales para cada muestra, los pronósticos y cálculo de las estadisticas de performance.
for(j in 1:numero_de_muestras){
Datos_Entrenamiento<- BostonHousing[muestras[[j]], ]
Datos_Prueba<- BostonHousing[-muestras[[j]], ]
Modelos_Entrenamiento[[j]]<-lm(formula = medv~.,data=Datos_Entrenamiento)
Pronostico_Prueba[[j]]<-Modelos_Entrenamiento[[j]] %>% predict(Datos_Prueba)
Resultados_Performance_data_entrenamiento[[j]]<-data.frame( 
           R2 = R2(Modelos_Entrenamiento[[j]]$fitted.values,
                    Datos_Entrenamiento$medv),
            RMSE = RMSE(Modelos_Entrenamiento[[j]]$fitted.values,
                        Datos_Entrenamiento$medv),
            MAE = MAE(Modelos_Entrenamiento[[j]]$fitted.values,
                      Datos_Entrenamiento$medv),
            MAPE= MAPE(Modelos_Entrenamiento[[j]]$fitted.values,
                       Datos_Entrenamiento$medv)*100,
            THEIL =TheilU(Modelos_Entrenamiento[[j]]$fitted.values,
                         Datos_Entrenamiento$medv,type = 1),
            Um =Um(Modelos_Entrenamiento[[j]]$fitted.values,
                         Datos_Entrenamiento$medv),
            Us=Us(Modelos_Entrenamiento[[j]]$fitted.values,
                         Datos_Entrenamiento$medv),
            Uc=Uc(Modelos_Entrenamiento[[j]]$fitted.values,
                         Datos_Entrenamiento$medv)
            )
Resultados_Performance[[j]]<-data.frame( 
            R2 = R2(Pronostico_Prueba[[j]], Datos_Prueba$medv),
            RMSE = RMSE(Pronostico_Prueba[[j]], Datos_Prueba$medv),
            MAE = MAE(Pronostico_Prueba[[j]], Datos_Prueba$medv),
            MAPE= MAPE(Pronostico_Prueba[[j]], Datos_Prueba$medv)*100,
            THEIL=TheilU(Pronostico_Prueba[[j]], Datos_Prueba$medv,
                         type = 1), 
            Um=Um(Pronostico_Prueba[[j]], Datos_Prueba$medv),
            Us=Us(Pronostico_Prueba[[j]], Datos_Prueba$medv),
            Uc=Uc(Pronostico_Prueba[[j]], Datos_Prueba$medv)
            )
}

2.3. Resultados de la simulación

bind_rows(Resultados_Performance_data_entrenamiento) %>% 
  stargazer(title = "Medidas de Performance Datos del Modelo",
            type = "html",
            digits = 3,
            summary.stat = c("n", "mean", "sd", "min", "p25", "p75", "max"))
Medidas de Performance Datos del Modelo
Statistic N Mean St. Dev. Min Pctl(25) Pctl(75) Max
R2 1,000 0.743 0.013 0.713 0.734 0.751 0.794
RMSE 1,000 4.653 0.141 4.177 4.565 4.759 4.948
MAE 1,000 3.265 0.095 2.905 3.204 3.332 3.512
MAPE 1,000 16.387 0.464 14.813 16.085 16.718 17.691
THEIL 1,000 0.096 0.003 0.087 0.095 0.099 0.102
Um 1,000 0.000 0.000 0 0 0 0
Us 1,000 0.074 0.004 0.058 0.072 0.077 0.085
Uc 1,000 0.928 0.004 0.918 0.925 0.931 0.945 
bind_rows(Resultados_Performance) %>% 
  stargazer(title = "Medidas de Performance Simulación",
            type = "html",
            digits = 3,
            summary.stat = c("n", "mean", "sd", "min", "p25", "p75", "max"))
Medidas de Performance Simulación
Statistic N Mean St. Dev. Min Pctl(25) Pctl(75) Max
R2 1,000 0.723 0.056 0.452 0.690 0.764 0.840
RMSE 1,000 4.862 0.575 3.465 4.450 5.226 6.961
MAE 1,000 3.411 0.281 2.633 3.216 3.598 4.492
MAPE 1,000 17.197 1.618 12.875 16.066 18.262 23.137
THEIL 1,000 0.101 0.012 0.073 0.092 0.108 0.148
Um 1,000 0.011 0.016 0.000 0.001 0.015 0.205
Us 1,000 0.081 0.066 0.00000 0.027 0.122 0.333
Uc 1,000 0.918 0.066 0.667 0.875 0.971 1.010

Parte 2

Simulación del pronóstico con los datos del ejercicio de estimadores HAC

1. Importación de datos

options(scipen = 999999)
library(lmtest)
library(stargazer)
library(equatiomatic)
load("C:/Users/crist/Downloads/Cristali Dayamari Aguilar Zacarias - smoke.RData")
data(data)
modelo_smoke<-lm(cigs~ cigpric + lcigpric + income + lincome + age +agesq + educ + white + restaurn, data = data)
extract_eq(modelo_smoke,wrap = TRUE)

\[ \begin{aligned} \operatorname{cigs} &= \alpha + \beta_{1}(\operatorname{cigpric}) + \beta_{2}(\operatorname{lcigpric}) + \beta_{3}(\operatorname{income})\ + \\ &\quad \beta_{4}(\operatorname{lincome}) + \beta_{5}(\operatorname{age}) + \beta_{6}(\operatorname{agesq}) + \beta_{7}(\operatorname{educ})\ + \\ &\quad \beta_{8}(\operatorname{white}) + \beta_{9}(\operatorname{restaurn}) + \epsilon \end{aligned} \]

coeftest(modelo_smoke)
## 
## t test of coefficients:
## 
##                   Estimate     Std. Error t value     Pr(>|t|)    
## (Intercept)  340.804374604  260.015587269  1.3107     0.190334    
## cigpric        2.002267667    1.492831189  1.3413     0.180220    
## lcigpric    -115.273464445   85.424315195 -1.3494     0.177585    
## income        -0.000046194    0.000133491 -0.3460     0.729402    
## lincome        1.404061178    1.708165841  0.8220     0.411340    
## age            0.778359013    0.160555612  4.8479 0.0000015001 ***
## agesq         -0.009150353    0.001749292 -5.2309 0.0000002158 ***
## educ          -0.494780616    0.168180198 -2.9420     0.003356 ** 
## white         -0.531051635    1.460721806 -0.3636     0.716287    
## restaurn      -2.644241351    1.129998690 -2.3400     0.019528 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

2. SIMULACIÓN

2.1. Definición de funciones para el cálculo de Theil y su descomposición

#Bias Proportion
Um <-function(pronosticado,observado){
  library(DescTools)
  ((mean(pronosticado)-mean(observado))^2)/MSE(pronosticado,observado) 
}
#Variance Proportion
Us <-function(pronosticado,observado){
  library(DescTools)
  ((sd(pronosticado)-sd(observado))^2)/MSE(pronosticado,observado)
}
#Covariance Proportion
Uc <-function(pronosticado,observado){
  library(DescTools)
  (2*(1-cor(pronosticado,observado))*sd(pronosticado)*sd(observado))/MSE(pronosticado,observado)}

#Coeficiente U de Theil (también aparece en la librería "DescTools")
THEIL_U <-function(pronosticado,observado){
   library(DescTools)
  RMSE(pronosticado,observado)/(sqrt(mean(pronosticado^2))+sqrt(mean(observado^2)))
}

2.2. Script de simulación

options(scipen = 999999)
library(dplyr) 
library(caret) 
library(DescTools) 
library(stargazer)
set.seed(50) 
numero_de_muestras<-500 

# Se crea la lista con las 500 muestras (indica la posición de la fila en cada data frame)
muestras<- data$cigs %>%
  createDataPartition(p = 0.75,
                      times = numero_de_muestras,
                      list = TRUE)

# Listas vacias, que contendran los datos de entrenamiento, los pronosticos para los datos de prueba, y para las estadisticas de cada muestra
Modelos_Entrenamiento<-vector(mode = "list",
                              length = numero_de_muestras)
Pronostico_Prueba<-vector(mode = "list",
                              length = numero_de_muestras)
Resultados_Performance_data_entrenamiento<-vector(mode = "list",
                              length = numero_de_muestras)
Resultados_Performance<-vector(mode = "list",
                              length = numero_de_muestras)

#Estimación de los modelos lineales para cada muestra, los pronósticos y cálculo de las estadisticas de performance.
for(j in 1:numero_de_muestras){
Datos_Entrenamiento<- data[muestras[[j]], ]
Datos_Prueba<- data[-muestras[[j]], ]
Modelos_Entrenamiento[[j]]<-lm(formula = cigs~ cigpric + lcigpric + income + lincome + age +agesq + educ + white + restaurn, data = Datos_Entrenamiento)
Pronostico_Prueba[[j]]<-Modelos_Entrenamiento[[j]] %>% predict(Datos_Prueba)
Resultados_Performance_data_entrenamiento[[j]]<-data.frame( 
           R2 = R2(Modelos_Entrenamiento[[j]]$fitted.values,
                    Datos_Entrenamiento$cigs),
            RMSE = RMSE(Modelos_Entrenamiento[[j]]$fitted.values,
                        Datos_Entrenamiento$cigs),
            MAE = MAE(Modelos_Entrenamiento[[j]]$fitted.values,
                      Datos_Entrenamiento$cigs),
            MAPE= MAPE(Modelos_Entrenamiento[[j]]$fitted.values,
                       Datos_Entrenamiento$cigs)*100,
            THEIL =TheilU(Modelos_Entrenamiento[[j]]$fitted.values,
                         Datos_Entrenamiento$cigs,type = 1),
            Um =Um(Modelos_Entrenamiento[[j]]$fitted.values,
                         Datos_Entrenamiento$cigs),
            Us=Us(Modelos_Entrenamiento[[j]]$fitted.values,
                         Datos_Entrenamiento$cigs),
            Uc=Uc(Modelos_Entrenamiento[[j]]$fitted.values,
                         Datos_Entrenamiento$cigs)
            )
Resultados_Performance[[j]]<-data.frame( 
            R2 = R2(Pronostico_Prueba[[j]], Datos_Prueba$cigs),
            RMSE = RMSE(Pronostico_Prueba[[j]], Datos_Prueba$cigs),
            MAE = MAE(Pronostico_Prueba[[j]], Datos_Prueba$cigs),
            MAPE= MAPE(Pronostico_Prueba[[j]], Datos_Prueba$cigs)*100,
            THEIL=TheilU(Pronostico_Prueba[[j]], Datos_Prueba$cigs,
                         type = 1), 
# Creamos THIL_U
            Um=Um(Pronostico_Prueba[[j]], Datos_Prueba$cigs),
            Us=Us(Pronostico_Prueba[[j]], Datos_Prueba$cigs),
            Uc=Uc(Pronostico_Prueba[[j]], Datos_Prueba$cigs)
            )
}

2.3. Resultados de la simulación

bind_rows(Resultados_Performance_data_entrenamiento) %>% 
  stargazer(title = "Medidas de Performance Datos del Modelo",
            type = "html",
            digits = 3,
            summary.stat = c("n", "mean", "sd", "min", "p25", "p75", "max"))
Medidas de Performance Datos del Modelo
Statistic N Mean St. Dev. Min Pctl(25) Pctl(75) Max
R2 500 0.058 0.007 0.040 0.054 0.063 0.083
RMSE 500 13.304 0.198 12.629 13.175 13.439 13.828
MAE 500 10.562 0.133 10.068 10.474 10.654 10.969
MAPE 500 Inf.000 Inf Inf Inf Inf
THEIL 500 0.522 0.006 0.505 0.517 0.526 0.541
Um 500 0.000 0.000 0 0 0 0
Us 500 0.613 0.020 0.554 0.601 0.625 0.669
Uc 500 0.389 0.020 0.332 0.376 0.401 0.448 
bind_rows(Resultados_Performance) %>% 
  stargazer(title = "Medidas de Performance Simulación",
            type = "html",
            digits = 3,
            summary.stat = c("n", "mean", "sd", "min", "p25", "p75", "max"))
Medidas de Performance Simulación
Statistic N Mean St. Dev. Min Pctl(25) Pctl(75) Max
R2 500 0.039 0.018 0.002 0.027 0.051 0.099
RMSE 500 13.479 0.595 11.734 13.064 13.853 15.361
MAE 500 10.728 0.290 9.834 10.528 10.904 11.626
MAPE 500 Inf.000 Inf Inf Inf Inf
THEIL 500 0.528 0.013 0.491 0.519 0.536 0.564
Um 500 0.002 0.003 0.00000 0.0003 0.003 0.021
Us 500 0.597 0.051 0.476 0.563 0.636 0.751
Uc 500 0.405 0.051 0.253 0.367 0.441 0.529