#EJERCICIO 1

##1)ESTIMANDO MODELO PARA EXPLICAR LOS SALARIOS

#importacion de datos
options(scipen = 999999)
load("C:/doc R/datos_parcial.RData")
library(lmtest)
modelo.wage<- lm(total~expend+ratio+salary+takers, data = sat)
summary(modelo.wage)
## 
## Call:
## lm(formula = total ~ expend + ratio + salary + takers, data = sat)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -90.531 -20.855  -1.746  15.979  66.571 
## 
## Coefficients:
##              Estimate Std. Error t value             Pr(>|t|)    
## (Intercept) 1045.9715    52.8698  19.784 < 0.0000000000000002 ***
## expend         4.4626    10.5465   0.423                0.674    
## ratio         -3.6242     3.2154  -1.127                0.266    
## salary         1.6379     2.3872   0.686                0.496    
## takers        -2.9045     0.2313 -12.559 0.000000000000000261 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 32.7 on 45 degrees of freedom
## Multiple R-squared:  0.8246, Adjusted R-squared:  0.809 
## F-statistic: 52.88 on 4 and 45 DF,  p-value: < 0.00000000000000022

##2)VERIFICANDO LA CAPACIDAD PREDICTIVA DEL MODELO

###definicion de las funciones para el calculo

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

###ejecutando simulacion

library(dplyr)
library(magrittr)
library(caret) #Realizar muestreo
library(DescTools) # Funciones para calcular Performance
library(stargazer) # Formato, y Resumen estadistico de las simulaciones

set.seed(35)
numero_de_muestras<-150


muestras<- sat$total %>%
  createDataPartition(p = 0.8,                        
                      times = numero_de_muestras,
                      list = TRUE) 


# Listas para datos de entrenamiento
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<- sat[muestras[[j]], ]
Datos_Prueba<- sat[-muestras[[j]], ]
Modelos_Entrenamiento[[j]]<-lm(formula =total~expend+ratio+salary+takers, 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$total),
            RMSE = RMSE(Modelos_Entrenamiento[[j]]$fitted.values,
                        Datos_Entrenamiento$wage),
            MAE = MAE(Modelos_Entrenamiento[[j]]$fitted.values,
                      Datos_Entrenamiento$total),
            MAPE= MAPE(Modelos_Entrenamiento[[j]]$fitted.values,
                       Datos_Entrenamiento$total)*100,
            THEIL=TheilU(Modelos_Entrenamiento[[j]]$fitted.values,
                         Datos_Entrenamiento$wage,type = 1),
            Um=Um(Modelos_Entrenamiento[[j]]$fitted.values,
                         Datos_Entrenamiento$total),
            Us=Us(Modelos_Entrenamiento[[j]]$fitted.values,
                         Datos_Entrenamiento$total),
            Uc=Uc(Modelos_Entrenamiento[[j]]$fitted.values,
                         Datos_Entrenamiento$total)
            )
Resultados_Performance[[j]]<-data.frame( 
            R2 = R2(Pronostico_Prueba[[j]], Datos_Prueba$total),
            RMSE = RMSE(Pronostico_Prueba[[j]], Datos_Prueba$total),
            MAE = MAE(Pronostico_Prueba[[j]], Datos_Prueba$total),
            MAPE= MAPE(Pronostico_Prueba[[j]], Datos_Prueba$total)*100,
            THEIL=TheilU(Pronostico_Prueba[[j]], Datos_Prueba$total,
                         type = 1), 
            Um=Um(Pronostico_Prueba[[j]], Datos_Prueba$total),
            Us=Us(Pronostico_Prueba[[j]], Datos_Prueba$total),
            Uc=Uc(Pronostico_Prueba[[j]], Datos_Prueba$total)
            )
} 
# Resultados performance Data entrenamiento
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 Max
R2 150 0.831 0.019 0.800 0.881
MAE 150 23.638 1.548 19.172 27.036
MAPE 150 2.455 0.158 2.007 2.807
Um 150 0.000 0.000 0 0
Us 150 0.047 0.006 0.032 0.057
Uc 150 0.977 0.006 0.967 0.992 
# Resultados perfmormance
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 Max
R2 150 0.806 0.110 0.458 0.974
RMSE 150 34.465 9.848 12.123 56.590
MAE 150 27.305 7.408 9.242 45.964
MAPE 150 2.837 0.753 0.961 4.676
THEIL 150 0.018 0.005 0.006 0.029
Um 150 0.116 0.139 0.00000 0.681
Us 150 0.129 0.155 0.00000 0.810
Uc 150 0.881 0.213 0.310 1.139

##3)error porcentual maximo esperado de preiccipn MAPE = 2.807

##4)error absoluto minimo esperado de prediccion MAE = 19.172

##5)El porcentaje de la varianza de la variable endógena que es explicada por los pronósticos del modelo está entre 0.458 (valor mínimo) y 0.974 (valor máximo) y la media del R cuadraro es 0.806, entonces, en promedio la varianza de la variable endógena es explicada en 80.6%.

##6)El coheficiente de sesgo (Um) en promedio es de 0.116 por tanto el valor pronosticado y el real es bastante bajo, el coheficiente de varianza (Us) en promedio es 0.129 en este caso el modelo explica muy bien la varianza entre los datos reales y estimados y el coheficiente de covarianza (Uc) en promedio es de 0.881 que representa un correlación alta entre el valor pronosticado por el modelo y el real EL MODELO SI ES RECOMENDABLE.

#EJERCICIO 2 ESIMAR MODELO SARIMA

library(forecast)
library(wooldridge)
data("AirPassengers")
AirPassengers<-na.omit(AirPassengers)
Yt<-AirPassengers
AirPassengers%>%autoplot(main = "grafico AirPassengers",xlab = "años/meses",ylab = "pasajeros")

III. Validación Cruzada (Cross Validated)

library(forecast)
library(dplyr)
library(tsibble)
library(fable)
library(fabletools)

#convertimos
Yt<-Yt %>% 
  as_tsibble() %>% 
  rename(AirPassengers=value)

#generador de la dta (roja)
data.cross.validation<-Yt %>% 
  as_tsibble() %>% 
  stretch_tsibble(.init = 60,.step = 1)

#generador de pronostico (verde) y simulacion 
TSCV<-data.cross.validation %>% 
model(ARIMA(AirPassengers ~ pdq(0, 1, 0) + PDQ(1, 1, 1))) %>% 
forecast(h=1) %>% accuracy(Yt)

print(TSCV) 
## # A tibble: 1 x 10
##   .model                 .type    ME  RMSE   MAE    MPE  MAPE  MASE RMSSE   ACF1
##   <chr>                  <chr> <dbl> <dbl> <dbl>  <dbl> <dbl> <dbl> <dbl>  <dbl>
## 1 ARIMA(AirPassengers ~~ Test  0.451  13.8  10.4 0.0390  2.93 0.325 0.381 -0.285