Con los datos utilizados en la presentación de pronóstico y simulación:
options(scipen = 999999)
library(lmtest)
library(stargazer)
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)
stargazer(modelo_boston,title = "Modelo estimado", type = "text")##
## Modelo estimado
## ===============================================
## Dependent variable:
## ---------------------------
## medv
## -----------------------------------------------
## crim -0.108***
## (0.033)
##
## zn 0.046***
## (0.014)
##
## indus 0.021
## (0.061)
##
## chas1 2.687***
## (0.862)
##
## nox -17.767***
## (3.820)
##
## rm 3.810***
## (0.418)
##
## age 0.001
## (0.013)
##
## dis -1.476***
## (0.199)
##
## rad 0.306***
## (0.066)
##
## tax -0.012***
## (0.004)
##
## ptratio -0.953***
## (0.131)
##
## b 0.009***
## (0.003)
##
## lstat -0.525***
## (0.051)
##
## Constant 36.459***
## (5.103)
##
## -----------------------------------------------
## Observations 506
## R2 0.741
## Adjusted R2 0.734
## Residual Std. Error 4.745 (df = 492)
## F Statistic 108.077*** (df = 13; 492)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01coeftest(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 ' ' 1library(stargazer)
X_m<-data.frame(crim=0.05,zn=15,indus=2,chas="0",nox=0.004,
rm=5,age=85,dis=5.56,rad=2,tax=300,ptratio=17,b=0.00005,lstat=5)
confidense<-c(0.95,0.99)
predict(object = modelo_boston,
newdata = X_m,
interval = "prediction",
level = confidense,
se.fit =TRUE)->predicciones
rownames(predicciones$fit)<-as.character(confidense*100)
colnames(predicciones$fit)<-c("Ym","Li","Ls")
stargazer(predicciones$fit,
title = "Pronósticos e intervalos de confianza",
type = "text") 95 26.116 15.558 36.673 99 26.116 12.221 40.010 ———————–
library(forecast)
library(kableExtra)
X_m<-data.frame(crim=0.05,zn=15,indus=2,chas="0",nox=0.004,
rm=5,age=85,dis=5.56,rad=2,tax=300,ptratio=17,b=0.00005,lstat=5)
confidense<-c(0.95,0.99)
pronosticos<-forecast(object = modelo_boston,
level = confidense,
newdata = X_m,ts = FALSE)
kable(pronosticos,
caption = "Pronóstico e intervalos de confianza:",
digits = 2,format = "html") | Point Forecast | Lo 95 | Hi 95 | Lo 99 | Hi 99 |
|---|---|---|---|---|
| 26.12 | 15.56 | 36.67 | 12.22 | 40.01 |
#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)))
}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) # 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 ;)bind_rows(Resultados_Performance_data_entrenamiento) %>%
stargazer(title = "Medidas de Performance Datos del Modelo",
type = "text",
digits = 3)##
## Medidas de Performance Datos del Modelo
## =============================================
## Statistic N Mean St. Dev. Min Max
## ---------------------------------------------
## R2 1,000 0.743 0.013 0.713 0.794
## RMSE 1,000 4.653 0.141 4.177 4.948
## MAE 1,000 3.265 0.095 2.905 3.512
## MAPE 1,000 16.387 0.464 14.813 17.691
## THEIL 1,000 0.096 0.003 0.087 0.102
## Um 1,000 0.000 0.000 0 0
## Us 1,000 0.074 0.004 0.058 0.085
## Uc 1,000 0.928 0.004 0.918 0.945
## ---------------------------------------------Interpretacion de performance.
En promedio el 74.3% de la varianza de la endogena, es explicada por los regresores del modelo, con un valor porcentual minimo explicado de 71.3% y un valor maximo de 79.4%.
El error porcentual promedio del modelo es de 1638.7%, el error minimo es de 1481.3% y el maximo es de 1769.1%.
La distancia promedio entre la varianza de los datos reales y los pronosticados es de 2.9.
El valor medio de corrrelacion entre el valor real y el pronosticado es de 92.8%, con un valor minimo de 91.8% y un maximo de 94.5%
bind_rows(Resultados_Performance) %>%
stargazer(title = "Medidas de Performance Simulación",
type = "text",
digits = 3)##
## Medidas de Performance Simulación
## ==============================================
## Statistic N Mean St. Dev. Min Max
## ----------------------------------------------
## R2 1,000 0.723 0.056 0.452 0.840
## RMSE 1,000 4.862 0.575 3.465 6.961
## MAE 1,000 3.411 0.281 2.633 4.492
## MAPE 1,000 17.197 1.618 12.875 23.137
## THEIL 1,000 0.101 0.012 0.073 0.148
## Um 1,000 0.011 0.016 0.000 0.205
## Us 1,000 0.081 0.066 0.00000 0.333
## Uc 1,000 0.918 0.066 0.667 1.010
## ----------------------------------------------En ambos casos utilice la semilla aleatoria de 50
load("C:/Users/MINEDUCYT/Desktop/ECONOMETRIA/Archivos R/-Template- smoke.RData")
options(scipen = 99999999)
library(lmtest)
library(dplyr)
modelo_Simulacion<-lm(cigs~cigpric+lcigpric+income+lincome+age+agesq+educ+white+restaurn, data = data)
coeftest(modelo_Simulacion)##
## 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 ' ' 1library(stargazer)
stargazer(modelo_Simulacion,title = "Modelo estimado", type = "text")##
## Modelo estimado
## ===============================================
## Dependent variable:
## ---------------------------
## cigs
## -----------------------------------------------
## cigpric 2.002
## (1.493)
##
## lcigpric -115.273
## (85.424)
##
## income -0.00005
## (0.0001)
##
## lincome 1.404
## (1.708)
##
## age 0.778***
## (0.161)
##
## agesq -0.009***
## (0.002)
##
## educ -0.495***
## (0.168)
##
## white -0.531
## (1.461)
##
## restaurn -2.644**
## (1.130)
##
## Constant 340.804
## (260.016)
##
## -----------------------------------------------
## Observations 807
## R2 0.055
## Adjusted R2 0.044
## Residual Std. Error 13.413 (df = 797)
## F Statistic 5.169*** (df = 9; 797)
## ===============================================
## Note: *p<0.1; **p<0.05; ***p<0.01library(stargazer)
X_m<-data.frame(cigpric=0.05,lcigpric=15,income=2,lincome=0,age=0.004,
agesq=5,educ=85,white=5.56,restaurn=2)
confidense<-c(0.95,0.99)
predict(object = modelo_Simulacion,
newdata = X_m,
interval = "prediction",
level = confidense,
se.fit =TRUE)->predicciones
rownames(predicciones$fit)<-as.character(confidense*100)
colnames(predicciones$fit)<-c("Ym","Li","Ls")
stargazer(predicciones$fit,
title = "Pronósticos e intervalos de confianza",
type = "text") 95 -1,438.538 -3,443.844 566.769 99 -1,438.538 -4,076.271 1,199.196 ———————————-
library(forecast)
library(kableExtra)
X_m<-data.frame(cigpric=0.05,lcigpric=15,income=2,lincome=0,age=0.004,
agesq=5,educ=85,white=5.56,restaurn=2)
confidense<-c(0.95,0.99)
pronosticos<-forecast(object = modelo_Simulacion,
level = confidense,
newdata = X_m,ts = FALSE)
kable(pronosticos,
caption = "Pronóstico e intervalos de confianza:",
digits = 2,format = "html") | Point Forecast | Lo 95 | Hi 95 | Lo 99 | Hi 99 |
|---|---|---|---|---|
| -1438.54 | -3443.84 | 566.77 | -4076.27 | 1199.2 |
#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)))
}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) # 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<-500 # Numero de muestras que se optendran del data frame
# 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_Entrenamiento2<-vector(mode = "list",
length = numero_de_muestras)
Pronostico_Prueba2<-vector(mode = "list",
length = numero_de_muestras)
Resultados_Performance_data_entrenamiento2<-vector(mode = "list",
length = numero_de_muestras)
Resultados_Performance2<-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_simulacion<- data[muestras[[j]], ]
Datos_Prueba2<- data[-muestras[[j]], ]
Modelos_Entrenamiento2[[j]]<-lm(formula = cigs~.,data=Datos_simulacion)
Pronostico_Prueba2[[j]]<-Modelos_Entrenamiento2[[j]] %>% predict(Datos_Prueba2)
Resultados_Performance_data_entrenamiento2[[j]]<-data.frame(
R2 = R2(Modelos_Entrenamiento2[[j]]$fitted.values,
Datos_simulacion$cigs),
RMSE = RMSE(Modelos_Entrenamiento2[[j]]$fitted.values,
Datos_simulacion$cigs),
MAE = MAE(Modelos_Entrenamiento2[[j]]$fitted.values,
Datos_simulacion$cigs),
MAPE= MAPE(Modelos_Entrenamiento2[[j]]$fitted.values,
Datos_simulacion$cigs)*100,
THEIL=TheilU(Modelos_Entrenamiento2[[j]]$fitted.values,
Datos_simulacion$cigs,type = 1),
Um=Um(Modelos_Entrenamiento2[[j]]$fitted.values,
Datos_simulacion$cigs),
Us=Us(Modelos_Entrenamiento2[[j]]$fitted.values,
Datos_simulacion$cigs),
Uc=Uc(Modelos_Entrenamiento2[[j]]$fitted.values,
Datos_simulacion$cigs)
)
Resultados_Performance2[[j]]<-data.frame(
R2 = R2(Pronostico_Prueba2[[j]], Datos_Prueba2$cigs),
RMSE = RMSE(Pronostico_Prueba2[[j]], Datos_Prueba2$cigs),
MAE = MAE(Pronostico_Prueba2[[j]], Datos_Prueba2$cigs),
MAPE= MAPE(Pronostico_Prueba2[[j]], Datos_Prueba2$cigs)*100,
THEIL=TheilU(Pronostico_Prueba2[[j]], Datos_Prueba2$cigs,
type = 1), # También se puede usar la función que creamos: THEIL_U
Um=Um(Pronostico_Prueba2[[j]], Datos_Prueba2$cigs),
Us=Us(Pronostico_Prueba2[[j]], Datos_Prueba2$cigs),
Uc=Uc(Pronostico_Prueba2[[j]], Datos_Prueba2$cigs)
)
} #No olvidar este corchete ;)bind_rows(Resultados_Performance_data_entrenamiento2) %>%
stargazer(title = "Medidas de Performance Datos del Modelo",
type = "text",
digits = 3,summary.stat = c("n","mean","sd","p25","p75","min","max"))##
## Medidas de Performance Datos del Modelo
## ==============================================================
## Statistic N Mean St. Dev. Pctl(25) Pctl(75) Min Max
## --------------------------------------------------------------
## R2 500 0.058 0.007 0.054 0.063 0.040 0.083
## RMSE 500 13.304 0.198 13.175 13.439 12.629 13.828
## MAE 500 10.562 0.133 10.474 10.654 10.068 10.969
## MAPE 500 Inf.000 Inf Inf Inf Inf
## THEIL 500 0.522 0.006 0.517 0.526 0.505 0.541
## Um 500 0.000 0.000 0 0 0 0
## Us 500 0.613 0.020 0.601 0.625 0.554 0.669
## Uc 500 0.389 0.020 0.376 0.401 0.332 0.448
## --------------------------------------------------------------En promedio el 50.3% de la varianza de la endogena es explicada por los regresores del modelo, con un valor porcentual minimo de 4% y un maximo de 8.3%.
La distancia promedio entre la varianza de los datos reales y los pronosticos es de 10.06.
El valor medio de correlacion entre el valor real y el pronosticado es de 38.9%, con un valor minimo de 33.2% y un valor maximo de 44.8%.
bind_rows(Resultados_Performance2) %>%
stargazer(title = "Medidas de Performance Simulación",
type = "text",
digits = 3,summary.stat = c("n","mean","sd","p25","p75","min","max"))##
## Medidas de Performance Simulación
## ===============================================================
## Statistic N Mean St. Dev. Pctl(25) Pctl(75) Min Max
## ---------------------------------------------------------------
## R2 500 0.039 0.018 0.027 0.051 0.002 0.099
## RMSE 500 13.479 0.595 13.064 13.853 11.734 15.361
## MAE 500 10.728 0.290 10.528 10.904 9.834 11.626
## MAPE 500 Inf.000 Inf Inf Inf Inf
## THEIL 500 0.528 0.013 0.519 0.536 0.491 0.564
## Um 500 0.002 0.003 0.0003 0.003 0.00000 0.021
## Us 500 0.597 0.051 0.563 0.636 0.476 0.751
## Uc 500 0.405 0.051 0.367 0.441 0.253 0.529
## ---------------------------------------------------------------En promedio el 3.9% de la varianza de la endogena, es explicada por los regresores del modelo, con un valor porcentual minimo explicado de 0.2% y un valor maximo de 9.9%.
La distancia promedio entre la varianza de los datos reales y los pronosticados es de 9.83.
El valor medio de corrrelacion entre el valor real y el pronosticado es de 40.5%, con un valor minimo de 25.3% y un maximo de 52.9%