PIB trimestral precios corrientes SLV 2009-2021

Carga de datos

library(readxl)
library(forecast)
datos_PIB <-
  read_excel("PIB_trimestral.xlsx",
             col_types = c("skip", "numeric"))

datos_PIB_ts <- ts(data = datos_PIB,
                    start = c(2009, 1),
                    frequency = 4)

Componente de tendencia Tt

ma2_12_PIB <- ma(datos_PIB_ts, 4, centre = TRUE)
autoplot(datos_PIB_ts,main = "PIB, El Salvador 2009-2021",
           xlab = "Años/Meses",
           ylab = "Indice")+
  autolayer(ma2_12_PIB,series = "Tt")
## Warning: Removed 4 row(s) containing missing values (geom_path).

Factores estacionales

library(magrittr)
## 
## Attaching package: 'magrittr'
## The following object is masked from 'package:tidyr':
## 
##     extract
Yt_PIB <- datos_PIB_ts 
Tt_PIB <- ma2_12_PIB
SI_PIB <- Yt_PIB - Tt_PIB

St_PIB <- tapply(SI_PIB, cycle(SI_PIB), mean, na.rm = TRUE) 
St_PIB <- St_PIB - sum(St_PIB) / 4 

St_PIB <-
  rep(St_PIB, len = length(Yt_PIB)) %>% ts(start = c(2009, 1), frequency = 4) 

library(tsibble)
## Warning: package 'tsibble' was built under R version 4.2.1
## 
## 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.2.1
## Loading required package: fabletools
## Warning: package 'fabletools' was built under R version 4.2.1
## 
## Attaching package: 'fabletools'
## The following objects are masked from 'package:forecast':
## 
##     accuracy, forecast
library(ggplot2)
Yt_PIB %>% as_tsibble() %>%
  model(
    classical_decomposition(value, type = "additive")
  ) %>%
  components() %>%
  autoplot() +
  labs(title = "Descomposición Clásica Aditiva, PIB")+xlab("Años/Meses")
## Warning: Removed 2 row(s) containing missing values (geom_path).

Serie PIB-Metodología Box-Jenkins

library(TSstudio)
## Warning: package 'TSstudio' was built under R version 4.2.1
library(forecast)

ts_plot(Yt_PIB,Xtitle = "Años/Meses")

Identificación

Orden de integración

library(kableExtra)
## Warning: package 'kableExtra' was built under R version 4.2.1
## 
## Attaching package: 'kableExtra'
## The following object is masked from 'package:dplyr':
## 
##     group_rows
library(magrittr)
d_PIB<-ndiffs(Yt_PIB)
D_PIB<-nsdiffs(Yt_PIB)
ordenes_integracion<-c(d_PIB,D_PIB)
names(ordenes_integracion)<-c("d_PIB","D_PIB")
ordenes_integracion %>% kable(caption = "Ordenes de Integración") %>% kable_material()
Ordenes de Integración
x
d_PIB 1
D_PIB 1 
#Gráfico de la serie diferenciada
Yt_PIB %>% 
  diff(lag = 4,diffences=D_PIB) %>% 
  diff(diffences=d_PIB) %>% 
  ts_plot(title = "Yt estacionaria")

Verificar los valores para (p,q) & (P,Q)

##Se construyen las funciones de autocorrelación simple y parcial.

Yt_PIB %>% 
  diff(lag = 4,diffences=D_PIB) %>% 
  diff(diffences=d_PIB) %>% 
  ts_cor(lag.max = 12)

Estimación del modelo Sarima

SARIMA(0,1,0)(1,1,1)[4]

Usando Forecast

library(forecast)
library(ggthemes)
## Warning: package 'ggthemes' was built under R version 4.2.1
modelo_estimado_PIB <- Yt_PIB %>% 
  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 = 12)
## Warning: Ignoring 12 observations
## Warning: Ignoring 10 observations
## Warning: Ignoring 4 observations

Grafico comparativo

Yt_Sarima_PIB<-modelo_estimado_PIB$fitted
grafico_comparativo<-cbind(Yt_PIB,Yt_Sarima_PIB)
ts_plot(grafico_comparativo)

Verificación de sobre ajuste/sub ajuste

