G02_MR14076_P3_C
LISBETH ESMERALDA MELCHOR REYES
2022-07-15
EJERCICIO 1
1. Estimación del modelo
library(stargazer)
load ("C:/Users/liizm/Downloads/datos_parcial.rData")
Dataframe <- as.data.frame (sat)
Modelo_total <- lm(formula = total ~ expend + ratio + salary + takers,
data = Dataframe)2. Capacidad predictiva del modelo usando muestra de entrenamiento al 80% total de la data, 150 iteraciones y seed de 35
options(scipen = 999999)
library(dplyr)
library(caret)
library(DescTools)
library(stargazer)
set.seed(35)
numero_de_muestras<-150
muestras<- Dataframe$total %>%
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)2
#Simulación y Medidas de Desempeño.
#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)))
}# Códigos para hacer la Simulación:
options(scipen = 999999)
library(dplyr)
library(caret)
library(DescTools)
library(stargazer)
set.seed(35)
numero_de_muestras<-150
muestras<- Dataframe$total %>%
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<- Dataframe[muestras[[j]], ]
Datos_Prueba<- Dataframe[-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$total),
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$total,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
bind_rows(Resultados_Performance_data_entrenamiento) %>%
stargazer(title = "Medidas de Performance Datos del Modelo Total",
type = "text",
digits = 3)##
## Medidas de Performance Datos del Modelo Total
## ===========================================
## Statistic N Mean St. Dev. Min Max
## -------------------------------------------
## R2 150 0.831 0.019 0.800 0.881
## RMSE 150 30.449 1.893 25.018 33.464
## MAE 150 23.638 1.548 19.172 27.036
## MAPE 150 2.455 0.158 2.007 2.807
## THEIL 150 0.016 0.001 0.013 0.017
## 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
## -------------------------------------------bind_rows(Resultados_Performance) %>%
stargazer(title = "Medidas de Performance Simulación del Modelo Total",
type = "text",
digits = 3)##
## Medidas de Performance Simulación del Modelo Total
## ============================================
## 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
## --------------------------------------------EJERCICIO 2
library(dplyr)
library(tidyr)
library(forecast)
library(TSstudio)
AirPassengers.ts<- ts(data = AirPassengers,
start = c(1949, 1),
frequency = 12)
AirPassengers.ts %>%
autoplot(main= "Datos de vuelos de pasajeros de aerolineas 1949-1960",
xlab = "Años/Meses",
ylab = "Vuelos")
library(forecast)
library(TSstudio)
Yt <-AirPassengers.ts
ts_plot(Yt,Xtitle = "Años/Meses")Estimación del modelo HW
library(forecast)
ModeloHW<-HoltWinters(x = Yt,
seasonal = "multiplicative",
optim.start = c(0.9,0.9,0.9))
ModeloHW## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
##
## Call:
## HoltWinters(x = Yt, seasonal = "multiplicative", optim.start = c(0.9, 0.9, 0.9))
##
## Smoothing parameters:
## alpha: 0.2755935
## beta : 0.03269302
## gamma: 0.8707316
##
## Coefficients:
## [,1]
## a 469.3232968
## b 3.0215403
## s1 0.9464609
## s2 0.8829237
## s3 0.9717366
## s4 1.0304824
## s5 1.0476882
## s6 1.1805270
## s7 1.3590775
## s8 1.3331703
## s9 1.1083379
## s10 0.9868812
## s11 0.8361331
## s12 0.9209875PronosticosHW<-forecast(object = ModeloHW,h = 12,level = c(0.95))
PronosticosHW## Point Forecast Lo 95 Hi 95
## Jan 1961 447.0559 427.3061 466.8057
## Feb 1961 419.7123 399.1325 440.2920
## Mar 1961 464.8671 442.9629 486.7712
## Apr 1961 496.0840 472.8329 519.3351
## May 1961 507.5326 483.1374 531.9278
## Jun 1961 575.4509 548.7081 602.1936
## Jul 1961 666.5923 636.6287 696.5558
## Aug 1961 657.9137 627.1820 688.6454
## Sep 1961 550.3088 521.6397 578.9778
## Oct 1961 492.9853 465.0152 520.9554
## Nov 1961 420.2073 393.4687 446.9459
## Dec 1961 465.6345 443.3003 487.9687Gráfico de la serie original y del pronóstico.
PronosticosHW %>% autoplot()
Orden de integración de la serie
library(kableExtra)
library(magrittr)
d<-ndiffs(Yt)
D<-nsdiffs(Yt)
ordenes_integracion<-c(d,D)
names(ordenes_integracion)<-c("d","D")
ordenes_integracion %>% kable(caption = "Ordenes de Integración") %>% kable_material()| x | |
|---|---|
| d | 1 |
| D | 1 |
Gráfico de la serie diferenciada
Yt %>%
diff(lag = 12,diffences=D) %>%
diff(diffences=d) %>%
ts_plot(title = "Yt estacionaria")Valores para (p,q) & (P,Q)
Yt %>%
diff(lag = 12,diffences=D) %>%
diff(diffences=d) %>%
ts_cor(lag.max = 72)Estimación del modelo Sarima usando la líbreria forecast
library(forecast)
library(ggthemes)
modelo_estimado <- Yt %>%
Arima(order = c(0, 1, 0),
seasonal = c(1, 1, 1))
summary(modelo_estimado)## Series: .
## ARIMA(0,1,0)(1,1,1)[12]
##
## Coefficients:
## sar1 sma1
## -0.4451 0.2752
## s.e. 0.4412 0.4643
##
## sigma^2 = 149.5: log likelihood = -513.06
## AIC=1032.13 AICc=1032.32 BIC=1040.75
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.2422851 11.57129 8.298806 0.03881706 2.936113 0.2590923
## ACF1
## Training set -0.3067414modelo_estimado %>% autoplot(type="both")+theme_solarized()
modelo_estimado %>% check_res(lag.max = 72)Yt_Sarima<-modelo_estimado$fitted
Yt_HW<-PronosticosHW$fitted
grafico_comparativo<-cbind(Yt,Yt_Sarima,Yt_HW)
ts_plot(grafico_comparativo)Verificación de sobre ajuste/sub ajuste
library(tsibble)
library(feasts)
library(fable)
library(fabletools)
library(tidyr)
library(dplyr)
a<-Yt %>% as_tsibble() %>%
model(arima_original=ARIMA(value ~ pdq(0, 1, 0) + PDQ(1, 1, 1)),
arima_010_011 = ARIMA(value ~ pdq(0, 1, 0) + PDQ(0, 1, 1)),
arima_010_110 = ARIMA(value ~ pdq(0, 1, 0) + PDQ(1, 1, 0)),
arima_automatico=ARIMA(value,ic="bic",stepwise = FALSE)
)
print(a)## # A mable: 1 x 4
## arima_original arima_010_011 arima_010_110
## <model> <model> <model>
## 1 <ARIMA(0,1,0)(1,1,1)[12]> <ARIMA(0,1,0)(0,1,1)[12]> <ARIMA(0,1,0)(1,1,0)[12]>
## # ... with 1 more variable: arima_automatico <model>a %>% pivot_longer(everything(), names_to = "Model name",
values_to = "Orders") %>% glance() %>%
arrange(AICc) ->tabla
tabla## # A tibble: 4 x 9
## `Model name` .model sigma2 log_lik AIC AICc BIC ar_roots ma_roots
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <list> <list>
## 1 arima_automatico Orders 138. -508. 1020. 1020. 1026. <cpl [1]> <cpl [0]>
## 2 arima_010_110 Orders 149. -513. 1030. 1031. 1036. <cpl [12]> <cpl [0]>
## 3 arima_010_011 Orders 149. -513. 1031. 1031. 1037. <cpl [0]> <cpl [12]>
## 4 arima_original Orders 149. -513. 1032. 1032. 1041. <cpl [12]> <cpl [12]>library(forecast)
library(dplyr)
library(tsibble)
library(fable)
library(fabletools)
Yt<- Yt %>% as_tsibble() %>% rename(AirPassengers= value)
data.cross.validation<-Yt %>%
as_tsibble() %>%
stretch_tsibble(.init = 24,.step = 1)
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.194 13.1 9.92 -0.0727 3.12 0.310 0.362 -0.267