PIB trimestral precios corrientes SLV 2009-2021

Carga de datos

library(readxl)
library(forecast)
datos_PIB <-
  read_excel("C:/Users/8abla/OneDrive/Documents/EMA118/modelo sarima/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)
## 
## 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(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)
library(forecast)

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

Identificación

Orden de integración

library(kableExtra)
## 
## 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)
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)
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("C:/Users/8abla/OneDrive/Documents/EMA118/modelo sarima/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]>