Definición de funciones
#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)))
}
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)
)
}
Resultados
bind_rows(Resultados_Performance_data_entrenamiento) %>%
stargazer(title = "Medidas de Datos del Modelo",type = "html",digits = 3,summary.stat = c("n", "mean", "sd", "min", "p25", "p75","max"))
Medidas de 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
|
|
|