Se estimarán los modelos: Partiendo del modelo original SARIMA(0,1,0)(1,1,1)[4], se estima un nuevo modelo con P-1: SARIMA(0,1,0)(0,1,1)[4], y otro con Q-1: SARIMA(0,1,0)(1,1,0)[4]

library(tsibble)
library(feasts)
library(fable)
## Warning: package 'fable' was built under R version 4.2.1
library(fabletools)
library(tidyr)
library(dplyr)
a_PIB<-Yt_PIB %>% 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_PIB)
## # 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_PIB %>% pivot_longer(everything(), names_to = "Model name",
                         values_to = "Orders") %>% glance() %>% 
  arrange(AICc) ->tabla_PIB
tabla_PIB
## # 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]>

Validación cruzada

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

Yt_PIB<-Yt_PIB %>% as_tsibble() %>% rename(PIB=value)

data.cross.validation_PIB<-Yt_PIB %>% 
  as_tsibble() %>% 
  stretch_tsibble(.init = 20,.step = 1)

TSCV_PIB<-data.cross.validation_PIB %>% 
model(ARIMA(PIB ~ pdq(0, 1, 0) + PDQ(1, 1, 1))) %>% 
forecast(h=1) %>% accuracy(Yt_PIB)
## 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_PIB)
## # 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

Intepretación del MAPE: en la validación cruzada obtuvimos un MAPE de 3.21% teniendo una diferencia al anterior, que era de 2.05%. Podemos decir entonces que nuestro modelo sigue teniendo un gran poder predictivo.

IPC general SLV 2009-2022 (mayo)

Carga de datos

library(readxl)
library(forecast)
datos_IPC <-
  read_excel("datos_IPC.xlsx",
             col_types = c("skip", "numeric"))

datos_IPC_ts <- ts(data = datos_IPC,
                    start = c(2009, 1),
                    frequency = 12)

Componente de tendencia Tt

ma2_12_IPC <- ma(datos_IPC_ts, 12, centre = TRUE)
autoplot(datos_IPC_ts,main = "IPC, El Salvador 2009-2022",
           xlab = "Años/Meses",
           ylab = "Indice")+
  autolayer(ma2_12_IPC,series = "Tt")
## Warning: Removed 12 row(s) containing missing values (geom_path).

Factores estacionales

library(magrittr)
Yt_IPC <- datos_IPC_ts 
Tt_IPC <- ma2_12_IPC
SI_IPC <- Yt_IPC - Tt_IPC

St_IPC <- tapply(SI_IPC, cycle(SI_IPC), mean, na.rm = TRUE) 
St_IPC <- St_IPC - sum(St_IPC) / 12 

St_IPC <-
  rep(St_IPC, len = length(Yt_IPC)) %>% ts(start = c(2009, 1), frequency = 12) 

library(tsibble)
library(feasts)
library(ggplot2)
Yt_IPC %>% as_tsibble() %>%
  model(
    classical_decomposition(value, type = "additive")
  ) %>%
  components() %>%
  autoplot() +
  labs(title = "Descomposición Clásica Aditiva, IPC")+xlab("Años/Meses")
## Warning: Removed 6 row(s) containing missing values (geom_path).

Serie PIB-Metodología Box-Jenkins

library(TSstudio)
library(forecast)

ts_plot(Yt_IPC,Xtitle = "Años/Meses")

Identificación

Orden de integración

library(kableExtra)
library(magrittr)
d_IPC<-ndiffs(Yt_IPC)
D_IPC<-nsdiffs(Yt_IPC)
ordenes_integracion<-c(d_IPC,D_IPC)
names(ordenes_integracion)<-c("d_IPC","D_IPC")
ordenes_integracion %>% kable(caption = "Ordenes de Integración") %>% kable_material()
Ordenes de Integración
x
d_IPC 1
D_IPC 0 
#Gráfico de la serie diferenciada
Yt_IPC %>% 
  diff(lag = 12,diffences=D_IPC) %>% 
  diff(diffences=d_IPC) %>% 
  ts_plot(title = "Yt estacionaria")

Verificar los valores para (p,q) & (P,Q)

##Se construyen las funciones de autocorrelación simple y parcial.

Yt_IPC %>% 
  diff(lag = 12,diffences=D_IPC) %>% 
  diff(diffences=d_IPC) %>% 
  ts_cor(lag.max = 36)

