library(readxl)
## Warning: package 'readxl' was built under R version 4.1.3
library(tidyr)
## Warning: package 'tidyr' was built under R version 4.1.3
library(dplyr)
## Warning: package 'dplyr' was built under R version 4.1.3
##
## Attaching package: 'dplyr'
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
library(forecast)
## Warning: package 'forecast' was built under R version 4.1.3
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
IVAE_03_22 <- read_excel("C:/Users/erick/Desktop/IVAE.xlsx",
col_names = FALSE, skip = 6, n_max = 10)
## New names:
## * `` -> `...1`
## * `` -> `...2`
## * `` -> `...3`
## * `` -> `...4`
## * `` -> `...5`
## * `` -> `...6`
## * `` -> `...7`
## * `` -> `...8`
## * `` -> `...9`
## * `` -> `...10`
## * `` -> `...11`
## * `` -> `...12`
## * `` -> `...13`
## * `` -> `...14`
## * `` -> `...15`
## * `` -> `...16`
## * `` -> `...17`
## * `` -> `...18`
## * `` -> `...19`
## * `` -> `...20`
## * `` -> `...21`
## * `` -> `...22`
## * `` -> `...23`
## * `` -> `...24`
## * `` -> `...25`
## * `` -> `...26`
## * `` -> `...27`
## * `` -> `...28`
## * `` -> `...29`
## * `` -> `...30`
## * `` -> `...31`
## * `` -> `...32`
## * `` -> `...33`
## * `` -> `...34`
## * `` -> `...35`
## * `` -> `...36`
## * `` -> `...37`
## * `` -> `...38`
## * `` -> `...39`
## * `` -> `...40`
## * `` -> `...41`
## * `` -> `...42`
## * `` -> `...43`
## * `` -> `...44`
## * `` -> `...45`
## * `` -> `...46`
## * `` -> `...47`
## * `` -> `...48`
## * `` -> `...49`
## * `` -> `...50`
## * `` -> `...51`
## * `` -> `...52`
## * `` -> `...53`
## * `` -> `...54`
## * `` -> `...55`
## * `` -> `...56`
## * `` -> `...57`
## * `` -> `...58`
## * `` -> `...59`
## * `` -> `...60`
## * `` -> `...61`
## * `` -> `...62`
## * `` -> `...63`
## * `` -> `...64`
## * `` -> `...65`
## * `` -> `...66`
## * `` -> `...67`
## * `` -> `...68`
## * `` -> `...69`
## * `` -> `...70`
## * `` -> `...71`
## * `` -> `...72`
## * `` -> `...73`
## * `` -> `...74`
## * `` -> `...75`
## * `` -> `...76`
## * `` -> `...77`
## * `` -> `...78`
## * `` -> `...79`
## * `` -> `...80`
## * `` -> `...81`
## * `` -> `...82`
## * `` -> `...83`
## * `` -> `...84`
## * `` -> `...85`
## * `` -> `...86`
## * `` -> `...87`
## * `` -> `...88`
## * `` -> `...89`
## * `` -> `...90`
## * `` -> `...91`
## * `` -> `...92`
## * `` -> `...93`
## * `` -> `...94`
## * `` -> `...95`
## * `` -> `...96`
## * `` -> `...97`
## * `` -> `...98`
## * `` -> `...99`
## * `` -> `...100`
## * `` -> `...101`
## * `` -> `...102`
## * `` -> `...103`
## * `` -> `...104`
## * `` -> `...105`
## * `` -> `...106`
## * `` -> `...107`
## * `` -> `...108`
## * `` -> `...109`
## * `` -> `...110`
## * `` -> `...111`
## * `` -> `...112`
## * `` -> `...113`
## * `` -> `...114`
## * `` -> `...115`
## * `` -> `...116`
## * `` -> `...117`
## * `` -> `...118`
## * `` -> `...119`
## * `` -> `...120`
## * `` -> `...121`
## * `` -> `...122`
## * `` -> `...123`
## * `` -> `...124`
## * `` -> `...125`
## * `` -> `...126`
## * `` -> `...127`
## * `` -> `...128`
## * `` -> `...129`
## * `` -> `...130`
## * `` -> `...131`
## * `` -> `...132`
## * `` -> `...133`
## * `` -> `...134`
## * `` -> `...135`
## * `` -> `...136`
## * `` -> `...137`
## * `` -> `...138`
## * `` -> `...139`
## * `` -> `...140`
## * `` -> `...141`
## * `` -> `...142`
## * `` -> `...143`
## * `` -> `...144`
## * `` -> `...145`
## * `` -> `...146`
## * `` -> `...147`
## * `` -> `...148`
## * `` -> `...149`
## * `` -> `...150`
## * `` -> `...151`
## * `` -> `...152`
## * `` -> `...153`
## * `` -> `...154`
## * `` -> `...155`
## * `` -> `...156`
## * `` -> `...157`
## * `` -> `...158`
## * `` -> `...159`
## * `` -> `...160`
data.ivae<-pivot_longer(data = IVAE_03_22[1,],names_to = "vars",cols = 2:160,values_to = "indice") %>% select("indice")
data.ivae.ts<- data.ivae %>% ts(start = c(2009,1),frequency = 12)
library(forecast)
descompocicion_aditiva<-decompose(data.ivae.ts,type = "additive")
Yt<-descompocicion_aditiva$x
TC<-descompocicion_aditiva$trend
library(forecast)
#Estimar el modelo
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.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 pronóstico:
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
#Gráfico de la serie original y del pronóstico.
PronosticosHW %>% autoplot()
library(TSstudio)
## Warning: package 'TSstudio' was built under R version 4.1.3
library(forecast)
ts_plot(Yt,Xtitle = "Años/Meses")
library(kableExtra)
## Warning: package 'kableExtra' was built under R version 4.1.3
##
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
##
## group_rows
library(magrittr)
## Warning: package 'magrittr' was built under R version 4.1.3
##
## Attaching package: 'magrittr'
## The following object is masked from 'package:tidyr':
##
## extract
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()
Ordenes de Integración
|
|
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")
Yt %>%
diff(lag = 12,diffences=D) %>%
diff(diffences=d) %>%
ts_cor(lag.max = 36)
library(forecast)
library(ggthemes)
## Warning: package 'ggthemes' was built under R version 4.1.3
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
modelo_estimado %>% autoplot(type="both")+theme_solarized()
modelo_estimado %>% check_res(lag.max = 36)
## Warning: Ignoring 36 observations
## Warning: Ignoring 34 observations
## Warning: Ignoring 4 observations
Yt_Sarima<-modelo_estimado$fitted
Yt_HW<-PronosticosHW$fitted
grafico_comparativo<-cbind(Yt,Yt_Sarima,Yt_HW)
ts_plot(grafico_comparativo)
library(tsibble)
## Warning: package 'tsibble' was built under R version 4.1.3
##
## Attaching package: 'tsibble'
## The following objects are masked from 'package:base':
##
## intersect, setdiff, union
library(feasts)
## Warning: package 'feasts' was built under R version 4.1.3
## Loading required package: fabletools
## Warning: package 'fabletools' was built under R version 4.1.3
##
## Attaching package: 'fabletools'
## The following objects are masked from 'package:forecast':
##
## accuracy, forecast
library(fable)
## Warning: package 'fable' was built under R version 4.1.3
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 6.13 -347. 702. 702. 714. <cpl [1]> <cpl [12]>
## 2 arima_010_011 Orders 6.72 -352. 707. 707. 713. <cpl [0]> <cpl [12]>
## 3 arima_original Orders 6.78 -351. 708. 709. 717. <cpl [12]> <cpl [12]>
## 4 arima_010_110 Orders 7.99 -361. 726. 726. 731. <cpl [12]> <cpl [0]>
library(forecast)
ma2_12 <- ma(data.ivae.ts, 12, centre = TRUE)
library(magrittr)
Yt <- data.ivae.ts
Tt <- ma2_12
SI <- Yt - Tt
St <- tapply(SI, cycle(SI), mean, na.rm = TRUE)
#Los factores estacionales deben sumar "0" en el modelo aditivo
St <- St - sum(St) / 12
#Generar la serie de factores para cada valor de la serie original
St <-
rep(St, len = length(Yt)) %>% ts(start = c(2009, 1), frequency = 12)
It<-Yt-Tt-St
library(tsibble)
library(feasts)
library(ggplot2)
## Warning: package 'ggplot2' was built under R version 4.1.3
Yt %>% as_tsibble() %>%
model(
classical_decomposition(value, type = "additive")
) %>%
components() %>%
autoplot() +
labs(title = "Descomposición Clásica Aditiva, IVAE")+xlab("Años/Meses")
## Warning: Removed 6 row(s) containing missing values (geom_path).
library(forecast)
library(dplyr)
library(tsibble)
library(fable)
library(fabletools)
Yt<-Yt %>% as_tsibble() %>% rename(IVAE=value)
data.cross.validation<-Yt %>%
as_tsibble() %>%
stretch_tsibble(.init = 60,.step = 1)
TSCV<-data.cross.validation %>%
model(ARIMA(IVAE ~ pdq(0, 1, 0) + PDQ(1, 1, 1))) %>%
forecast(h=1) %>% accuracy(Yt)
## Warning: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing.
## 1 observation is missing at 2022 abr.
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.0506 2.94 2.09 0.00243 2.00 0.542 0.494 0.155
#El MAPE es de 2.004%, quiere decir que de cada 100 valores pronósticados, la distancia o el error cometido es de 2.004%
PIB TRIMESTRAL
PIBTRI <- read_excel("C:/Users/erick/Desktop/pibrimestral.xlsx",
col_types=c("skip", "numeric"))
data.pib.ts<- ts(data=PIBTRI, start = c(2009,1), frequency=4)
library(forecast)
descompocicion_aditiva2<-decompose(data.pib.ts,type = "additive")
Yt2<-descompocicion_aditiva2$x
TC2<-descompocicion_aditiva2$trend
library(forecast)
#Estimar el modelo
ModeloHW2<-HoltWinters(x = Yt2,
seasonal = "multiplicative",
optim.start = c(0.9,0.9,0.9))
ModeloHW2
## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
##
## Call:
## HoltWinters(x = Yt2, seasonal = "multiplicative", optim.start = c(0.9, 0.9, 0.9))
##
## Smoothing parameters:
## alpha: 0.8742302
## beta : 0
## gamma: 0
##
## Coefficients:
## [,1]
## a 7429.8128051
## b 48.3963750
## s1 0.9751646
## s2 1.0058150
## s3 0.9783521
## s4 1.0406683
#Generar el pronóstico:
PronosticosHW2<-forecast(object = ModeloHW2,h = 4,level = c(0.95))
PronosticosHW2
## Point Forecast Lo 95 Hi 95
## 2022 Q1 7292.485 6847.413 7737.557
## 2022 Q2 7570.373 6930.977 8209.768
## 2022 Q3 7411.019 6646.235 8175.804
## 2022 Q4 7933.429 7119.875 8746.982
#Gráfico de la serie original y del pronóstico.
PronosticosHW2 %>% autoplot()
library(TSstudio)
library(forecast)
ts_plot(Yt2,Xtitle = "Años/Meses")
library(kableExtra)
library(magrittr)
d2<-ndiffs(Yt2)
D2<-nsdiffs(Yt2)
ordenes_integracion2<-c(d2,D2)
names(ordenes_integracion2)<-c("d2","D2")
ordenes_integracion2 %>% kable(caption = "Ordenes de Integración") %>% kable_material()
Ordenes de Integración
|
|
x
|
|
d2
|
1
|
|
D2
|
1
|
#Gráfico de la serie diferenciada
Yt2 %>%
diff(lag = 4,diffences=D2) %>%
diff(diffences=d2) %>%
ts_plot(title = "Yt estacionaria")
Yt2 %>%
diff(lag = 4,diffences=D2) %>%
diff(diffences=d2) %>%
ts_cor(lag.max = 12)
library(forecast)
library(ggthemes)
modelo_estimado2 <- Yt2 %>%
Arima(order = c(0, 1, 0),
seasonal = c(1, 1, 2))
summary(modelo_estimado2)
## Series: .
## ARIMA(0,1,0)(1,1,2)[4]
##
## Coefficients:
## sar1 sma1 sma2
## -0.3494 -0.5042 0.0640
## s.e. 0.5359 0.5323 0.3901
##
## sigma^2 = 79136: log likelihood = -331.65
## AIC=671.31 AICc=672.26 BIC=678.71
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 6.615822 258.7679 121.0037 -0.002908518 2.055323 0.3682692
## ACF1
## Training set -0.1013843
modelo_estimado2 %>% autoplot(type="both")+theme_solarized()
modelo_estimado2 %>% check_res(lag.max = 12)
## Warning: Ignoring 12 observations
## Warning: Ignoring 10 observations
## Warning: Ignoring 4 observations
Yt_Sarima2<-modelo_estimado2$fitted
Yt_HW2<-PronosticosHW2$fitted
grafico_comparativo2<-cbind(Yt2,Yt_Sarima2,Yt_HW2)
ts_plot(grafico_comparativo2)
library(tsibble)
library(feasts)
library(fable)
library(fabletools)
library(tidyr)
library(dplyr)
a2<-Yt2 %>% as_tsibble() %>%
model(arima_original2=ARIMA(value ~ pdq(0, 1, 0) + PDQ(1, 1, 2)),
arima_010_0112 = ARIMA(value ~ pdq(0, 1, 0) + PDQ(0, 1, 1)),
arima_010_1102 = ARIMA(value ~ pdq(0, 1, 0) + PDQ(1, 1, 1)),
arima_automatico2=ARIMA(value,ic="bic",stepwise = FALSE)
)
print(a2)
## # A mable: 1 x 4
## arima_original2 arima_010_0112 arima_010_1102
## <model> <model> <model>
## 1 <ARIMA(0,1,0)(1,1,2)[4]> <ARIMA(0,1,0)(0,1,1)[4]> <ARIMA(0,1,0)(1,1,1)[4]>
## # ... with 1 more variable: arima_automatico2 <model>
a2 %>% pivot_longer(everything(), names_to = "Model name",
values_to = "Orders") %>% glance() %>%
arrange(AICc) ->tabla2
tabla2
## # 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_010_0112 Orders 75930. -332. 668. 669. 672. <cpl [0]> <cpl [4]>
## 2 arima_010_1102 Orders 77377. -332. 669. 670. 675. <cpl [4]> <cpl [4]>
## 3 arima_original2 Orders 79136. -332. 671. 672. 679. <cpl [4]> <cpl [8]>
## 4 arima_automatico2 Orders 62738. -334. 676. 677. 683. <cpl [1]> <cpl [4]>
library(forecast)
ma2_122 <- ma(data.pib.ts, 4, centre = TRUE)
library(magrittr)
Yt2 <- data.pib.ts
Tt2 <- ma2_122
SI2 <- Yt2 - Tt2
St2 <- tapply(SI2, cycle(SI2), mean, na.rm = TRUE)
#Los factores estacionales deben sumar "0" en el modelo aditivo
St2 <- St2 - sum(St2) / 4
#Generar la serie de factores para cada valor de la serie original
St2 <-
rep(St2, len = length(Yt2)) %>% ts(start = c(2009, 1), frequency = 4)
It2<-Yt2-Tt2-St2
library(tsibble)
library(feasts)
library(ggplot2)
Yt2 %>% as_tsibble() %>%
model(
classical_decomposition(value, type = "additive")
) %>%
components() %>%
autoplot() +
labs(title = "Descomposición Clásica Aditiva, IVAE")+xlab("Años/Meses")
## Warning: Removed 2 row(s) containing missing values (geom_path).
library(forecast)
library(dplyr)
library(tsibble)
library(fable)
library(fabletools)
Yt2<-Yt2 %>% as_tsibble() %>% rename(PIB=value)
data.cross.validation2<-Yt2 %>%
as_tsibble() %>%
stretch_tsibble(.init = 30,.step = 1)
TSCV2<-data.cross.validation2 %>%
model(ARIMA(PIB ~ pdq(0, 1, 0) + PDQ(1, 1, 2))) %>%
forecast(h=4) %>% accuracy(Yt2)
## Warning: It looks like you're trying to fully specify your ARIMA model but have not said if a constant should be included.
## You can include a constant using `ARIMA(y~1)` to the formula or exclude it by adding `ARIMA(y~0)`.
## Warning: 1 error encountered for ARIMA(PIB ~ pdq(0, 1, 0) + PDQ(1, 1, 2))
## [1] Could not find an appropriate ARIMA model.
## This is likely because automatic selection does not select models with characteristic roots that may be numerically unstable.
## For more details, refer to https://otexts.com/fpp3/arima-r.html#plotting-the-characteristic-roots
## Warning: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing.
## 4 observations are missing between 2022 Q1 and 2022 Q4
print(TSCV2)
## # 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(PIB ~ pdq(0, 1, 0~ Test 136. 725. 409. 1.50 6.34 1.24 1.53 0.681
#El MAPE es de 6.33%, quiere decir que de cada 100 valores pronósticados, la distancia o el error cometido es de 6.33%
IPC
IPC <- read_excel("C:/Users/erick/Desktop/IPC.xlsx",
col_types=c("skip", "numeric"))
data.ipc.ts<- ts(data=IPC, start = c(2009,1), frequency=12)
library(forecast)
descompocicion_aditiva3<-decompose(data.ipc.ts,type = "additive")
Yt3<-descompocicion_aditiva3$x
TC3<-descompocicion_aditiva3$trend
library(forecast)
#Estimar el modelo
ModeloHW3<-HoltWinters(x = Yt3,
seasonal = "multiplicative",
optim.start = c(0.9,0.9,0.9))
ModeloHW3
## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
##
## Call:
## HoltWinters(x = Yt3, seasonal = "multiplicative", optim.start = c(0.9, 0.9, 0.9))
##
## Smoothing parameters:
## alpha: 0.8353873
## beta : 0.07612657
## gamma: 1
##
## Coefficients:
## [,1]
## a 122.7690242
## b 0.4801476
## s1 1.0043291
## s2 1.0021040
## s3 0.9986014
## s4 0.9975846
## s5 0.9989586
## s6 0.9984650
## s7 0.9959074
## s8 0.9975598
## s9 1.0007205
## s10 1.0041558
## s11 1.0048925
## s12 1.0053839
#Generar el pronóstico:
PronosticosHW3<-forecast(object = ModeloHW3,h = 12,level = c(0.95))
PronosticosHW3
## Point Forecast Lo 95 Hi 95
## Jun 2022 123.7827 122.6365 124.9289
## Jul 2022 123.9896 122.4499 125.5294
## Aug 2022 124.0357 122.1459 125.9256
## Sep 2022 124.3884 122.1648 126.6120
## Oct 2022 125.0394 122.4870 127.5918
## Nov 2022 125.4570 122.5842 128.3298
## Dec 2022 125.6139 122.4280 128.7997
## Jan 2023 126.3012 122.7897 129.8128
## Feb 2023 127.1819 123.3369 131.0269
## Mar 2023 128.1007 123.9171 132.2842
## Apr 2023 128.6771 124.1616 133.1927
## May 2023 129.2228 123.9338 134.5118
#Gráfico de la serie original y del pronóstico.
PronosticosHW3 %>% autoplot()
library(TSstudio)
library(forecast)
ts_plot(Yt3,Xtitle = "Años/Meses")
library(kableExtra)
library(magrittr)
d3<-ndiffs(Yt3)
D3<-nsdiffs(Yt3)
ordenes_integracion3<-c(d3,D3)
names(ordenes_integracion3)<-c("d3","D3")
ordenes_integracion3 %>% kable(caption = "Ordenes de Integración") %>% kable_material()
Ordenes de Integración
|
|
x
|
|
d3
|
1
|
|
D3
|
0
|
#Gráfico de la serie diferenciada
Yt3 %>%
diff(lag = 12,diffences=D3) %>%
diff(diffences=d3) %>%
ts_plot(title = "Yt estacionaria")
Yt3 %>%
diff(lag = 12,diffences=D3) %>%
diff(diffences=d3) %>%
ts_cor(lag.max = 36)
library(forecast)
library(ggthemes)
modelo_estimado3 <- Yt3 %>%
Arima(order = c(1, 1, 1),
seasonal = c(1, 0, 3))
summary(modelo_estimado3)
## Series: .
## ARIMA(1,1,1)(1,0,3)[12]
##
## Coefficients:
## ar1 ma1 sar1 sma1 sma2 sma3
## 0.9289 -0.7894 0.9905 -0.9135 -0.1630 0.1131
## s.e. 0.0840 0.1307 0.0416 0.1244 0.1065 0.0880
##
## sigma^2 = 0.2175: log likelihood = -104.18
## AIC=222.37 AICc=223.1 BIC=243.89
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.05115671 0.4561007 0.3004523 0.04545039 0.2756045 0.1628142
## ACF1
## Training set 0.07731549
modelo_estimado3 %>% autoplot(type="both")+theme_solarized()
modelo_estimado3 %>% check_res(lag.max = 36)
## Warning: Ignoring 36 observations
## Warning: Ignoring 34 observations
## Warning: Ignoring 4 observations
Yt_Sarima3<-modelo_estimado3$fitted
Yt_HW3<-PronosticosHW3$fitted
grafico_comparativo3<-cbind(Yt3,Yt_Sarima3,Yt_HW3)
ts_plot(grafico_comparativo3)
library(tsibble)
library(feasts)
library(fable)
library(fabletools)
library(tidyr)
library(dplyr)
a3<-Yt3 %>% as_tsibble() %>%
model(arima_original3=ARIMA(value ~ pdq(1, 1, 1) + PDQ(1, 0, 3)),
arima_010_0113 = ARIMA(value ~ pdq(0, 1, 0) + PDQ(0, 1, 3)),
arima_010_1103 = ARIMA(value ~ pdq(0, 1, 0) + PDQ(1, 1, 2)),
arima_automatico3=ARIMA(value,ic="bic",stepwise = FALSE)
)
print(a3)
## # A mable: 1 x 4
## arima_original3 arima_010_0113
## <model> <model>
## 1 <ARIMA(1,1,1)(1,0,3)[12] w/ drift> <ARIMA(0,1,0)(0,1,3)[12]>
## # ... with 2 more variables: arima_010_1103 <model>, arima_automatico3 <model>
a3 %>% pivot_longer(everything(), names_to = "Model name",
values_to = "Orders") %>% glance() %>%
arrange(AICc) ->tabla3
tabla3
## # 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_original3 Orders 0.221 -103. 222. 223. 246. <cpl [13]> <cpl>
## 2 arima_automatico3 Orders 0.225 -105. 223. 224. 241. <cpl [27]> <cpl [0]>
## 3 arima_010_0113 Orders 0.225 -114. 235. 236. 247. <cpl [0]> <cpl>
## 4 arima_010_1103 Orders 0.230 -114. 236. 236. 248. <cpl [12]> <cpl>
library(forecast)
ma2_123 <- ma(data.ipc.ts, 12, centre = TRUE)
library(magrittr)
Yt3 <- data.ipc.ts
Tt3 <- ma2_123
SI3 <- Yt3 - Tt3
St3 <- tapply(SI3, cycle(SI3), mean, na.rm = TRUE)
#Los factores estacionales deben sumar "0" en el modelo aditivo
St3 <- St3 - sum(St3) / 12
#Generar la serie de factores para cada valor de la serie original
St3 <-
rep(St3, len = length(Yt3)) %>% ts(start = c(2009, 1), frequency = 12)
It3<-Yt3-Tt3-St3
library(tsibble)
library(feasts)
library(ggplot2)
Yt3 %>% as_tsibble() %>%
model(
classical_decomposition(value, type = "additive")
) %>%
components() %>%
autoplot() +
labs(title = "Descomposición Clásica Aditiva, IVAE")+xlab("Años/Meses")
## Warning: Removed 6 row(s) containing missing values (geom_path).
library(forecast)
library(dplyr)
library(tsibble)
library(fable)
library(fabletools)
Yt3<-Yt3 %>% as_tsibble() %>% rename(IPC=value)
data.cross.validation3<-Yt3 %>%
as_tsibble() %>%
stretch_tsibble(.init = 60,.step = 1)
TSCV3<-data.cross.validation3 %>%
model(ARIMA(IPC ~ pdq(0, 1, 0) + PDQ(0, 1, 0))) %>%
forecast(h=12) %>% accuracy(Yt3)
## Warning: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing.
## 12 observations are missing between 2022 jun. and 2023 may.
print(TSCV3)
## # 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(IPC ~ pdq(0, 1, 0~ Test 0.415 1.89 1.35 0.345 1.19 0.732 0.695 0.881
#El MAPE es de 1.19%, quiere decir que de cada 100 valores pronósticados, la distancia o el error cometido es de 1.19%
HIDROCARBUROS
hidro <- read_excel("C:/Users/erick/Desktop/hidrocarburos.xlsx",
col_types=c("skip", "numeric"))
data.hidro.ts<- ts(data=hidro, start = c(2017,1), frequency=12)
library(forecast)
descompocicion_aditiva4<-decompose(data.hidro.ts,type = "additive")
Yt4<-descompocicion_aditiva4$x
TC4<-descompocicion_aditiva4$trend
library(forecast)
#Estimar el modelo
ModeloHW4<-HoltWinters(x = Yt4,
seasonal = "multiplicative",
optim.start = c(0.9,0.9,0.9))
ModeloHW4
## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
##
## Call:
## HoltWinters(x = Yt4, seasonal = "multiplicative", optim.start = c(0.9, 0.9, 0.9))
##
## Smoothing parameters:
## alpha: 0.1273406
## beta : 0
## gamma: 0.4793705
##
## Coefficients:
## [,1]
## a 2.465101e+05
## b 1.239453e+03
## s1 9.247101e-01
## s2 8.984179e-01
## s3 8.791852e-01
## s4 8.941624e-01
## s5 9.460044e-01
## s6 1.049825e+00
## s7 9.469010e-01
## s8 9.667736e-01
## s9 9.762889e-01
## s10 1.117116e+00
## s11 9.509380e-01
## s12 8.557828e-01
#Generar el pronóstico:
PronosticosHW4<-forecast(object = ModeloHW4,h = 12,level = c(0.95))
PronosticosHW4
## Point Forecast Lo 95 Hi 95
## Jun 2022 229096.5 185614.4 272578.5
## Jul 2022 223696.1 179110.4 268281.8
## Aug 2022 219997.1 174364.1 265630.2
## Sep 2022 224853.1 177997.3 271708.9
## Oct 2022 239062.2 190615.5 287508.9
## Nov 2022 266599.7 215810.4 317389.1
## Dec 2022 241636.1 191288.5 291983.7
## Jan 2023 247905.6 196238.0 299573.1
## Feb 2023 251555.6 198735.7 304375.5
## Mar 2023 289226.4 232632.2 345820.6
## Apr 2023 247380.9 193379.5 301382.3
## May 2023 223687.5 193456.8 253918.2
#Gráfico de la serie original y del pronóstico.
PronosticosHW4 %>% autoplot()
library(TSstudio)
library(forecast)
ts_plot(Yt4,Xtitle = "Años/Meses")
library(kableExtra)
library(magrittr)
d4<-ndiffs(Yt4)
D4<-nsdiffs(Yt4)
ordenes_integracion4<-c(d4,D4)
names(ordenes_integracion4)<-c("d4","D4")
ordenes_integracion4 %>% kable(caption = "Ordenes de Integración") %>% kable_material()
Ordenes de Integración
|
|
x
|
|
d4
|
0
|
|
D4
|
0
|
#Gráfico de la serie diferenciada
Yt4 %>%
diff(lag = 12,diffences=D4) %>%
diff(diffences=d4) %>%
ts_plot(title = "Yt estacionaria")
Yt4 %>%
diff(lag = 12,diffences=D4) %>%
diff(diffences=d4) %>%
ts_cor(lag.max = 36)
library(forecast)
library(ggthemes)
modelo_estimado4 <- Yt4 %>%
Arima(order = c(1, 0, 1),
seasonal = c(1, 0, 0))
summary(modelo_estimado4)
## Series: .
## ARIMA(1,0,1)(1,0,0)[12] with non-zero mean
##
## Coefficients:
## ar1 ma1 sar1 mean
## 0.6883 -0.5363 -0.2097 215368.837
## s.e. 0.2890 0.3277 0.1494 4910.911
##
## sigma^2 = 1.058e+09: log likelihood = -765.8
## AIC=1541.6 AICc=1542.62 BIC=1552.48
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -237.9014 31510 25506.37 -2.377165 12.31299 0.6275001 -0.02020648
modelo_estimado4 %>% autoplot(type="both")+theme_solarized()
modelo_estimado4 %>% check_res(lag.max = 36)
## Warning: Ignoring 36 observations
## Warning: Ignoring 34 observations
## Warning: Ignoring 4 observations
Yt_Sarima4<-modelo_estimado4$fitted
Yt_HW4<-PronosticosHW4$fitted
grafico_comparativo4<-cbind(Yt4,Yt_Sarima4,Yt_HW4)
ts_plot(grafico_comparativo4)
library(tsibble)
library(feasts)
library(fable)
library(fabletools)
library(tidyr)
library(dplyr)
a4<-Yt4 %>% as_tsibble() %>%
model(arima_original4=ARIMA(value ~ pdq(1, 0, 1) + PDQ(1, 0, 0)),
arima_010_0114 = ARIMA(value ~ pdq(1, 0, 1) + PDQ(0, 0, 0)),
arima_010_1104 = ARIMA(value ~ pdq(1, 0, 1) + PDQ(1, 0, 0)),
arima_automatico4=ARIMA(value,ic="bic",stepwise = FALSE)
)
print(a4)
## # A mable: 1 x 4
## arima_original4 arima_010_0114
## <model> <model>
## 1 <ARIMA(1,0,1)(1,0,0)[12] w/ mean> <ARIMA(1,0,1) w/ mean>
## # ... with 2 more variables: arima_010_1104 <model>, arima_automatico4 <model>
a4 %>% pivot_longer(everything(), names_to = "Model name",
values_to = "Orders") %>% glance() %>%
arrange(AICc) ->tabla4
tabla4
## # 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_automatico4 Orders 1.11e9 -769. 1541. 1541. 1545. <cpl> <cpl>
## 2 arima_original4 Orders 1.06e9 -766. 1542. 1543. 1552. <cpl> <cpl>
## 3 arima_010_1104 Orders 1.06e9 -766. 1542. 1543. 1552. <cpl> <cpl>
## 4 arima_010_0114 Orders 1.11e9 -768. 1543. 1544. 1552. <cpl> <cpl>
library(forecast)
ma2_124 <- ma(data.hidro.ts, 12, centre = TRUE)
library(magrittr)
Yt4 <- data.hidro.ts
Tt4 <- ma2_124
SI4 <- Yt4 - Tt4
St4 <- tapply(SI4, cycle(SI4), mean, na.rm = TRUE)
#Los factores estacionales deben sumar "0" en el modelo aditivo
St4 <- St4 - sum(St4) / 12
#Generar la serie de factores para cada valor de la serie original
St4 <-
rep(St4, len = length(Yt4)) %>% ts(start = c(2017, 1), frequency = 12)
It4<-Yt4-Tt4-St4
library(tsibble)
library(feasts)
library(ggplot2)
Yt4 %>% as_tsibble() %>%
model(
classical_decomposition(value, type = "additive")
) %>%
components() %>%
autoplot() +
labs(title = "Descomposición Clásica Aditiva, IVAE")+xlab("Años/Meses")
## Warning: Removed 6 row(s) containing missing values (geom_path).
library(forecast)
library(dplyr)
library(tsibble)
library(fable)
library(fabletools)
Yt4<-Yt4 %>% as_tsibble() %>% rename(HIDRO=value)
data.cross.validation4<-Yt4 %>%
as_tsibble() %>%
stretch_tsibble(.init = 60,.step = 1)
TSCV4<-data.cross.validation4 %>%
model(ARIMA(HIDRO ~ pdq(1, 0, 1) + PDQ(1, 0, 0))) %>%
forecast(h=12) %>% accuracy(Yt4)
## Warning: The future dataset is incomplete, incomplete out-of-sample data will be treated as missing.
## 12 observations are missing between 2022 jun. and 2023 may.
print(TSCV4)
## # 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(HIDRO ~ pdq(1,~ Test 15649. 49078. 40075. 3.09 16.8 0.986 0.990 0.234
#El MAPE es de 16.78%, quiere decir que de cada 100 valores pronósticados, la distancia o el error cometido es de 16.78%