Actualización del script de simulación

Modelo Original

options(scipen = 99999)
load("C:/Users/Economia/Desktop/smoke.RData") # en otros ejercicios podría ser necesario importar un csv o un archivo de Excel
data<-data #aquí se puede cambiar el dataframe de manera sencilla
equation<-as.formula("cigs~cigpric+lcigpric+income+lincome+age+agesq+educ+white+restaurn") #Guardar la fórmula del modelo para no tener que estar escribiendola en cada ocasión
endogena<-data$cigs #aquí designamos la variable endógena del modelo, esto es necesario para el script de simulación.
Modelo_smoke<-lm(formula =equation,data = data)
stargazer::stargazer(Modelo_smoke,
                     title = "Modelo Estimado",
                     type = "html",
                     digits = 8)

Modelo Estimado

	Dependent variable:
	cigs
cigpric	2.00226800
	(1.49283100)
lcigpric	-115.27350000
	(85.42432000)
income	-0.00004619
	(0.00013349)
lincome	1.40406100
	(1.70816600)
age	0.77835900***
	(0.16055560)
agesq	-0.00915035***
	(0.00174929)
educ	-0.49478060***
	(0.16818020)
white	-0.53105160
	(1.46072200)
restaurn	-2.64424100**
	(1.12999900)
Constant	340.80440000
	(260.01560000)
Observations	807
R2	0.05515398
Adjusted R2	0.04448445
Residual Std. Error	13.41285000 (df = 797)
F Statistic	5.16929800*** (df = 9; 797)
Note:	p<0.1; p<0.05; p<0.01

#Funciones para la desigualdad de 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) 
library(caret)
library(DescTools) 
library(stargazer) 
set.seed(50) 
numero_de_muestras<-500
proporcion_entrenamiento<-0.75
# Creación de las muestras, aquí usamos la variable endógena que definimos con anterioridad 
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)
# Estimar los modelos de cada muestra y sus medidas de desempeño predictivo
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)
            )
}

Desempeño con los datos de entrenamiento

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

Desempeño con los datos de prueba (los simulados)

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