Carga de datos

IVAE_03_22 <-
  read_excel(
    "C:/Users/Keiry/Documents/Eco22/IVAE_22_03.xlsx",
    col_names = FALSE,
    skip = 6,
    n_max = 10
  )
IVAE_03_22[1:10, 1:10]

## # A tibble: 10 x 10
##    ...1                     ...2  ...3  ...4  ...5  ...6  ...7  ...8  ...9 ...10
##    <chr>                   <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl> <dbl>
##  1   1   IVAE               86.7  80.8  87.2  83.9  91.4  93.5  86.4  86.7  87.6
##  2   2   Agricultura, Gan~  86.8  70.2  65.2  67.4 138.  153.   75.6 137.  110. 
##  3   3   Índice de Produc~  86.6  85.4 100.   91.6  91.2  89.7  86.8  81.2  86.6
##  4   4   Construcción       66.6  76.9  82.6  82.6  69.0  80.8  82.5  67.3  74.5
##  5   5   Comercio, Transp~  91.0  74.8  77.8  80.0  89.6  90.6  81.9  82.3  80.1
##  6   6   Información y Co~  94.6  76.9  81.9  81.0  92.3  91.9  96.9  81.3 104. 
##  7   7   Actividades Fina~  97.4  85.7  94.9  88.9  92.6  91.6  92.1  94.2  92.6
##  8   8   Actividades Inmo~  96.8  95.2  94.6  93.9  94.0  94.6  95.3  95.9  96.4
##  9   9   Actividades Prof~  80.7  76.9  81.0  77.2  86.1  87.8  82.9  75.9  80.4
## 10   10   Actividades de ~  80.8  83.0  91.3  83.9  84.9  86.7  90.3  87.7  88.4

#descomposicion de serie
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)

Pronóstico de series temporales, enfoque moderno (estocástico)

ts.plot(data.ivae.ts, Xtitle = "Años/Meses")

Identificación

d <- ndiffs(data.ivae.ts)
D <- nsdiffs(data.ivae.ts)
ordenes_integracion <- c(d, D)
names(ordenes_integracion) <- c("d", "D")
ordenes_integracion %>% kable(caption = "Ordenes de Integración") %>% kable_material()

#grafico de la serie diferenciada
data.ivae.ts %>%
  diff(lag = 12, diffences = D) %>%
  diff(diffences = d) %>%
  ts_plot(title = "Yt estacionaria")

data.ivae.ts %>%
  diff(lag = 12, diffences = D) %>%
  diff(diffences = d) %>%
  ts_cor(lag.max = 36)

modelo_estimado <-
  data.ivae.ts %>% 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)

#Grafico
data.ivae.ts_Sarima <- modelo_estimado$fitted
#cargando Pronóstico HW
#Generar el pronostico
PronosticoHW2 <-
  hw(
    y = data.ivae.ts,
    h = 12,
    level = c(0.96, 0.99),
    seasonal = "multiplicative",
    initial = "optimal"
  )
data.ivae.ts_HW <- PronosticoHW2$fitted
grafico_comparativo <-
  cbind(data.ivae.ts, data.ivae.ts_Sarima, data.ivae.ts_HW)
ts_plot(grafico_comparativo)

Verificación de sobre ajuste/sub ajuste

a <-
  data.ivae.ts %>% 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]>

Validación Cruzada (Cross Validated)

Yt <- data.ivae.ts %>% 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)

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

Carga de datos

serie.ipc <-
  read_excel("C:/Users/Keiry/Documents/Eco22/IPC_2.xlsx",
             col_types = c("skip", "numeric"),
             skip = 5)

Descomposición de Serie Temporal

Serie temporal

data.ipc.ts <- serie.ipc %>% ts(start = c(2009, 1), frequency = 12)

Modelo aditivo

ma2_12 <- ma(data.ipc.ts, 12, centre = TRUE)
autoplot(data.ipc.ts,main = "IPC, El Salvador 2009-2022[mayo]",
           xlab = "Años/Meses",
           ylab = "Indice")+
  autolayer(ma2_12,series = "Tt")

Calculo de los factores estacionales

Yt <- data.ipc.ts #Serie original
Tt <- ma2_12 #Media móvil centrada (2x12-MA) como componente de Tendencia Ciclo
SI <- Yt - Tt #Diferencia que contiene componentes Estacional e Irregular

