Pronostico serie temporal IVAE 03-2022
ABEL ANTONIO BOLAÑOS PEÑA - BP18001 G03
2022-06-22
IMPORTANDO EL MODELO
library(readxl)
library(tidyr)
library(dplyr)
library(forecast)
IVAE_2022 <- read_excel("C:/doc R/IVAE_2022.xlsx",col_names =FALSE, skip = 6, n_max = 10)
data.ivae<-pivot_longer(data = IVAE_2022[1,],names_to = "vars",cols = 2:160,values_to = "indice") %>% select("indice")
serie.ivae.ts<-data.ivae %>% ts(start = c(2009,1),frequency = 12)
Yt <- serie.ivae.ts
serie.ivae.ts %>%
autoplot(main="IVAE ENERO 2009- MARZ0 2022",
xlab="Años/meses",
ylab="Indice")
I. PRONOSTICO DE SERIES TEMPORALES, ENFOQUE DETERMINISTICO (CLASICO)*
Pronostico modelo de Holt Winters
Usando Stats y Forecast
MULTIPLICATIVO
library(forecast)
#ESTIMAR EL MODELO
modelohw<-HoltWinters(x = Yt,
seasonal= "multiplicative",
optim.start=c(0.9,0.9,0.9)) #OPCIONAL:bUSCA establecer los valores quinceales para la busqueda de parametros que busca HW
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.8472467
## beta : 0
## gamma: 1
##
## Coefficients:
## [,1]
## a 118.3719784
## b 0.1600947
## s1 0.9529095
## s2 1.0192349
## s3 1.0440936
## s4 0.9962326
## s5 1.0009929
## s6 0.9838373
## s7 0.9554392
## s8 1.0145286
## s9 1.0872315
## s10 0.9748891
## s11 0.9626701
## s12 1.0059813#GENERAR EL PRONOSTICO
pronosticoshw<-forecast(object = modelohw,h = 12,level = c(0.95))
pronosticoshw## Point Forecast Lo 95 Hi 95
## Apr 2022 112.9503 107.46492 118.4358
## May 2022 120.9752 113.57246 128.3779
## Jun 2022 124.0929 115.22236 132.9634
## Jul 2022 118.5640 108.86876 128.2592
## Aug 2022 119.2908 108.50116 130.0804
## Sep 2022 117.4038 105.81297 128.9947
## Oct 2022 114.1680 101.97016 126.3657
## Nov 2022 121.3911 107.66977 135.1125
## Dec 2022 130.2643 114.88360 145.6450
## Jan 2023 116.9603 102.34981 131.5708
## Feb 2023 115.6485 100.48405 130.8129
## Mar 2023 121.0126 83.19827 158.8270#GRAFICO D ELA SERIE ORIGINAL Y DEL PRONOSTICO
pronosticoshw %>% autoplot(main="Pronostico HW multiplicativo", xlab="Años/meses")
ADITIVO
library(forecast)
#ESTIMAR EL MODELO
modelohw.ad<-HoltWinters(x = Yt,
seasonal= "additive",
optim.start=c(0.9,0.9,0.9)) #OPCIONAL:bUSCA establecer los valores quinceales para la busqueda de parametros que busca HW
modelohw.ad## Holt-Winters exponential smoothing with trend and additive seasonal component.
##
## Call:
## HoltWinters(x = Yt, seasonal = "additive", optim.start = c(0.9, 0.9, 0.9))
##
## Smoothing parameters:
## alpha: 0.9142733
## beta : 0
## gamma: 1
##
## Coefficients:
## [,1]
## a 118.0938359
## b 0.1600947
## s1 -3.9032233
## s2 3.0467916
## s3 3.9273013
## s4 -1.6040517
## s5 -0.5787492
## s6 -1.8288955
## s7 -4.0808499
## s8 2.2490593
## s9 9.6275698
## s10 -2.7962367
## s11 -4.4557494
## s12 0.9861641#GENERAR EL PRONOSTICO
pronosticoshw.ad<-forecast(object = modelohw.ad,h = 12,level = c(0.95))
pronosticoshw.ad## Point Forecast Lo 95 Hi 95
## Apr 2022 114.3507 109.10096 119.6005
## May 2022 121.4608 114.34766 128.5740
## Jun 2022 122.5014 113.92039 131.0825
## Jul 2022 117.1302 107.29801 126.9623
## Aug 2022 118.3156 107.37443 129.2567
## Sep 2022 117.2255 105.27790 129.1731
## Oct 2022 115.1336 102.25799 128.0093
## Nov 2022 121.6237 107.88248 135.3648
## Dec 2022 129.1623 114.60695 143.7176
## Jan 2023 116.8985 101.57229 132.2248
## Feb 2023 115.3991 99.33889 131.4594
## Mar 2023 121.0011 104.23903 137.7632#GRAFICO D ELA SERIE ORIGINAL Y DEL PRONOSTICO
pronosticoshw.ad %>% autoplot(main="Pronostico HW aditivo", xlab="Años/meses")
Usando forecast(aproximacion por espacios de los estados ETS)
MULTIPLICATIVO
library(forecast)
#generar el pronostico
pronosticoshw2.mul<-hw(y = Yt,
h = 12,
level = c(0.95),
seasonal = "multiplicative",
initial = "optimal")
pronosticoshw2.mul## Point Forecast Lo 95 Hi 95
## Apr 2022 111.9184 105.09260 118.7442
## May 2022 119.3664 110.67455 128.0583
## Jun 2022 122.8422 112.63807 133.0464
## Jul 2022 118.4498 107.52066 129.3789
## Aug 2022 120.0189 107.93209 132.1056
## Sep 2022 118.4260 105.56932 131.2828
## Oct 2022 114.3087 101.05402 127.5633
## Nov 2022 120.2837 105.49291 135.0745
## Dec 2022 127.9776 111.38423 144.5709
## Jan 2023 114.4558 98.88074 130.0309
## Feb 2023 113.5923 97.43199 129.7527
## Mar 2023 119.3883 101.68918 137.0873#grafico de la serie original y del pronostico
pronosticoshw2.mul %>% autoplot(main ="Pronostico aproximacion por espacios multiplicativo",xlab="Años/meses")
ADITIVO
library(forecast)
#ESTIMAR EL MODELO
modelohw.ad2<-HoltWinters(x = Yt,
seasonal= "additive",
optim.start=c(0.9,0.9,0.9))
modelohw.ad2## Holt-Winters exponential smoothing with trend and additive seasonal component.
##
## Call:
## HoltWinters(x = Yt, seasonal = "additive", optim.start = c(0.9, 0.9, 0.9))
##
## Smoothing parameters:
## alpha: 0.9142733
## beta : 0
## gamma: 1
##
## Coefficients:
## [,1]
## a 118.0938359
## b 0.1600947
## s1 -3.9032233
## s2 3.0467916
## s3 3.9273013
## s4 -1.6040517
## s5 -0.5787492
## s6 -1.8288955
## s7 -4.0808499
## s8 2.2490593
## s9 9.6275698
## s10 -2.7962367
## s11 -4.4557494
## s12 0.9861641#GENERAR EL PRONOSTICO
pronosticoshw.ad.2<-forecast(object = modelohw.ad,h = 12,level = c(0.95))
pronosticoshw.ad.2## Point Forecast Lo 95 Hi 95
## Apr 2022 114.3507 109.10096 119.6005
## May 2022 121.4608 114.34766 128.5740
## Jun 2022 122.5014 113.92039 131.0825
## Jul 2022 117.1302 107.29801 126.9623
## Aug 2022 118.3156 107.37443 129.2567
## Sep 2022 117.2255 105.27790 129.1731
## Oct 2022 115.1336 102.25799 128.0093
## Nov 2022 121.6237 107.88248 135.3648
## Dec 2022 129.1623 114.60695 143.7176
## Jan 2023 116.8985 101.57229 132.2248
## Feb 2023 115.3991 99.33889 131.4594
## Mar 2023 121.0011 104.23903 137.7632#GRAFICO D ELA SERIE ORIGINAL Y DEL PRONOSTICO
pronosticoshw.ad.2 %>% autoplot(main="Pronostico HW aditivo", xlab="Años/meses")
#GENERAR PRONOSTICOS PARA 2022 COMPLETO USANDO UN MODELO SARIMA #II. Pronóstico de series temporales, enfoque moderno (estocástico)
##Usando la metodlogia paso a paso mostrada en clases
Metodología Box-Jenkins
library(TSstudio)
library(forecast)
ts_plot(Yt,Xtitle = "Años/Meses")Identificacion
#ORDEN DE INTEGRACION
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")#VERIFICAR LOS VALORES para (p,q) & (P,Q)
Yt %>%
diff(lag = 12,diffences=D) %>%
diff(diffences=d) %>%
ts_cor(lag.max = 36)SARIMA (0,1,0)(1,1,1)
Estimación del modelo
USANDO 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.1022 -0.7232
## s.e. 0.1266 0.1135
##
## sigma^2 = 6.779: log likelihood = -351.22
## AIC=708.43 AICc=708.6 BIC=717.38
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.05202814 2.477766 1.68543 0.02782791 1.663587 0.4366552
## ACF1
## Training set 0.06837002#grafico de las raices (estabilidad/convertibilidad)
modelo.estimado %>% autoplot(type="both")+theme_solarized()
#verificacion
modelo.estimado %>% check_res(lag.max = 36)#grafico comparativo
Yt.Sarima<-modelo.estimado$fitted
Yt<-pronosticoshw$fitted
grafico.comparativo<-cbind(Yt,Yt.Sarima,Yt)
ts_plot(grafico.comparativo)Verificación de sobre ajuste/sub ajuste
SARIMA(0,1,0)(1,1,1) nuevos P-1 SARIMA(0,1,0)(0,1,1) & Q-1SARIMA(0,1,0)(1,1,0)
library(tsibble)
library(feasts)
library(fable)
library(fabletools)
library(tidyr)
library(dplyr)
#generando modelo automatico
a<-Yt %>% as_tsibble() %>%
model(arima_original=ARIMA(value ~ pdq(0, 1, 0) + PDQ(1, 1, 1)), #original
arima_010_011 = ARIMA(value ~ pdq(0, 1, 0) + PDQ(0, 1, 1)), #sin componente autoregresivo estacional
arima_010_110 = ARIMA(value ~ pdq(0, 1, 0) + PDQ(1, 1, 0)), #sin componente autoregresivo de media movil
arima_automatico=ARIMA(value,ic="bic",stepwise = FALSE) #modelo automatico
)
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>#ordenar recuerdo al criterio
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 5.93 -321. 651. 652. 666. <cpl [13]> <cpl [12]>
## 2 arima_original Orders 6.45 -325. 655. 655. 664. <cpl [12]> <cpl [12]>
## 3 arima_010_011 Orders 6.68 -328. 659. 659. 665. <cpl [0]> <cpl [12]>
## 4 arima_010_110 Orders 7.15 -331. 666. 666. 672. <cpl [12]> <cpl [0]>III. Validación Cruzada (Cross Validated)
library(forecast)
library(dplyr)
library(tsibble)
library(fable)
library(fabletools)
#convertimos
Yt<-Yt %>%
as_tsibble() %>%
rename(IVAE=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(IVAE ~ 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(IVAE ~ pdq(0~ Test -0.0247 2.98 2.06 -0.0556 1.97 0.514 0.481 0.0851EL MODELO TIENE BUEN PODER DE PREDICCION