Ejercicio de Pronóstico y Simulación
Parte 1.
Ejercicio de Pronóstico y Simulación
Simulación de la presentación.
Estimación del modelo.
options(scipen = 999999)
library(lmtest)
library(stargazer)
library(equatiomatic) # optativo remotes::install_github("datalorax/equatiomatic")
library(mlbench) #esta librería tiene el data frame BostonHousing
data(BostonHousing)
#Modelo estimado medv~. indica "medv" en función del resto de variables del dataframe
modelo_boston<-lm(formula = medv~.,data=BostonHousing)
extract_eq(modelo_boston,wrap = TRUE) #optativo\[ \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} \]
##
## 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 ' ' 1Simulación en R.
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)))
}Script de Simulación.
options(scipen = 999999)
library(dplyr)
library(caret) # Permite Realizar muestreo sobre los data frame
library(DescTools) # Contiene las funciones para calcular las medidas de performance
library(stargazer) # Para dar formato, y obtener resumen estadistico de las simulaciones
set.seed(50) # Permite fijar la semilla aleatoria, para reproducir los resultados obtenidos en esta clase
numero_de_muestras<-1000 # Numero de muestras que se optendran del data frame
# 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), # También se puede usar la función que creamos: THEIL_U
Um=Um(Pronostico_Prueba[[j]], Datos_Prueba$medv),
Us=Us(Pronostico_Prueba[[j]], Datos_Prueba$medv),
Uc=Uc(Pronostico_Prueba[[j]], Datos_Prueba$medv)
)
} #No olvidar este corchete ;)Resultados de la Simulación.
bind_rows(Resultados_Performance_data_entrenamiento) %>%
stargazer(title = "Medidas de Performance Datos del Modelo",
type = "text",
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 Simulacion",
type = "text",
digits = 3,
summary.stat = c("n", "mean", "sd", "min", "p25", "p75", "max"))##
## Medidas de Performance Simulacion
## ================================================================
## 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 con los datos del ejercicio de estimadores HAC.
Estimación del modelo.
options(scipen = 999999)
library(stargazer)
load("Salvador Antonio Figueroa Gonzalez - smoke.RData")
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 = "text",
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.01Simulación en R.
Definición de funciones para el cálculo de Theil y su descomposición
#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)))
}Script de Simulación.
options(scipen = 999999) #No mostrar notación cientifica.
library(dplyr) # Para manejo de datos y activar el operador "pipe" %>%
library(caret) # Permite Realizar muestreo sobre los data frame
library(DescTools)
library(stargazer)
set.seed(50)
numero_de_muestras_2<-500
proporcion_entrenamiento_2<-0.75
# Creación de las muestras, aquí usamos la variable endógena que definimos con anterioridad
muestras_2<- endogena %>%
createDataPartition(p = proporcion_entrenamiento_2,
times = numero_de_muestras_2,
list = TRUE)
#Listas vacias para la simulación
Modelos_Entrenamiento_2<-vector(mode = "list",
length = numero_de_muestras_2)
Pronostico_Prueba_2<-vector(mode = "list",
length = numero_de_muestras_2)
Resultados_Performance_data_entrenamiento_2<-vector(mode = "list",
length = numero_de_muestras_2)
Resultados_Performance_2<-vector(mode = "list",
length = numero_de_muestras_2)
# Estimar los modelos de cada muestra y sus medidas de desempeño predictivo
for(j in 1:numero_de_muestras_2){
Datos_Entrenamiento_2<- data[muestras_2[[j]], ]
Datos_Prueba_2<- data[-muestras_2[[j]], ]
Modelos_Entrenamiento_2[[j]]<-lm(formula = equation,data=Datos_Entrenamiento_2)
Pronostico_Prueba_2[[j]]<-Modelos_Entrenamiento_2[[j]] %>% predict(Datos_Prueba_2)
Fe<-Modelos_Entrenamiento_2[[j]]$fitted.values
Ye<-Datos_Entrenamiento_2$cigs
Resultados_Performance_data_entrenamiento_2[[j]]<-data.frame(
R2_2 = R2(Fe,Ye),
RMSE_2 = RMSE(Fe,Ye),
MAE_2 = MAE(Fe,Ye),
MAPE_2= MAPE(Fe,Ye)*100,
THEIL_2=TheilU(Fe,Ye,type = 1),
Um_2=Um(Fe,Ye),
Us_2=Us(Fe,Ye),
Uc_2=Uc(Fe,Ye)
)
Fp<-Pronostico_Prueba_2[[j]]
Yp<-Datos_Prueba_2$cigs
Resultados_Performance_2[[j]]<-data.frame(
R2_2 = R2(Fp,Yp ),
RMSE_2 = RMSE(Fp, Yp),
MAE_2 = MAE(Fp,Yp),
MAPE_2= MAPE(Fp,Yp)*100,
THEIL_2=TheilU(Fp,Yp,type = 1),
Um_2=Um(Fp,Yp),
Us_2=Us(Fp,Yp),
Uc_2=Uc(Fp,Yp)
)
}Resultados de la Simulación.
library(dplyr)
bind_rows(Resultados_Performance_data_entrenamiento_2) %>%
stargazer(title= "Medidas de Performance Datos del Modelo",
type = "text",
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_2 500 0.058 0.007 0.040 0.054 0.063 0.083
## RMSE_2 500 13.304 0.198 12.629 13.175 13.439 13.828
## MAE_2 500 10.562 0.133 10.068 10.474 10.654 10.969
## MAPE_2 500 Inf.000 Inf Inf Inf Inf
## THEIL_2 500 0.522 0.006 0.505 0.517 0.526 0.541
## Um_2 500 0.000 0.000 0 0 0 0
## Us_2 500 0.613 0.020 0.554 0.601 0.625 0.669
## Uc_2 500 0.389 0.020 0.332 0.376 0.401 0.448
## --------------------------------------------------------------Para el valor del R2, en promedio el modelo explica el 5.8% de la varianza de la variable endógena ante cambios en las variables exógenas, como mínimo se explica en 0.7% de la varianza y como máximo se puede esperar un poder explicativo del 8.3% de la varianza de la variable endógena. Para el RMSE en promedio hay una distancia de 13 unidades entre el dato real y el dato pronosticado, existe una diferencia bastante grande el modelo predice que una persona puede fumar entre 0 a 13 cigarrillos. Mientras que para el MAE entre el dato real y el dato pronosticado por el modelo hay una diferencia en promedio de 10. Para el MAPE, el error porcentual promedio del modelo es infinito. Us+Um es mayor a Uc hay evidencias de que el modelo es inconsistente para predecir y posee alta volatilidad, Uc = 0.389, hay baja correlación entre los datos pronosticados y los datos reales. De acuerdo al modelo la máxima correlación esperada es de 0.448 por lo que el pronóstico del modelo puede no estar correcto por lo que no vale la pena para predecir.
bind_rows(Resultados_Performance_2) %>%
stargazer(title = "Medidas de Performance Simulación",
type = "text",
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_2 500 0.039 0.018 0.002 0.027 0.051 0.099
## RMSE_2 500 13.479 0.595 11.734 13.064 13.853 15.361
## MAE_2 500 10.728 0.290 9.834 10.528 10.904 11.626
## MAPE_2 500 Inf.000 Inf Inf Inf Inf
## THEIL_2 500 0.528 0.013 0.491 0.519 0.536 0.564
## Um_2 500 0.002 0.003 0.00000 0.0003 0.003 0.021
## Us_2 500 0.597 0.051 0.476 0.563 0.636 0.751
## Uc_2 500 0.405 0.051 0.253 0.367 0.441 0.529
## ---------------------------------------------------------------Para el valor del R2, en promedio el modelo explica el 3.8% de la varianza de la variable endógena ante cambios en las variables exógenas, como mínimo se explica en 1.8% de la varianza y como máximo se puede esperar un poder explicativo del 9.9% de la varianza de la variable endógena. Para el RMSE en promedio hay una distancia de 13 unidades entre el dato real y el dato pronosticado, existe una diferencia bastante grande el modelo predice que una persona puede fumar entre 0 a 13 cigarrillos. Mientras que para el MAE entre el dato real y el dato pronosticado por el modelo hay una diferencia en promedio de 10. Para el MAPE, el error porcentual promedio del modelo es infinito. Us+Um es mayor a Uc hay evidencias de que el modelo es inconsistente para predecir y posee alta volatilidad, Uc = 0.405, hay baja correlación entre los datos pronosticados y los datos reales. De acuerdo al modelo la máxima correlación esperada es de 0.529 por lo que el pronóstico del modelo puede no estar correcto por lo que no vale la pena para predecir.