St <- tapply(SI, cycle(SI), mean, na.rm = TRUE) #Promediando los resultados de cada mes
#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) 
autoplot(St,
         main = "Factores Estacionales",
         xlab = "Años/Meses",
         ylab = "Factor Estacional") 

Componente irregular

It<-Yt-Tt-St
autoplot(It,
         main = "Componente Irregular",
         xlab = "Aos/Meses",
         ylab = "It")

## Descomposicion usando libreria TStudio

ts_decompose(data.ipc.ts,type = "additive",showline = TRUE)

ts_seasonal(data.ipc.ts,type = "box",title = "Análisis de Valores Estacionales")

Pronóstico Enfoque clásico

#Estimar el modelo
ModeloHW <- HoltWinters(x = data.ipc.ts,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 = data.ipc.ts, 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

PronosticosHW<-forecast(object = ModeloHW, h = 12, level = c(0.90,0.95))
PronosticosHW

##          Point Forecast    Lo 90    Hi 90    Lo 95    Hi 95
## Jun 2022       123.7827 122.8208 124.7446 122.6365 124.9289
## Jul 2022       123.9896 122.6975 125.2818 122.4499 125.5294
## Aug 2022       124.0357 122.4498 125.6217 122.1459 125.9256
## Sep 2022       124.3884 122.5223 126.2546 122.1648 126.6120
## Oct 2022       125.0394 122.8974 127.1814 122.4870 127.5918
## Nov 2022       125.4570 123.0461 127.8680 122.5842 128.3298
## Dec 2022       125.6139 122.9402 128.2875 122.4280 128.7997
## Jan 2023       126.3012 123.3543 129.2482 122.7897 129.8128
## Feb 2023       127.1819 123.9551 130.4087 123.3369 131.0269
## Mar 2023       128.1007 124.5897 131.6116 123.9171 132.2842
## Apr 2023       128.6771 124.8876 132.4667 124.1616 133.1927
## May 2023       129.2228 124.7841 133.6615 123.9338 134.5118

#Grafico de la serie original y del pronostico
PronosticosHW %>% autoplot()

Pronostico enfoque moderno

ts_plot(data.ipc.ts,Xtitle = "Años/Meses")

Orden de integracion

d<-ndiffs(data.ipc.ts)
D<-nsdiffs(data.ipc.ts)
ordenes_integracion<-c(d,D)
names(ordenes_integracion)<-c("d","D")
ordenes_integracion %>% kable(caption = "Ordenes de Integración") %>% kable_material()

#Grafico de serie diferenciada
data.ipc.ts %>% 
  diff(lag = 12,diffences=D) %>% 
  diff(diffences=d) %>% 
  ts_plot(title = "Serie estacionaria IPC")

Verificando los valores para (p,q) (P,Q)

data.ipc.ts %>% 
  diff(lag = 12,diffences=D) %>% 
  diff(diffences=d) %>% 
  ts_cor(lag.max = 36)

Estimacion del modelo

modelo_estimado <- data.ipc.ts %>% 
  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.1112  -1.0000
## s.e.  0.0950   0.1139
## 
## sigma^2 = 0.2312:  log likelihood = -114.66
## AIC=235.33   AICc=235.49   BIC=244.32
## 
## Training set error measures:
##                      ME      RMSE       MAE        MPE      MAPE      MASE
## Training set 0.02246671 0.4579215 0.3010507 0.01772013 0.2729747 0.1631385
##                   ACF1
## Training set 0.2438271

modelo_estimado %>% autoplot(type="both")+theme_solarized()

modelo_estimado %>% check_res(lag.max = 36)

data_ipc_Sarima<-modelo_estimado$fitted
data_ipc_HW<-PronosticosHW$fitted
grafico_comparativo<-cbind(data.ipc.ts,data_ipc_Sarima,data_ipc_HW)
ts_plot(grafico_comparativo)

Verificacion de ajuste

a_ipc<-data.ipc.ts %>% 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_ipc)

## # 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_ipc %>% pivot_longer(everything(), names_to = "Model name",
                         values_to = "Orders") %>% glance() %>% 
  arrange(AICc) ->tabla_1
tabla_1

## # 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  0.225   -105.  223.  224.  241. <cpl [27]> <cpl [0]> 
## 2 arima_010_011    Orders  0.228   -115.  235.  235.  241. <cpl [0]>  <cpl [12]>
## 3 arima_original   Orders  0.231   -115.  235.  235.  244. <cpl [12]> <cpl [12]>
## 4 arima_010_110    Orders  0.348   -132.  268.  268.  274. <cpl [12]> <cpl [0]>