Estimación del modelo Sarima

SARIMA(1,1,1)(1,0,3)[12]

Usando Forecast

library(forecast)
library(ggthemes)
modelo_estimado_IPC <- Yt_IPC %>% 
  Arima(order = c(1, 1, 1),
        seasonal = c(1, 1, 3))

summary(modelo_estimado_IPC)
## Series: . 
## ARIMA(1,1,1)(1,1,3)[12] 
## 
## Coefficients:
##          ar1      ma1    sar1     sma1    sma2    sma3
##       0.8881  -0.7426  0.3184  -1.2531  0.1176  0.1360
## s.e.  0.1928   0.2594  0.8327   0.9348  0.8045  0.1308
## 
## sigma^2 = 0.2114:  log likelihood = -108.14
## AIC=230.29   AICc=231.09   BIC=251.27
## 
## Training set error measures:
##                      ME     RMSE       MAE       MPE      MAPE     MASE
## Training set 0.01597872 0.431751 0.2803503 0.0130507 0.2553127 0.151921
##                    ACF1
## Training set 0.06957202
modelo_estimado_IPC %>% autoplot(type="both")+theme_solarized()

modelo_estimado_IPC %>% check_res(lag.max = 36)
## Warning: Ignoring 36 observations
## Warning: Ignoring 34 observations
## Warning: Ignoring 4 observations

Grafico comparativo

Yt_Sarima_IPC<-modelo_estimado_IPC$fitted
grafico_comparativo<-cbind(Yt_IPC,Yt_Sarima_IPC)
ts_plot(grafico_comparativo)

Verificación de sobre ajuste/sub ajuste

Se estimarán los modelos: Partiendo del modelo original SARIMA(1,1,1)(1,0,3)[12], se estima un nuevo modelo con P-1: SARIMA(1,1,1)(0,0,3)[12], y otro con Q-1: SARIMA(1,1,1)(1,0,2)[12]

library(tsibble)
library(feasts)
library(fable)
library(fabletools)
library(tidyr)
library(dplyr)
a_IPC<-Yt_IPC %>% as_tsibble() %>% 
  model(arima_original=ARIMA(value ~ pdq(1, 1, 1) + PDQ(1, 0, 3)),
        arima_010_011 = ARIMA(value ~ pdq(1, 1, 1) + PDQ(0, 0, 3)),
        arima_010_110 = ARIMA(value ~ pdq(1, 1, 1) + PDQ(1, 0, 2)),
        arima_automatico=ARIMA(value,ic="bic",stepwise = FALSE)
  )
print(a_IPC)
## # A mable: 1 x 4
##                       arima_original                      arima_010_011
##                              <model>                            <model>
## 1 <ARIMA(1,1,1)(1,0,3)[12] w/ drift> <ARIMA(1,1,1)(0,0,3)[12] w/ drift>
## # … with 2 more variables: arima_010_110 <model>, arima_automatico <model>
a_IPC %>% pivot_longer(everything(), names_to = "Model name",
                         values_to = "Orders") %>% glance() %>% 
  arrange(AICc) ->tabla_IPC
tabla_IPC
## # 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.219   -103.  220.  221.  241. <cpl [1]>  <cpl [37]>
## 2 arima_010_110    Orders  0.219   -103.  220.  221.  241. <cpl [13]> <cpl [25]>
## 3 arima_original   Orders  0.221   -103.  222.  223.  246. <cpl [13]> <cpl [37]>
## 4 arima_automatico Orders  0.225   -105.  223.  224.  241. <cpl [27]> <cpl [0]>
library(forecast)
library(dplyr)
library(tsibble)
library(fable)
library(fabletools)

Yt_IPC<-Yt_IPC %>% as_tsibble() %>% rename(IPC=value)

data.cross.validation_IPC<-Yt_IPC %>% 
  as_tsibble() %>% 
  stretch_tsibble(.init = 60,.step = 1)

TSCV_IPC<-data.cross.validation_IPC %>% 
model(ARIMA(IPC ~ pdq(1, 1, 1) + PDQ(1, 0, 3))) %>% 
forecast(h=1) %>% accuracy(Yt_IPC)
## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf

