estimando el modelo

options(scipen = 999999)
library(lmtest)
library(readxl)
ventas_empresa <- read_excel("C:/Users/manue/Desktop/Practicas_ECMA/datos_Rdata/Datos_Examen2/ventas_empresa.xlsx")
library(stargazer)
library(mlbench) #esta librería tiene el data frame BostonHousing
data(ventas_empresa)
#Modelo estimado medv~. indica "medv" en función del resto de variables del dataframe
modelo_ventas<-lm(formula = V~.,data=ventas_empresa)
coeftest(modelo_ventas)
## 
## t test of coefficients:
## 
##              Estimate Std. Error t value    Pr(>|t|)    
## (Intercept) 107.44351   18.05749  5.9501 0.000008084 ***
## C             0.92257    0.22273  4.1420   0.0005047 ***
## P             0.95018    0.15585  6.0969 0.000005859 ***
## M             1.29779    0.43073  3.0130   0.0068718 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Predicción usando “predict” de “R” base

library(stargazer)
#Data para la predicción X'm
X_m<-data.frame(C=0.05,P=15,M=2)
# Intervalos de Confianza del 95% y del 99%
confidense<-c(0.95,0.99)
#Predicción usando predict
predict(object = modelo_ventas,
           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 = "html")
Pronósticos e intervalos de confianza
Ym Li Ls
95 124.338 81.623 167.053
99 124.338 66.072 182.603

Predicción usando librería “forecast”

library(forecast)
library(kableExtra)
#Data para la predicción X'm
X_m<-data.frame(C=0.05,P=15,M=2)
#Nivel de confianza para el intervalo de confianza
confidense<-c(0.95,0.99)

#Realizando el pronóstico con forecast
pronosticos<-forecast(object = modelo_ventas,
         level = confidense,
         newdata = X_m,ts = FALSE)
kable(pronosticos,
      caption = "Pronóstico e intervalos de confianza:",
      digits = 2,format = "html") #Poner results='asis' en opciones del chunk
Pronóstico e intervalos de confianza:
Point Forecast Lo 95 Hi 95 Lo 99 Hi 99
124.34 81.62 167.05 66.07 182.6

Simulació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)))
}
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<-100 # 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<- ventas_empresa$V %>%
  createDataPartition(p = 0.7,
                      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<- ventas_empresa[muestras[[j]], ]
Datos_Prueba<- ventas_empresa[-muestras[[j]], ]
Modelos_Entrenamiento[[j]]<-lm(formula = V~.,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$V),
            RMSE = RMSE(Modelos_Entrenamiento[[j]]$fitted.values,
                        Datos_Entrenamiento$V),
            MAE = MAE(Modelos_Entrenamiento[[j]]$fitted.values,
                      Datos_Entrenamiento$V),
            MAPE= MAPE(Modelos_Entrenamiento[[j]]$fitted.values,
                       Datos_Entrenamiento$V)*100,
            THEIL=TheilU(Modelos_Entrenamiento[[j]]$fitted.values,
                         Datos_Entrenamiento$V,type = 1),
            Um=Um(Modelos_Entrenamiento[[j]]$fitted.values,
                         Datos_Entrenamiento$V),
            Us=Us(Modelos_Entrenamiento[[j]]$fitted.values,
                         Datos_Entrenamiento$V),
            Uc=Uc(Modelos_Entrenamiento[[j]]$fitted.values,
                         Datos_Entrenamiento$V)
            )
Resultados_Performance[[j]]<-data.frame( 
            R2 = R2(Pronostico_Prueba[[j]], Datos_Prueba$V),
            RMSE = RMSE(Pronostico_Prueba[[j]], Datos_Prueba$V),
            MAE = MAE(Pronostico_Prueba[[j]], Datos_Prueba$V),
            MAPE= MAPE(Pronostico_Prueba[[j]], Datos_Prueba$V)*100,
            THEIL=TheilU(Pronostico_Prueba[[j]], Datos_Prueba$V,
                         type = 1), # También se puede usar la función que creamos: THEIL_U
            Um=Um(Pronostico_Prueba[[j]], Datos_Prueba$V),
            Us=Us(Pronostico_Prueba[[j]], Datos_Prueba$V),
            Uc=Uc(Pronostico_Prueba[[j]], Datos_Prueba$V)
            )
} #No olvidar este corchete ;)

Resultados de la Simulación

bind_rows(Resultados_Performance_data_entrenamiento) %>% 
  stargazer(title = "Medidas de Performance Datos del Modelo",
            type = "html",
            digits = 3)
Medidas de Performance Datos del Modelo
Statistic N Mean St. Dev. Min Pctl(25) Pctl(75) Max
R2 100 0.981 0.003 0.973 0.979 0.983 0.987
RMSE 100 8.480 0.493 7.400 8.164 8.889 9.315
MAE 100 7.273 0.512 6.035 6.947 7.661 8.321
MAPE 100 1.130 0.076 0.949 1.085 1.187 1.285
THEIL 100 0.006 0.0004 0.006 0.006 0.007 0.007
Um 100 0.000 0.000 0 0 0 0
Us 100 0.005 0.001 0.004 0.005 0.006 0.007
Uc 100 1.047 0.001 1.045 1.047 1.048 1.049 
bind_rows(Resultados_Performance) %>% 
  stargazer(title = "Medidas de Performance Simulación",
            type = "html",
            digits = 3)
Medidas de Performance Simulación
Statistic N Mean St. Dev. Min Pctl(25) Pctl(75) Max
R2 100 0.979 0.021 0.883 0.971 0.991 0.999
RMSE 100 10.191 2.570 4.404 8.214 11.777 15.779
MAE 100 9.069 2.553 3.541 7.164 10.821 14.628
MAPE 100 1.398 0.386 0.573 1.127 1.677 2.247
THEIL 100 0.008 0.002 0.003 0.006 0.009 0.012
Um 100 0.258 0.247 0.0002 0.045 0.412 0.931
Us 100 0.258 0.236 0.00000 0.062 0.395 1.036
Uc 100 0.731 0.373 0.058 0.382 1.044 1.333