Validacion Cruzada

Yt_ipc<-data.ipc.ts %>% as_tsibble() %>% rename(IPC=value)

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

TSCV_ipc<-data.cross.validation %>% 
model(ARIMA(IPC ~ pdq(0, 1, 0) + PDQ(1, 1, 1))) %>% 
forecast(h=1) %>% accuracy(Yt_ipc)

print(TSCV_ipc)

## # 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,~ Test  0.0284 0.473 0.333 0.0200 0.296 0.181 0.173 0.324

Carga de datos

serie.pib <-
  read_excel("C:/Users/Keiry/Documents/Eco22/PIB_2.xlsx",
             col_types = c("skip", "numeric"),
             skip = 5)

Descomposición de Serie Temporal

Serie temporal

data.pib.ts <- serie.pib %>% ts(start = c(2009, 1), frequency = 4)

Modelo aditivo

ma2_12_1 <- ma(data.pib.ts, 4, centre = TRUE)

Calculo de los factores estacionales

Yt_pib <- data.pib.ts #Serie original
Tt_pib <- ma2_12_1 #Media móvil centrada (2x12-MA) como componente de Tendencia Ciclo
SI_1 <- Yt_pib - Tt_pib #Diferencia que contiene componentes Estacional e Irregular

St_pib <- tapply(SI_1, cycle(SI_1), mean, na.rm = TRUE) #Promediando los resultados de cada mes
#Los factores estacionales deben sumar "0" en el modelo aditivo
St_pib <- St_pib - sum(St_pib) / 4 
#Generar la serie de factores para cada valor de la serie original
St_pib <-
  rep(St_pib, len = length(Yt_pib)) %>% ts(start = c(2009, 1), frequency = 4) 
autoplot(St_pib,
         main = "Factores Estacionales",
         xlab = "Años/Meses",
         ylab = "Factor Estacional") 

Componente irregular

It_pib<-Yt_pib-Tt_pib-St_pib
autoplot(It_pib,
         main = "Componente Irregular",
         xlab = "Aos/Meses",
         ylab = "It_pib")

Descomposicion aditiva

descomposicion_aditiva_pib <- decompose(data.pib.ts, type = "additive")
autoplot(descomposicion_aditiva_pib,main="Descomposición Aditiva",xlab="Años/Meses")

Descomposicion usando libreria TStudio

ts_decompose(data.pib.ts,type = "additive",showline = TRUE)

ts_seasonal(data.pib.ts,type = "box",title = "Análisis de Valores Estacionales")

Pronóstico Enfoque clásico

#Estimar el modelo
ModeloHW_pib <- HoltWinters(x = data.pib.ts,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 = data.pib.ts, 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

PronosticosHW_pib<-forecast(object = ModeloHW_pib, h = 4, level = c(0.90,0.95))
PronosticosHW_pib

##         Point Forecast    Lo 90    Hi 90    Lo 95    Hi 95
## 2022 Q1       7292.485 6918.969 7666.001 6847.413 7737.557
## 2022 Q2       7570.373 7033.775 8106.971 6930.977 8209.768
## 2022 Q3       7411.019 6769.192 8052.847 6646.235 8175.804
## 2022 Q4       7933.429 7250.673 8616.184 7119.875 8746.982

#Grafico de la serie original y del pronostico
PronosticosHW_pib %>% autoplot()

Pronostico enfoque moderno

ts_plot(data.pib.ts,Xtitle = "Años/Meses")

Orden de integracion

d_pib<-ndiffs(data.pib.ts)
D_pib<-nsdiffs(data.pib.ts)
ordenes_integracion_pib<-c(d_pib,D_pib)
names(ordenes_integracion_pib)<-c("d_","D")
ordenes_integracion_pib %>% kable(caption = "Ordenes de Integración") %>% kable_material()

#Grafico de serie diferenciada
data.pib.ts %>% 
  diff(lag = 4,diffences=D_pib) %>% 
  diff(diffences=d_pib) %>% 
  ts_plot(title = "Serie estacionaria PIB")

Verificando los valores para (p,q) (P,Q)

data.pib.ts %>% 
  diff(lag = 4,diffences=D_pib) %>% 
  diff(diffences=d_pib) %>% 
  ts_cor(lag.max = 36)

Estimacion del modelo

modelo_estimado_pib <- data.pib.ts %>% 
  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 = 36)

