Estimadores HAC

UNIVERSIDAD DE EL SALVADOR

FACULTAD DE CIENCIAS ECONÓMICAS

ESCUELA DE ECONOMÍA

CICLO I - 2023

Asignatura

Econometría

Grupo teórico

GT-03

Docente

MSF.Carlos Ademir Peréz

Estudiante:

Ascencio Orellana, Fernando Javier AO21017

Ciudad universitaria, 23 de junio de 2023

Cargar los datos

options(scipen = 999999)
library(lmtest)
library(stargazer)
library(equatiomatic) 
library(wooldridge)
data <- wooldridge::smoke
equation<-as.formula("cigs~cigpric+lcigpric+income+lincome+age+agesq+educ+white+restaurn") 
endogena<-data$cigs 
Modelo_estimado<-lm(formula =equation,data = data)
stargazer::stargazer(Modelo_estimado,
                     title = "Modelo Estimado",
                     type = "html",
                     digits = 7)
Modelo Estimado
Dependent variable:
cigs
cigpric 2.0022680
(1.4928310)
lcigpric -115.2735000
(85.4243200)
income -0.0000462
(0.0001335)
lincome 1.4040610
(1.7081660)
age 0.7783590***
(0.1605556)
agesq -0.0091504***
(0.0017493)
educ -0.4947806***
(0.1681802)
white -0.5310516
(1.4607220)
restaurn -2.6442410**
(1.1299990)
Constant 340.8044000
(260.0156000)
Observations 807
R2 0.0551540
Adjusted R2 0.0444844
Residual Std. Error 13.4128500 (df = 797)
F Statistic 5.1692980*** (df = 9; 797)
Note: p<0.1; p<0.05; p<0.01 
extract_eq(Modelo_estimado,wrap = TRUE)

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

Funciones para la desigualdad Theil

Um<-function(pronosticado,observado){
  library(DescTools)
  ((mean(pronosticado)-mean(observado))^2)/MSE(pronosticado,observado) 
}
Us<-function(pronosticado,observado){
  library(DescTools)
  ((sd(pronosticado)-sd(observado))^2)/MSE(pronosticado,observado)
}
Uc<-function(pronosticado,observado){
  library(DescTools)
  (2*(1-cor(pronosticado,observado))*sd(pronosticado)*sd(observado))/MSE(pronosticado,observado)}
THEIL_U<-function(pronosticado,observado){
   library(DescTools)
  RMSE(pronosticado,observado)/(sqrt(mean(pronosticado^2))+sqrt(mean(observado^2)))
}

Simulación

#Script de Simulación versión 2
options(scipen = 999999) 
library(dplyr) 
## 
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(caret)
## Loading required package: ggplot2
## Loading required package: lattice
library(DescTools) 
## 
## Attaching package: 'DescTools'
## The following objects are masked from 'package:caret':
## 
##     MAE, RMSE
library(stargazer) 
set.seed(50) 
numero_de_muestras<-500
proporcion_entrenamiento<-0.75
muestras<- endogena %>%
  createDataPartition(p = proporcion_entrenamiento,
                      times = numero_de_muestras,
                      list = TRUE)

#Listas vacias para la simulación
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)

for(j in 1:numero_de_muestras){
Datos_Entrenamiento<- data[muestras[[j]], ]
Datos_Prueba<- data[-muestras[[j]], ]
Modelos_Entrenamiento[[j]]<-lm(formula = equation,data=Datos_Entrenamiento)
Pronostico_Prueba[[j]]<-Modelos_Entrenamiento[[j]] %>% predict(Datos_Prueba)
Fe<-Modelos_Entrenamiento[[j]]$fitted.values
Ye<-Datos_Entrenamiento$cigs
Resultados_Performance_data_entrenamiento[[j]]<-data.frame( 
            R2 = R2(Fe,Ye),
            RMSE = RMSE(Fe,Ye),
            MAE = MAE(Fe,Ye),
            MAPE= MAPE(Fe,Ye)*100,
            THEIL=TheilU(Fe,Ye,type = 1),
            Um=Um(Fe,Ye),
            Us=Us(Fe,Ye),
            Uc=Uc(Fe,Ye)
            )
Fp<-Pronostico_Prueba[[j]]
Yp<-Datos_Prueba$cigs
Resultados_Performance[[j]]<-data.frame( 
            R2 = R2(Fp,Yp ),
            RMSE = RMSE(Fp, Yp),
            MAE = MAE(Fp,Yp),
            MAPE= MAPE(Fp,Yp)*100,
            THEIL=TheilU(Fp,Yp,type = 1), 
            Um=Um(Fp,Yp),
            Us=Us(Fp,Yp),
            Uc=Uc(Fp,Yp)
            )
}

Resultados de la simulación

library(dplyr)
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