## Warning in max(PDQ$Q): ningun argumento finito para max; retornando -Inf
## Warning in sqrt(diag(best$var.coef)): Se han producido NaNs

## Warning in sqrt(diag(best$var.coef)): Se han producido NaNs
## 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 in sqrt(diag(best$var.coef)): Se han producido NaNs
## 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)`.
## 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)`.
## 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: 52 errors (1 unique) encountered for ARIMA(IPC ~ pdq(1, 1, 1) + PDQ(1, 0, 3))
## [52] 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. 
## 1 observation is missing at 2022 jun.
print(TSCV_IPC)
## # 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(1, 1,… Test  0.0411 0.349 0.265 0.0330 0.232 0.143 0.128 0.398

Intepretación del MAPE: en la validación cruzada obtuvimos un MAPE de 0.23% teniendo una diferencia al anterior, que era de 0.25%. Podemos decir entonces que nuestro modelo sigue teniendo un gran poder predictivo.

4. Importación de Hidrocarburos SLV 2017-2022 (mayo)

Carga de datos

library(readxl)
library(forecast)
datos_comercio <-
  read_excel("datos_comercio.xlsx",
             col_types = c("skip", "numeric"))

datos_comercio_ts <- ts(data = datos_comercio,
                    start = c(2017, 1),
                    frequency = 12)

Componente de tendencia Tt

ma2_12_comercio <- ma(datos_comercio_ts, 12, centre = TRUE)
autoplot(datos_comercio_ts,main = "Importación de Hidrocarburos, El Salvador 2017-2022",
           xlab = "Años/Meses",
           ylab = "Indice")+
  autolayer(ma2_12_comercio,series = "Tt")
## Warning: Removed 12 row(s) containing missing values (geom_path).

Factores estacionales

library(magrittr)
Yt_comercio <- datos_comercio_ts 
Tt_comercio <- ma2_12_comercio
SI_comercio <- Yt_comercio - Tt_comercio

St_comercio <- tapply(SI_comercio, cycle(SI_comercio), mean, na.rm = TRUE) 
St_comercio <- St_comercio - sum(St_comercio) / 12 

St_comercio <-
  rep(St_comercio, len = length(Yt_comercio)) %>% ts(start = c(2017, 1), frequency = 12) 

library(tsibble)
library(feasts)
library(ggplot2)
Yt_comercio %>% as_tsibble() %>%
  model(
    classical_decomposition(value, type = "additive")
  ) %>%
  components() %>%
  autoplot() +
  labs(title = "Descomposición Clásica Aditiva de la Importación de Hidrocarburos")+xlab("Años/Meses")
## Warning: Removed 6 row(s) containing missing values (geom_path).

Serie PIB-Metodología Box-Jenkins

library(TSstudio)
library(forecast)

ts_plot(Yt_comercio,Xtitle = "Años/Meses")

Identificación

Orden de integración

library(kableExtra)
library(magrittr)
d_comercio<-ndiffs(Yt_comercio)
D_comercio<-nsdiffs(Yt_comercio)
ordenes_integracion<-c(d_comercio,D_comercio)
names(ordenes_integracion)<-c("d_comercio","D_comercio")
ordenes_integracion %>% kable(caption = "Ordenes de Integración") %>% kable_material()
Ordenes de Integración
x
d_comercio 0
D_comercio 0 
#Gráfico de la serie diferenciada
Yt_comercio %>% 
  diff(lag = 12,diffences=D_comercio) %>% 
  diff(diffences=d_comercio) %>% 
  ts_plot(title = "Yt estacionaria")

Verificar los valores para (p,q) & (P,Q)

##Se construyen las funciones de autocorrelación simple y parcial.

Yt_comercio %>% 
  diff(lag = 12,diffences=D_comercio) %>% 
  diff(diffences=d_comercio) %>% 
  ts_cor(lag.max = 36)

Estimación del modelo Sarima

SARIMA(1,0,1)(1,0,1)[12]

Usando Forecast

library(forecast)
library(ggthemes)
modelo_estimado_comercio <- Yt_comercio %>% 
  Arima(order = c(1, 0, 1),
        seasonal = c(1, 0, 1))