data_pib_Sarima<-modelo_estimado_pib$fitted
data_pib_HW<-PronosticosHW_pib$fitted
grafico_comparativo_pib<-cbind(data.pib.ts,data_pib_Sarima,data_pib_HW)
ts_plot(grafico_comparativo_pib)

Verificacion de ajuste

a_pib<-data.pib.ts %>% 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_2
tabla_2

## # 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_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]>

Validacion Cruzada

Yt_pib<-data.pib.ts %>% as_tsibble() %>% rename(PIB=value)

data.cross.validation<-Yt_pib %>% 
  as_tsibble() %>% 
  stretch_tsibble(.init = 30,.step = 1)

TSCV_pib<-data.cross.validation %>% 
model(ARIMA(PIB ~ pdq(0, 1, 0) + PDQ(1, 1, 2))) %>% 
forecast(h=4) %>% accuracy(Yt_pib)

print(TSCV_pib)

## # 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

Carga de datos

serie.hidro <-
  read_excel("C:/Users/Keiry/Documents/Eco22/Hidrocarburos_2.xlsx",
             col_types = c("skip", "numeric"),
             skip = 5)

Descomposición de Serie Temporal

Serie temporal

data.hidro.ts <- serie.hidro %>% ts(start = c(2017, 1), frequency = 12)

Modelo aditivo

ma2_12_2 <- ma(data.hidro.ts, 12, centre = TRUE)

Calculo de los factores estacionales

Yt_hidro <- data.hidro.ts #Serie original
Tt_hidro <- ma2_12_2 #Media móvil centrada (2x12-MA) como componente de Tendencia Ciclo
SI_2 <- Yt_hidro - Tt_hidro #Diferencia que contiene componentes Estacional e Irregular

St_hidro <- tapply(SI_2, cycle(SI_2), mean, na.rm = TRUE) #Promediando los resultados de cada mes
#Los factores estacionales deben sumar "0" en el modelo aditivo
St_hidro <- St_hidro - sum(St_hidro) / 12 
#Generar la serie de factores para cada valor de la serie original
St_hidro <-
  rep(St_hidro, len = length(Yt_hidro)) %>% ts(start = c(2017, 1), frequency = 12) 
autoplot(St_hidro,
         main = "Factores Estacionales",
         xlab = "Años/Meses",
         ylab = "Factor Estacional") 

Componente irregular

It_hidro<-Yt_hidro-Tt_hidro-St_hidro
autoplot(It_hidro,
         main = "Componente Irregular",
         xlab = "Aos/Meses",
         ylab = "It_hidro")

Descomposicion aditiva

descomposicion_aditiva_hidro <- decompose(data.hidro.ts, type = "additive")
autoplot(descomposicion_aditiva_hidro,main="Descomposición Aditiva",xlab="Años/Meses")

Descomposicion usando libreria TStudio

ts_decompose(data.hidro.ts,type = "additive",showline = TRUE)

ts_seasonal(data.hidro.ts,type = "box",title = "Análisis de Valores Estacionales")

Pronóstico Enfoque clásico

#Estimar el modelo
ModeloHW_hidro <- HoltWinters(x = data.hidro.ts,seasonal = "multiplicative",optim.start = c(0.9,0.9,0.9))
ModeloHW_hidro

## Holt-Winters exponential smoothing with trend and multiplicative seasonal component.
## 
## Call:
## HoltWinters(x = data.hidro.ts, seasonal = "multiplicative", optim.start = c(0.9,     0.9, 0.9))
## 
## Smoothing parameters:
##  alpha: 0.1356726
##  beta : 0
##  gamma: 0.4908596
## 
## Coefficients:
##             [,1]
## a   2.550794e+05
## b   1.239453e+03
## s1  9.784035e-01
## s2  9.249193e-01
## s3  8.986508e-01
## s4  8.789500e-01
## s5  8.949436e-01
## s6  9.477674e-01
## s7  1.050046e+00
## s8  9.488853e-01
## s9  9.671890e-01
## s10 9.777809e-01
## s11 1.117223e+00
## s12 9.485809e-01

PronosticosHW_hidro<-forecast(object = ModeloHW_hidro, h = 12, level = c(0.90,0.95))
PronosticosHW_hidro

