Metodología Box-Jenkins
Importacion de la data
library(readxl)
library(forecast)
PIBTPC<- read_excel("~/Econometria/tableExportPIB.xlsx",
col_types = c("skip", "numeric"), skip = 5)
## Serie
serie.PIB.ts<- ts(data = PIBTPC, start = c(2009,1),frequency = 4)
Yt_1 <- serie.PIB.ts
serie.PIB.ts %>%
autoplot(main = "PIB trimestral precios corrientes SLV 2009-2021",
xlab = "años/ trimestres",
ylab = "Indice")
Generacion de Grafica general
library(TSstudio)
library(forecast)
ts_plot(Yt_1,Xtitle= "Años/Meses")
Orden de Integración
library(kableExtra)
library(magrittr)
d_1<-ndiffs(Yt_1)
D_1<-nsdiffs(Yt_1)
ordenes_integracion<-c(d_1,D_1)
names(ordenes_integracion)<-c("d","D")
ordenes_integracion %>% kable(caption = "Ordenes de Integracion") %>% kable_material()
| x | |
|---|---|
| d | 1 |
| D | 1 |
#Gráfico de la serie diferenciada
Yt_1 %>%
diff(lag=4,diffences=D_1) %>%
diff(diffences=d_1) %>%
ts_plot(title= "Yt estacionaria")
Se construyen las funciones de autocorrelación simple y parcial.
Yt_1 %>%
diff(lag=4,diffences=D) %>%
diff(diffences=d) %>%
ts_cor(lag.max = 12)
Usando forecast
library(forecast)
library(ggthemes)
modelo_estimado_PIB<-Yt_1 %>%
Arima(order = c(0, 1, 0),
seasonal = c(1, 1, 1))
summary(modelo_estimado_PIB)
## Series: .
## ARIMA(0,1,0)(1,1,1)[4]
##
## Coefficients:
## sar1 sma1
## -0.3980 -0.4550
## s.e. 0.3732 0.3688
##
## sigma^2 = 77377: log likelihood = -331.67
## AIC=669.34 AICc=669.9 BIC=674.89
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 6.530071 258.7683 120.9372 -0.004468918 2.054645 0.3680668
## ACF1
## Training set -0.09806694
modelo_estimado_PIB %>% autoplot(type="both")+theme_solarized()
modelo_estimado_PIB %>% check_res(lag.max = 4)
## Warning: Ignoring 4 observations
## Warning: Ignoring 4 observations
## Warning: Ignoring 2 observations
Importacion de PronosticosHW de la practica anterior # Estimación del modelo
library(forecast)
#Estimar el modelo
modeloHW_PIB<-HoltWinters(x = Yt_1,
seasonal = "multiplicative",
optim.start = c(0.9,0.9,0.9))
modeloHW_PIB
## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
##
## Call:
## HoltWinters(x = Yt_1, 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 pronostico:
pronosticosMW_PIB<-forecast(object = modeloHW_PIB,h=4, level=c(0.95))
pronosticosMW_PIB
## 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
Yt_Sarima_PIB<-modelo_estimado_PIB$fitted
Yt_HW_PIB<-pronosticosMW_PIB$fitted
grafico_comparativo_1<-cbind(Yt_1, Yt_Sarima_PIB,Yt_HW_PIB)
ts_plot(grafico_comparativo_1)
library(tsibble)
##
## Attaching package: 'tsibble'
## The following objects are masked from 'package:base':
##
## intersect, setdiff, union
library(feasts)
## Loading required package: fabletools
##
## Attaching package: 'fabletools'
## The following objects are masked from 'package:forecast':
##
## accuracy, forecast
library(fable)
library(fabletools)
library(tidyr)
##
## Attaching package: 'tidyr'
## The following object is masked from 'package:magrittr':
##
## extract
library(dplyr)
##
## Attaching package: 'dplyr'
## The following object is masked from 'package:kableExtra':
##
## group_rows
## The following objects are masked from 'package:stats':
##
## filter, lag
## The following objects are masked from 'package:base':
##
## intersect, setdiff, setequal, union
a<-Yt_1 %>% 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)[4]> <ARIMA(0,1,0)(0,1,1)[4]> <ARIMA(0,1,0)(1,1,0)[4]>
## # … 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 × 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_011 Orders 75930. -332. 668. 669. 672. <cpl [0]> <cpl [4]>
## 2 arima_010_110 Orders 78594. -333. 669. 670. 673. <cpl [4]> <cpl [0]>
## 3 arima_original Orders 77377. -332. 669. 670. 675. <cpl [4]> <cpl [4]>
## 4 arima_automatico Orders 62738. -334. 676. 677. 683. <cpl [1]> <cpl [4]>
library(forecast)
library(dplyr)
library(tsibble)
library(fable)
library(fabletools)
#convertimos
Yt<-Yt_1 %>%
as_tsibble() %>%
rename(PIB=value)
#generador de la dta (roja)
data.cross.validation.PIB<-Yt %>%
as_tsibble() %>%
stretch_tsibble(.init = 20,.step = 1)
#generador de pronostico (verde) y simulacion
TSCV<-data.cross.validation.PIB %>%
model(ARIMA(PIB ~ 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 Q1
print(TSCV)
## # A tibble: 1 × 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, … Test 30.1 370. 197. 0.228 3.21 0.600 0.780 0.0422
Importar la data
library(readxl)
library(forecast)
IPC<- read_excel("~/Econometria/tableExportIPC.xlsx",col_types=c("skip","numeric"),
skip = 5)
## Serie
serie.IPC.ts<- ts(data =IPC, start = c(2009,1),
frequency = 12)
Yt.2 <- serie.IPC.ts
serie.IPC.ts %>%
autoplot(main = "IPC general SLV 2009-2022 (mayo)",
xlab = "años/ trimestres",
ylab = "Indice")
Generacion de Grafica general
library(TSstudio)
library(forecast)
ts_plot(Yt.2,Xtitle= "Años/Meses")
Orden de Integración
library(kableExtra)
library(magrittr)
d_3<-ndiffs(Yt.2)
D_3<-nsdiffs(Yt.2)
ordenes_integracion<-c(d_3,D_3)
names(ordenes_integracion)<-c("d","D")
ordenes_integracion %>% kable(caption = "Ordenes de Integracion") %>% kable_material()
| x | |
|---|---|
| d | 1 |
| D | 0 |
#Gráfico de la serie diferenciada
Yt.2 %>%
diff(lag=12,diffences=D_3) %>%
diff(diffences=d_3) %>%
ts_plot(title= "Yt estacionaria")
Se construyen las funciones de autocorrelación simple y parcial.
Yt.2 %>%
diff(lag=12,diffences=D_2) %>%
diff(diffences=d_2) %>%
ts_cor(lag.max = 36)
Usando forecast
library(forecast)
library(ggthemes)
modelo_estimado_IPC<-Yt.2 %>%
Arima(order = c(1, 1, 1),
seasonal = c(1, 0, 2))
summary(modelo_estimado_IPC)
## Series: .
## ARIMA(1,1,1)(1,0,2)[12]
##
## Coefficients:
## ar1 ma1 sar1 sma1 sma2
## 0.9414 -0.7965 0.0702 0.0478 -0.0678
## s.e. 0.0810 0.1355 0.6798 0.6773 0.1110
##
## sigma^2 = 0.224: log likelihood = -105.16
## AIC=222.31 AICc=222.86 BIC=240.76
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.0566041 0.464391 0.3076996 0.05030787 0.2821278 0.1667415
## ACF1
## Training set 0.08593583
#grafico de las raices (estabilidad/convertibilidad)
modelo_estimado_IPC %>% autoplot(type="both")+theme_solarized()
#verificacion
modelo_estimado_IPC %>% check_res(lag.max = 36)
## Warning: Ignoring 36 observations
## Warning: Ignoring 34 observations
## Warning: Ignoring 4 observations
library(forecast)
#Estimar el modelo
modeloHW_IPC<-HoltWinters(x = Yt.2,
seasonal = "multiplicative",
optim.start = c(0.9,0.9,0.9))
modeloHW_IPC
## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
##
## Call:
## HoltWinters(x = Yt.2, 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
pronosticosMW_IPC<-forecast(object = modeloHW_IPC,h=12, level=c(0.95,0.70))
pronosticosMW_IPC
## Point Forecast Lo 95 Hi 95 Lo 70 Hi 70
## Jun 2022 123.7827 122.6365 124.9289 123.1766 124.3888
## Jul 2022 123.9896 122.4499 125.5294 123.1754 124.8038
## Aug 2022 124.0357 122.1459 125.9256 123.0364 125.0351
## Sep 2022 124.3884 122.1648 126.6120 123.2126 125.5643
## Oct 2022 125.0394 122.4870 127.5918 123.6897 126.3891
## Nov 2022 125.4570 122.5842 128.3298 123.9379 126.9762
## Dec 2022 125.6139 122.4280 128.7997 123.9292 127.2985
## Jan 2023 126.3012 122.7897 129.8128 124.4443 128.1582
## Feb 2023 127.1819 123.3369 131.0269 125.1487 129.2152
## Mar 2023 128.1007 123.9171 132.2842 125.8884 130.3129
## Apr 2023 128.6771 124.1616 133.1927 126.2893 131.0650
## May 2023 129.2228 123.9338 134.5118 126.4259 132.0196
Yt_Sarima_IPC<-modelo_estimado_IPC$fitted
Yt_HW_IPC<-pronosticosMW_IPC$fitted
grafico_comparativo_2_IPC<-cbind(Yt.2, Yt_Sarima_IPC,Yt_HW_IPC)
ts_plot(grafico_comparativo_2_IPC)
library(tsibble)
library(feasts)
library(fable)
library(fabletools)
library(tidyr)
library(dplyr)
#generando modelo automatico
a_3<-Yt.2 %>% as_tsibble() %>%
model(arima_original=ARIMA(value ~ pdq(1, 1, 1) + PDQ(1, 0, 2)),
arima_010_011 = ARIMA(value ~ pdq(1, 1, 1) + PDQ(0, 0, 2)),
arima_010_110 = ARIMA(value ~ pdq(1, 1, 1) + PDQ(1, 0, 1)),
arima_automatico=ARIMA(value,ic="bic",stepwise = FALSE)
)
print(a_3)
## # A mable: 1 x 4
## arima_original arima_010_011
## <model> <model>
## 1 <ARIMA(1,1,1)(1,0,2)[12] w/ drift> <ARIMA(1,1,1)(0,0,2)[12] w/ drift>
## # … with 2 more variables: arima_010_110 <model>, arima_automatico <model>
#ordenar
a_3 %>% pivot_longer(everything(), names_to = "Model name",
values_to = "Orders") %>% glance() %>%
arrange(AICc) ->tabla_3
tabla_3
## # A tibble: 4 × 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_011 Orders 0.218 -103. 218. 218. 236. <cpl [1]> <cpl [25]>
## 2 arima_010_110 Orders 0.219 -103. 218. 219. 237. <cpl [13]> <cpl [13]>
## 3 arima_original Orders 0.219 -103. 220. 221. 241. <cpl [13]> <cpl [25]>
## 4 arima_automatico Orders 0.225 -105. 223. 224. 241. <cpl [27]> <cpl [0]>
library(forecast)
library(dplyr)
library(tsibble)
library(fable)
library(fabletools)
#convertimos
Yt<-Yt.2 %>% as_tsibble() %>% rename(IPC=value)
#generador de la dta (roja)
data.cross.validation.IPC<-Yt %>%
as_tsibble() %>%
stretch_tsibble(.init = 60,.step = 1)
#generador de pronostico (verde) y simulacion
TSCV_3<-data.cross.validation.IPC %>%
model(ARIMA(IPC ~ 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 jun.
print(TSCV_3)
## # A tibble: 1 × 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,… Test 0.0284 0.473 0.333 0.0200 0.296 0.181 0.173 0.324
Importamos la data
library(readxl)
library(forecast)
HCB <- read_excel("~/Econometria/tableExportHCB.xlsx",col_type= c("skip", "numeric"),
skip = 7)
## Serie
serie.HCB.ts<- ts(data = HCB, start = c(2017,1),
frequency = 12)
Yt_4<-serie.HCB.ts
serie.HCB.ts %>%
autoplot(main = "Importación de Hidrocarburos SLV 2017-2022 (mayo)",
xlab = "años/ trimestres",
ylab = "Indice")
library(TSstudio)
library(forecast)
ts_plot(Yt_4,Xtitle= "Años/Meses")
Orden de Integración
library(kableExtra)
library(magrittr)
d_4<-ndiffs(Yt_4)
D_4<-nsdiffs(Yt_4)
ordenes_integracion<-c(d_4,D_4)
names(ordenes_integracion)<-c("d","D")
ordenes_integracion %>% kable(caption = "Ordenes de Integracion") %>% kable_material()
| x | |
|---|---|
| d | 0 |
| D | 0 |
#Gráfico de la serie diferenciada
Yt_4 %>%
diff(lag=12,diffences=D_4) %>%
diff(diffences=d_4) %>%
ts_plot(title= "Yt estacionaria")
Yt_4 %>%
diff(lag=12,diffences=D_4) %>%
diff(diffences=d_4) %>%
ts_cor(lag.max = 36)
Usando forecast
library(forecast)
library(ggthemes)
modelo_estimado_HCB<-Yt_4 %>%
Arima(order = c(1, 0, 1),
seasonal = c(1, 0, 0))
summary(modelo_estimado_HCB)
## Series: .
## ARIMA(1,0,1)(1,0,0)[12] with non-zero mean
##
## Coefficients:
## ar1 ma1 sar1 mean
## 0.6652 -0.4851 -0.2022 215790.730
## s.e. 0.2952 0.3398 0.1519 5147.504
##
## sigma^2 = 1.042e+09: log likelihood = -741.69
## AIC=1493.39 AICc=1494.44 BIC=1504.1
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -428.6553 31240.12 25005.59 -2.421527 12.06669 0.6217076
## ACF1
## Training set -0.007961808
#grafico de las raices (estabilidad/convertibilidad)
modelo_estimado_HCB %>% autoplot(type="both")+theme_solarized()
#verificacion
modelo_estimado_HCB %>% check_res(lag.max = 36)
## Warning: Ignoring 36 observations
## Warning: Ignoring 34 observations
## Warning: Ignoring 4 observations
library(forecast)
#Estimar el modelo
modeloHW_M_HC<-HoltWinters(x = Yt_4,
seasonal = "multiplicative",
optim.start = c(0.9,0.9,0.9))
modeloHW_M_HC
## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
##
## Call:
## HoltWinters(x = Yt_4, seasonal = "multiplicative", optim.start = c(0.9, 0.9, 0.9))
##
## Smoothing parameters:
## alpha: 0.1580326
## beta : 0
## gamma: 0.4946113
##
## Coefficients:
## [,1]
## a 2.501763e+05
## b 1.683628e+03
## s1 9.111372e-01
## s2 8.894296e-01
## s3 8.823914e-01
## s4 8.896954e-01
## s5 9.398533e-01
## s6 1.039846e+00
## s7 9.375516e-01
## s8 9.546773e-01
## s9 9.615739e-01
## s10 1.097096e+00
## s11 9.326123e-01
## s12 8.393040e-01
#Generar el pronostico:
pronosticosMW_M_HC<-forecast(object = modeloHW_M_HC,h=12, level=c(0.95, 0.90))
pronosticosMW_M_HC
## Point Forecast Lo 95 Hi 95 Lo 90 Hi 90
## Apr 2022 229478.9 182755.6 276202.2 190267.5 268690.3
## May 2022 225509.1 177130.5 273887.7 184908.5 266109.7
## Jun 2022 225210.3 175204.2 275216.3 183243.9 267176.6
## Jul 2022 228572.3 176863.3 280281.4 185176.7 271967.9
## Aug 2022 243040.8 189049.4 297032.2 197729.8 288351.8
## Sep 2022 270649.0 213366.1 327932.0 222575.6 318722.4
## Oct 2022 245602.6 188944.1 302261.0 198053.3 293151.9
## Nov 2022 251696.2 193221.1 310171.3 202622.3 300770.0
## Dec 2022 255133.4 195092.9 315173.9 204745.8 305520.9
## Jan 2023 292938.5 227754.8 358122.3 238234.6 347642.4
## Feb 2023 250589.4 188983.5 312195.2 198888.1 302290.6
## Mar 2023 226930.8 188989.5 264872.2 195089.4 258772.2
#Gráfico comparativo
#Gráfico comparativo
Yt_Sarima_M_HC<-modelo_estimado_HCB$fitted
Yt_HW_M_HC<-pronosticosMW_M_HC$fitted
grafico_comparativo_4<-cbind(Yt_4, Yt_Sarima_M_HC,Yt_HW_M_HC)
ts_plot(grafico_comparativo_4)
library(tsibble)
library(feasts)
library(fable)
library(fabletools)
library(tidyr)
library(dplyr)
#generando modelo automatico
a_4<-Yt_4 %>% as_tsibble() %>%
model(arima_original=ARIMA(value ~ pdq(1, 0, 1) + PDQ(1, 0, 0)),
arima_010_011 = ARIMA(value ~ pdq(1, 1, 1) + PDQ(1, 0, 0)),
arima_010_110 = ARIMA(value ~ pdq(1, 0, 1) + PDQ(1, 1, 0)),
arima_automatico=ARIMA(value,ic="bic",stepwise = FALSE)
#modelo automatico
)
## 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_010_011
## [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
print(a_4)
## # A mable: 1 x 4
## arima_original arima_010_011 arima_010_110
## <model> <model> <model>
## 1 <ARIMA(1,0,1)(1,0,0)[12] w/ mean> <NULL model> <ARIMA(1,0,1)(1,1,0)[12]>
## # … with 1 more variable: arima_automatico <model>
#ordenar
a_4 %>% pivot_longer(everything(), names_to = "Model name",
values_to = "Orders") %>% glance() %>%
arrange(AICc) ->tabla_4
tabla_4
## # A tibble: 3 × 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_110 Orders 1.47e9 -612. 1232. 1233. 1240. <cpl> <cpl>
## 2 arima_automatico Orders 1.10e9 -745. 1493. 1494. 1498. <cpl> <cpl>
## 3 arima_original Orders 1.04e9 -742. 1493. 1494. 1504. <cpl> <cpl>
library(forecast)
library(dplyr)
library(tsibble)
library(fable)
library(fabletools)
#convertimos
Yt<-Yt_4 %>%
as_tsibble() %>%
rename(M_HC= value)
#generador de la dta (roja)
data.cross.validation.M_HC<-Yt %>%
as_tsibble() %>%
stretch_tsibble(.init = 60,.step = 1)
#generador de pronostico (verde) y simulacion
TSCV_4<-data.cross.validation.M_HC %>%
model(ARIMA(M_HC ~ 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_4)
## # A tibble: 1 × 10
## .model .type ME RMSE MAE MPE MAPE MASE RMSSE ACF1
## <chr> <chr> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
## 1 ARIMA(M_HC ~ pdq(0… Test -35010. 41442. 35010. -16.9 16.9 0.870 0.839 -0.194