summary(modelo_estimado_comercio)
## Series: . 
## ARIMA(1,0,1)(1,0,1)[12] with non-zero mean 
## 
## Coefficients:
##          ar1      ma1     sar1    sma1        mean
##       0.7511  -0.5712  -0.9386  0.8262  215475.626
## s.e.  0.1966   0.2292   0.3393  0.5689    6060.192
## 
## sigma^2 = 1.01e+09:  log likelihood = -765.03
## AIC=1542.06   AICc=1543.51   BIC=1555.11
## 
## Training set error measures:
##                    ME     RMSE     MAE       MPE     MAPE      MASE        ACF1
## Training set 75.04379 30528.35 24731.4 -2.144265 11.95625 0.6084344 -0.04328773
modelo_estimado_comercio %>% autoplot(type="both")+theme_solarized()

modelo_estimado_comercio %>% check_res(lag.max = 36)
## Warning: Ignoring 36 observations
## Warning: Ignoring 34 observations
## Warning: Ignoring 4 observations

Grafico comparativo

Yt_Sarima_comercio<-modelo_estimado_comercio$fitted
grafico_comparativo<-cbind(Yt_comercio,Yt_Sarima_comercio)
ts_plot(grafico_comparativo)

Verificación de sobre ajuste/sub ajuste

Se estimarán los modelos: Partiendo del modelo original SARIMA(1,0,1)(1,0,1)[12], se estima un nuevo modelo con P-1: SARIMA(1,0,1)(0,0,1)[12], y otro con Q-1: SARIMA(1,0,1)(1,0,0)[12]

library(tsibble)
library(feasts)
library(fable)
library(fabletools)
library(tidyr)
library(dplyr)
a_comercio<-Yt_comercio %>% as_tsibble() %>% 
  model(arima_original=ARIMA(value ~ pdq(1, 0 ,1) + PDQ(1, 0, 1)),
        arima_010_011 = ARIMA(value ~ pdq(1, 0, 1) + PDQ(0, 0, 1)),
        arima_010_110 = ARIMA(value ~ pdq(1, 0, 1) + PDQ(1, 0, 0)),
        arima_automatico=ARIMA(value,ic="bic",stepwise = FALSE)
  )
## 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_original
## [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_comercio)
## # A mable: 1 x 4
##   arima_original                     arima_010_011
##          <model>                           <model>
## 1   <NULL model> <ARIMA(1,0,1)(0,0,1)[12] w/ mean>
## # … with 2 more variables: arima_010_110 <model>, arima_automatico <model>
a_comercio %>% pivot_longer(everything(), names_to = "Model name",
                         values_to = "Orders") %>% glance() %>% 
  arrange(AICc) ->tabla_comercio
tabla_comercio
## # 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_automatico Orders     1.11e9   -769. 1541. 1541. 1545. <cpl>    <cpl>   
## 2 arima_010_110    Orders     1.06e9   -766. 1542. 1543. 1552. <cpl>    <cpl>   
## 3 arima_010_011    Orders     1.10e9   -767. 1544. 1545. 1555. <cpl>    <cpl>
library(forecast)
library(dplyr)
library(tsibble)
library(fable)
library(fabletools)

Yt_comercio<-Yt_comercio %>% as_tsibble() %>% rename(comercio=value)

data.cross.validation_comercio<-Yt_comercio %>% 
  as_tsibble() %>% 
  stretch_tsibble(.init = 60,.step = 1)

TSCV_comercio<-data.cross.validation_comercio %>% 
model(ARIMA(comercio ~ pdq(1, 0, 1) + PDQ(1, 0, 1))) %>% 
forecast(h=1) %>% accuracy(Yt_comercio)
## 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)`.
## 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)`.
## 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)`.
## 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: 4 errors (1 unique) encountered for ARIMA(comercio ~ pdq(1, 0, 1) + PDQ(1, 0, 1))
## [4] 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. 
## 1 observation is missing at 2022 jun.
print(TSCV_comercio)
## # 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(comercio ~ pdq… Test  23348. 32825. 23348.  8.95  8.95 0.574 0.662  -0.5

Intepretación del MAPE: en la validación cruzada obtuvimos un MAPE de 8.94% teniendo una diferencia al anterior, que era de 11.95%. Podemos decir entonces que nuestro modelo pierde poder predictivo.