##          Point Forecast    Lo 90    Hi 90    Lo 95    Hi 95
## May 2022       250783.2 213969.3 287597.1 206916.8 294649.7
## Jun 2022       238220.6 200453.3 275987.9 193218.1 283223.1
## Jul 2022       232568.8 193871.3 271266.3 186457.9 278679.7
## Aug 2022       228559.7 188968.7 268150.7 181384.1 275735.2
## Sep 2022       233827.9 193112.7 274543.1 185312.7 282343.0
## Oct 2022       248804.2 206570.1 291038.3 198479.1 299129.3
## Nov 2022       276955.4 232471.1 321439.7 223949.1 329961.7
## Dec 2022       251449.9 207524.2 295375.5 199109.2 303790.5
## Jan 2023       257499.0 212393.7 302604.4 203752.7 311245.4
## Feb 2023       261530.9 215370.6 307691.2 206527.5 316534.2
## Mar 2023       300212.8 250497.5 349928.1 240973.3 359452.3
## Apr 2023       256072.1 226710.1 285434.1 221085.1 291059.1

#Grafico de la serie original y del pronostico
PronosticosHW_hidro %>% autoplot()

Pronostico enfoque moderno

ts_plot(data.hidro.ts,Xtitle = "Años/Meses")

Orden de integracion

d_hidro<-ndiffs(data.hidro.ts)
D_hidro<-nsdiffs(data.hidro.ts)
ordenes_integracion_hidro<-c(d_hidro,D_hidro)
names(ordenes_integracion_hidro)<-c("d","D")
ordenes_integracion_hidro %>% kable(caption = "Ordenes de Integración") %>% kable_material()

#Grafico de serie diferenciada
data.hidro.ts %>% 
  diff(lag = 12,diffences=D_hidro) %>% 
  diff(diffences=d_hidro) %>% 
  ts_plot(title = "Serie estacionaria Hidrocarburos")

Verificando los valores para (p,q) (P,Q)

data.hidro.ts %>% 
  diff(lag = 12,diffences=D_hidro) %>% 
  diff(diffences=d_hidro) %>% 
  ts_cor(lag.max = 36)

Estimacion del modelo

modelo_estimado_hidro <- data.hidro.ts %>% 
  Arima(order = c(0, 1, 0),
        seasonal = c(1, 1, 1))

summary(modelo_estimado_hidro)

## Series: . 
## ARIMA(0,1,0)(1,1,1)[12] 
## 
## Coefficients:
##          sar1     sma1
##       -0.2376  -0.6941
## s.e.   0.2483   0.3848
## 
## sigma^2 = 1.936e+09:  log likelihood = -622.57
## AIC=1251.14   AICc=1251.65   BIC=1256.94
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 1657.922 38496.49 27289.55 0.09035022 12.73477 0.6830186
##                    ACF1
## Training set -0.4642396

modelo_estimado_hidro %>% autoplot(type="both")+theme_solarized()

modelo_estimado_hidro %>% check_res(lag.max = 36)

data_hidro_Sarima<-modelo_estimado_hidro$fitted
data_hidro_HW<-PronosticosHW_hidro$fitted
grafico_comparativo_hidro<-cbind(data.hidro.ts,data_hidro_Sarima,data_hidro_HW)
ts_plot(grafico_comparativo_hidro)

Verificacion de ajuste

a_hidro<-data.hidro.ts %>% 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_hidro)

## # 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_hidro %>% pivot_longer(everything(), names_to = "Model name",
                         values_to = "Orders") %>% glance() %>% 
  arrange(AICc) ->tabla_3
tabla_3

## # 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_011    Orders     1.65e9   -623. 1250. 1250. 1254. <cpl>    <cpl>   
## 2 arima_original   Orders     1.94e9   -623. 1251. 1252. 1257. <cpl>    <cpl>   
## 3 arima_010_110    Orders     2.34e9   -624. 1253. 1253. 1257. <cpl>    <cpl>   
## 4 arima_automatico Orders     1.10e9   -757. 1517. 1517. 1521. <cpl>    <cpl>

Validacion Cruzada

Yt_hidro<-data.hidro.ts %>% as_tsibble() %>% rename(Hidrocarburos=value)

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

TSCV_hidro<-data.cross.validation %>% 
model(ARIMA(Hidrocarburos ~ pdq(0, 1, 0) + PDQ(1, 1, 1))) %>% 
forecast(h=12) %>% accuracy(Yt_hidro)

print(TSCV_hidro)

## # 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(Hidrocarburos ~~ Test  6704. 52159. 46080.  1.26  18.0  1.15  1.07 0.339