Pronóstico y Simulación

Verónica Nayeli Miranda Mejía

2024-06-15

UNIVERSIDAD DE EL SALVADOR

FACULTAD DE CIENCIAS ECONÓMICAS

ESCUELA DE ECONOMÍA

ECONOMETRÍA

TEMA:

“Ejercicio de Pronóstico y Simulación”

DOCENTE:

MSF. Carlos Ademir Pérez Alas.

ESTUDIANTE:

Verónica Nayeli Miranda Mejia

Grupo teórico

GT-03

Ciudad Universitaria, Sábado 15 de junio de 2024

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} \]

coeftest(modelo_boston)
## 
## t test of coefficients:
## 
##                Estimate  Std. Error  t value  Pr(>|t|)    
## (Intercept)  3.6459e+01  5.1035e+00   7.1441 3.283e-12 ***
## crim        -1.0801e-01  3.2865e-02  -3.2865 0.0010868 ** 
## zn           4.6420e-02  1.3727e-02   3.3816 0.0007781 ***
## indus        2.0559e-02  6.1496e-02   0.3343 0.7382881    
## chas1        2.6867e+00  8.6158e-01   3.1184 0.0019250 ** 
## nox         -1.7767e+01  3.8197e+00  -4.6513 4.246e-06 ***
## rm           3.8099e+00  4.1793e-01   9.1161 < 2.2e-16 ***
## age          6.9222e-04  1.3210e-02   0.0524 0.9582293    
## dis         -1.4756e+00  1.9945e-01  -7.3980 6.013e-13 ***
## rad          3.0605e-01  6.6346e-02   4.6129 5.071e-06 ***
## tax         -1.2335e-02  3.7605e-03  -3.2800 0.0011116 ** 
## ptratio     -9.5275e-01  1.3083e-01  -7.2825 1.309e-12 ***
## b            9.3117e-03  2.6860e-03   3.4668 0.0005729 ***
## lstat       -5.2476e-01  5.0715e-02 -10.3471 < 2.2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Pronóstico en R

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

library(stargazer)
#Data para la predicción X'm
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)
# Intervalos de Confianza del 95% y del 99%
confidense<-c(0.95,0.99)
#Predicción usando predict
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 = "Pronosticos e Intervalos de Confianza",
          type = "html") #Poner results='asis' en opciones del chunk
Pronosticos e Intervalos de Confianza
Ym Li Ls
95 26.116 15.558 36.673
99 26.116 12.221 40.010

Interpretación

Al utilizar este modelo, en un intervalo de confianza del 95% la valoración de la casa sería de $261,000.00 dólares, asimismo puede tener un valor mínimo de $155,000.00 y un valor máximo de $366,000.00. Mientras que para un intervalo de confianza del 99%, la estimación puntual no cambia, pero si cambia el intervalo de confianza siendo más ancho, donde el valor minimo sería de $122,000.00 y el valor máximo de $40,000.00.

Predicción usando librería “forecast”

library(forecast)
library(kableExtra)
#Data para la predicción X'm
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)
#Nivel de confianza para el intervalo de confianza
confidense<-c(0.95,0.99)

#Realizando el pronóstico con forecast
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") #Poner results='asis' en opciones del chunk
Pronóstico e intervalos de confianza:
Point Forecast Lo 95 Hi 95 Lo 99 Hi 99
26.12 15.56 36.67 12.22 40.01

Interpretación

Al utilizar este modelo, en un intervalo de confianza del 95% la valoración de la casa sería de $261,000.00 dólares, asimismo puede tener un valor mínimo de $155,000.00 y un valor máximo de $366,000.00. Mientras que para un intervalo de confianza del 99%, la estimación puntual no cambia, pero si cambia el intervalo de confianza siendo más ancho, donde el valor minimo sería de $122,000.00 y el valor máximo de $40,000.00.

Ejemplo Simulación en R

Definición de funciones para el cálculo Theil y su descomposición, debe estar instalada la librería DescTools

#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) #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 ;)

Resultados de la Simulación

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 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

Interpretación

Al observar Theil, se muestra que tiene un valor finito menor que uno, es decir que tiene un buen comportamiento.

La proporción de sesgo en datos del modelo estimado presenta ceros.

El 74.3% de la varianza en Y, es explicada por los regresores.

La explicación de la varianza podria estar entre 71.3% y 79.4%.

Según MAE, el modelo se equivoca en $3,200.00 dólares al momento de la evaluación.

Según MAPE el error mínimo esperado es de 14.81 para la data de estrenamiento y el máximo es de 17.70

El pronóstico de varianza esta entre un mínimo de 0.05 y un máximo de 0.08.

El pronóstico de covarianza esta alto, se encuentra entre un mínimo de 0.918 y un máximo de 0.945.

La suma de proporción de sesgo y de varianza da menos de 0.20, lo que significa que es un modelo bastante bueno.

bind_rows(Resultados_Performance) %>%
  stargazer(title = "Medidas de Performance Simulacion",
            type = "html",
            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

Interpretación

Según MAE, el modelo se equivoca en $3,400.00 dólares al momento de la evaluación.

Según MAPE el error mínimo esperado es de 12.875 para la data de estrenamiento y el máximo es de 23.137

La diferencia de correlación es de 0.01

La distancia promedio entre las medias de los datos observados y proyectados del modelo es de 0.011.

Las distancias en las volatilidades es pequeña, en general la distancia entre sus medias y las volatilidades, es inferior al 0.20 por lo tanto se garantiza que los pronósticos que se hicieron al inicio son bastantes certeros.

Con los datos del ejercicio de estimadores HAC:

Realice una simulación del pronóstico, usando como datos de entrenamiento el 75% de los datos originales, considere 500 replicas y analice los resultados del desempeño del modelo con los datos “desconocidos”. En ambos casos utilice la semilla aleatoria de 50

Estimación del Modelo

options(scipen = 999999)
library(stargazer)
load("C:/Users/User/Desktop/documents_rstudio/periodo_3/Verónica Nayeli Miranda Mejia - smoke.RData")
ata<-data
equation<-as.formula("cigs~cigpric+lcigpric+income+lincome+age+agesq+educ+white+restaurn") 
endogena<-data$cigs 
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.01

Simulación en R

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 = "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_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

Interpretación

Al observar Theil, se muestra que tiene un valor finito menor que uno, es decir que tiene un buen comportamiento.

La proporción de sesgo en datos del modelo estimado presenta ceros.

El 0.05 de la varianza en Y, es explicada por los regresores.

La explicación de la varianza podria estar entre 4% y 8%.

Según MAE, el modelo se equivoca en $10.56 dólares al momento de la evaluación.

Según MAPE el error mínimo esperado es infinito para la data de estrenamiento y el máximo es también infinito.

El pronóstico de varianza esta entre un mínimo de 0.55 y un máximo de 0.66.

El pronóstico de covarianza esta bajo, se encuentra entre un mínimo de 0.332 y un máximo de 0.44.

La suma de proporción de sesgo y de varianza no da menos de 0.20, lo que significa que no es un modelo bastante bueno.

bind_rows(Resultados_Performance_2) %>% 
  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_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

Interpretación

Según MAE, el modelo se equivoca en $10.73 dólares al momento de la evaluación.

Según MAPE el error mínimo y máximo esta en inf.

La diferencia de correlación es de 0.016

La distancia promedio entre las medias de los datos observados y proyectados del modelo es de 0.002.

Las distancias en las volatilidades no es pequeña, en general la distancia entre sus medias y las volatilidades, es superior al 0.20 por lo tanto se garantiza que los pronósticos que se hicieron al inicio no son bastantes certeros.