Ejercicio de Pronóstico Modelos SARIMA y Holt-Winters

UNIVERSIDAD DE EL SALVADOR
FACULTAD DE CIENCIAS ECONÓMICAS
ESCUELA DE ECONOMÍA
CICLO I - 2024
Escudo UES
DESCOMPOSICION DE SERIES TEMPORALES
Asignatura:
Econometría
Grupo teórico:
Gt 03
Docente:
MSF. Carlos Ademir Pérez Alas
Estudiante:
Vanessa Iveth López González LG20034
Ciudad Universitaria, 26 de Junio de 2024

PIB trimestral precios corrientes SLV 2009-2021

SE HACE LA CARGA DE DATOS Y CONVERSIÓN DE OBJETO A SERIE TEMPORAL

library(readxl)
PIB_trimestral <- read_excel("C:/Users/iveth/OneDrive/Documents/PIB_trimestral.xlsx", 
    col_types = c("skip", "numeric"), skip = 5)
names(PIB_trimestral)<-c("PIB_trimestral")

# convertir en objeto de Serie Temporal (ts)
PIB_trimestral.ts<-ts(data = PIB_trimestral$PIB_trimestral,start = c(2009,1),frequency = 4) |> print()
##         Qtr1    Qtr2    Qtr3    Qtr4
## 2009 4227.33 4436.39 4312.85 4625.05
## 2010 4375.00 4579.71 4498.45 4994.77
## 2011 4727.87 5163.09 5032.40 5360.43
## 2012 5176.52 5401.29 5225.97 5582.38
## 2013 5208.33 5550.62 5428.02 5803.99
## 2014 5434.08 5701.62 5589.75 5868.02
## 2015 5617.57 5870.22 5828.24 6122.21
## 2016 5695.11 6181.81 5962.28 6352.22
## 2017 5934.20 6271.01 6150.91 6623.07
## 2018 6159.12 6529.22 6480.41 6852.10
## 2019 6401.11 6713.27 6691.85 7074.90
## 2020 6462.10 5378.00 6118.58 6962.52
## 2021 6736.65 7203.03 7236.93 7866.54

1. Generar pronósticos para 2022 completo (diciembre) usando un modelo SARIMA:

a) Usando la metodología paso a paso mostrada en Clase. (manual)

Calculo de las diferencias

library(forecast)
library(kableExtra)
d<-ndiffs(PIB_trimestral.ts)
D<-nsdiffs(PIB_trimestral.ts)
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

Encontramos p,q y P,Q

Respuestas: p=0 y q=0 P=2 y Q=0

SARIMA(p,d,q)(P,D,Q) = SARIMA(0,1,0)(2,1,0)[4]

library(forecast)
PIB_trimestral.arima.manual<-Arima(y = PIB_trimestral.ts,order = c(0,1,0), seasonal = c(2,1,0)) |> print()
## Series: PIB_trimestral.ts 
## ARIMA(0,1,0)(2,1,0)[4] 
## 
## Coefficients:
##          sar1     sar2
##       -0.8752  -0.3604
## s.e.   0.1439   0.2166
## 
## sigma^2 = 75962:  log likelihood = -331.39
## AIC=668.77   AICc=669.33   BIC=674.32

Pronostico manual

library(forecast)

#Tabla
pronostico.arima.manual<-forecast(object = PIB_trimestral.arima.manual,h = 4,level = c(0.95)) |> print()
##         Point Forecast    Lo 95    Hi 95
## 2022 Q1       7360.342 6820.153 7900.532
## 2022 Q2       6972.921 6208.978 7736.865
## 2022 Q3       7350.701 6415.065 8286.337
## 2022 Q4       8001.800 6921.421 9082.180
#Gráfica
PIB_trimestral.arima.manual |> forecast(h = 4,level = c(.95)) |> autoplot()

b) Usando la generación automática disponible en forecast.

library(forecast)
library(TSstudio)
PIB_trimestral.arima.automatico<-auto.arima(y = PIB_trimestral.ts)
summary(PIB_trimestral.arima.automatico)
## Series: PIB_trimestral.ts 
## ARIMA(1,0,0)(2,1,0)[4] with drift 
## 
## Coefficients:
##          ar1     sar1     sar2    drift
##       0.5874  -0.8878  -0.3696  53.5077
## s.e.  0.1216   0.1435   0.2155   9.6501
## 
## sigma^2 = 62491:  log likelihood = -332.83
## AIC=675.66   AICc=677.09   BIC=685.02
## 
## Training set error measures:
##                    ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 1.201582 229.9509 105.3609 -0.0516886 1.814615 0.3266292
##                    ACF1
## Training set 0.04733266

Pronostico automatico

library(forecast)

#Tabla
pronostico.arima.automatico<-forecast(object = PIB_trimestral.arima.automatico,h = 4,level = c(.95)) |> print()
##         Point Forecast    Lo 95    Hi 95
## 2022 Q1       7190.474 6700.516 7680.431
## 2022 Q2       6698.563 6130.341 7266.786
## 2022 Q3       7020.812 6427.981 7613.642
## 2022 Q4       7636.634 7035.547 8237.721
#Gráfica
PIB_trimestral.arima.automatico |> forecast(h = 4,level = c(.95)) |> autoplot()

2. Generar pronósticos para 2022 completo (diciembre) usando un modelo de Holt-Winters:

Usando Stats y forecast

library(forecast)

Yt <- PIB_trimestral.ts #Serie original

#Estimar el modelo
Modelo_HW<-HoltWinters(x = Yt,
                      seasonal = "multiplicative",
                      optim.start = c(0.9,0.9,0.9))
Modelo_HW
## 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.824158
##  beta : 0
##  gamma: 0
## 
## Coefficients:
##            [,1]
## a  7531.9276190
## b    48.1826250
## s1    0.9754034
## s2    1.0055201
## s3    0.9784087
## s4    1.0406678

Generar pronostico

Pronosticos_HW<-forecast(object = Modelo_HW,h = 12,level = c(0.95))
Pronosticos_HW
##         Point Forecast    Lo 95    Hi 95
## 2022 Q1       7393.665 6979.857 7807.474
## 2022 Q2       7670.402 7076.084 8264.719
## 2022 Q3       7510.730 6799.632 8221.828
## 2022 Q4       8038.803 7282.403 8795.203
## 2023 Q1       7581.655 6760.761 8402.549
## 2023 Q2       7864.196 6922.198 8806.194
## 2023 Q3       7699.300 6693.620 8704.979
## 2023 Q4       8239.372 7169.660 9309.083
## 2024 Q1       7769.645 6684.981 8854.309
## 2024 Q2       8057.990 6865.721 9250.260
## 2024 Q3       7887.869 6656.153 9119.584
## 2024 Q4       8439.940 7129.816 9750.064

Gráfico de la serie original y del pronóstico

Pronosticos_HW %>% autoplot()

Usando Forecast (Aproximación por Espacios de los Estados ETS )

Generar pronóstico

library(forecast)

Pronosticos_HW2<-hw(y = Yt,
                   h = 12,
                   level = c(0.95),
                   seasonal = "multiplicative",
                   initial = "optimal")
Pronosticos_HW2
##         Point Forecast    Lo 95     Hi 95
## 2022 Q1       7369.533 6773.650  7965.416
## 2022 Q2       7646.908 6848.590  8445.226
## 2022 Q3       7615.998 6677.689  8554.307
## 2022 Q4       8173.968 7036.020  9311.917
## 2023 Q1       7673.978 6497.332  8850.624
## 2023 Q2       7959.583 6638.054  9281.112
## 2023 Q3       7924.258 6516.656  9331.859
## 2023 Q4       8501.498 6900.278 10102.717
## 2024 Q1       7978.423 6396.112  9560.733
## 2024 Q2       8272.257 6554.329  9990.186
## 2024 Q3       8232.517 6450.310 10014.724
## 2024 Q4       8829.027 6844.057 10813.997

Gráfico de la serie original y del pronóstico

Pronosticos_HW2 %>% autoplot()

3. Realice una validación cruzada para los modelos estimados. Interprete el MAPE

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 = 40,.step = 1)

TSCV1<-data.cross.validation %>% 
model(ARIMA(IVAE ~ pdq(0, 1, 0) + PDQ(2, 1, 0))) %>% 
forecast(h=1) %>% accuracy(Yt)

print(TSCV1)
## # 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(IVAE ~ pdq(0, 1… Test   108.  603.  401.  1.11  6.36  1.24  1.30 -0.0245

Interpretación del MAPE: 6.356362= en promedio, los pronósticos del modelo tienen un error absoluto del 6.356362% en relación con los valores reales de la serie temporal.

Indice de precios al consumidor general SLV 2009-2022 (mayo)

SE HACE LA CARGA DE DATOS Y CONVERSIÓN DE OBJETO A SERIE TEMPORAL

library(readxl)

IPCg_mensual <- read_excel("C:/Users/iveth/OneDrive/Documents/PIB general.xlsx", 
    col_types = c("skip", "numeric"), skip = 5)
names(IPCg_mensual)<-c("IPCg_mensual")

# convertir en objeto de Serie Temporal (ts)
IPCg_mensual.ts<-ts(data = IPCg_mensual$IPCg_mensual,start = c(2009,1),frequency = 12) |> print()
##         Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct
## 2009  99.65  99.53 100.00 100.08 100.36 100.40  99.99  99.79  99.73  99.03
## 2010 100.44 100.57 100.89 100.71 100.50 100.96 100.98 100.81 101.11 101.78
## 2011 102.77 102.96 103.63 106.71 107.23 107.29 107.57 107.69 107.40 107.32
## 2012 107.65 108.00 108.16 108.83 108.53 107.94 107.61 107.80 108.24 108.33
## 2013 108.59 109.05 109.54 108.85 108.69 108.91 108.78 108.87 109.06 108.91
## 2014 109.51 109.72 109.99 109.47 109.72 110.15 110.77 111.04 110.91 110.97
## 2015 108.69 108.56 109.10 109.11 109.33 109.24 109.16 108.82 108.41 110.76
## 2016 110.67 110.37 110.32 110.05 110.13 110.24 110.12 109.85 109.51 109.79
## 2017 110.39 110.69 110.92 111.00 111.19 111.26 111.24 111.10 111.22 111.36
## 2018 111.96 112.05 111.93 111.97 112.11 112.26 112.42 112.71 112.76 113.02
## 2019 112.24 112.44 112.69 112.87 113.01 112.85 112.56 112.16 111.99 112.04
## 2020 112.15 112.00 112.09 111.69 111.94 112.59 112.49 111.82 111.56 111.81
## 2021 112.49 113.19 114.08 114.81 114.84 115.51 116.36 116.63 117.10 117.95
## 2022 119.79 120.73 121.71 122.32 123.43                                   
##         Nov    Dec
## 2009 100.41 100.00
## 2010 102.24 102.13
## 2011 107.48 107.29
## 2012 108.18 108.13
## 2013 109.00 108.98
## 2014 110.40 109.50
## 2015 110.69 110.61
## 2016 109.78 109.58
## 2017 111.62 111.81
## 2018 112.82 112.30
## 2019 112.17 112.29
## 2020 111.98 112.20
## 2021 118.92 119.06
## 2022

1. Generar pronósticos para 2022 completo (diciembre) usando un modelo SARIMA:

a) Usando la metodología paso a paso mostrada en Clase. (manual)

Calculo de las diferencias

library(forecast)
library(kableExtra)

d<-ndiffs(IPCg_mensual.ts)
D<-nsdiffs(IPCg_mensual.ts)

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 0

Encontramos p,q y P,Q

library(TSstudio)

IPCg_mensual.ts_diff<-IPCg_mensual.ts |>  
  diff(lag = 12,diffences=D) |> 
  diff(diffences=d)
IPCg_mensual.ts_diff |> 
  ts_cor(lag.max = 3*12) # si fueran datos mensuales sería 3X12

Respuestas: p=1 y q=0 P=2 y Q=0

SARIMA(p,d,q)(P,D,Q) = SARIMA(0,1,0)(1,0,0)[12]

library(forecast)
IPCg_mensual.arima.manual<-Arima(y = IPCg_mensual.ts,order = c(1,1,0), seasonal = c(1,0,0)) |> print()
## Series: IPCg_mensual.ts 
## ARIMA(1,1,0)(1,0,0)[12] 
## 
## Coefficients:
##          ar1    sar1
##       0.3059  0.1928
## s.e.  0.0763  0.0849
## 
## sigma^2 = 0.2266:  log likelihood = -107.5
## AIC=221   AICc=221.15   BIC=230.23

Pronostico manual

library(forecast)

#Tabla
pronostico.arima.manual<-forecast(object = IPCg_mensual.arima.manual,h = 12,level = c(0.95)) |> print()
##          Point Forecast    Lo 95    Hi 95
## Jun 2022       123.8969 122.9640 124.8298
## Jul 2022       124.1641 122.6296 125.6986
## Aug 2022       124.2478 122.2330 126.2625
## Sep 2022       124.3480 121.9326 126.7634
## Oct 2022       124.5148 121.7524 127.2773
## Nov 2022       124.7027 121.6311 127.7743
## Dec 2022       124.7300 121.3773 128.0826
## Jan 2023       124.8708 121.2588 128.4828
## Feb 2023       125.0520 121.1981 128.9059
## Mar 2023       125.2409 121.1593 129.3225
## Apr 2023       125.3585 121.0613 129.6557
## May 2023       125.5725 121.0700 130.0750
#Gráfica
IPCg_mensual.arima.manual |> forecast(h = 12,level = c(.95)) |> autoplot()

b) Usando la generación automática disponible en forecast.

library(forecast)
library(TSstudio)
IPCg_mensual.arima.automatico<-auto.arima(y = IPCg_mensual.ts)
summary(IPCg_mensual.arima.automatico)
## Series: IPCg_mensual.ts 
## ARIMA(0,1,1)(1,0,0)[12] with drift 
## 
## Coefficients:
##          ma1    sar1   drift
##       0.2321  0.1481  0.1557
## s.e.  0.0758  0.0861  0.0526
## 
## sigma^2 = 0.2203:  log likelihood = -104.62
## AIC=217.24   AICc=217.5   BIC=229.54
## 
## Training set error measures:
##                       ME      RMSE       MAE           MPE      MAPE    MASE
## Training set 0.001979736 0.4634631 0.3022345 -0.0008175962 0.2770172 0.16378
##                    ACF1
## Training set 0.01394664

Pronóstico automático

library(forecast)

#Tabla
pronostico.arima.automatico<-forecast(object = IPCg_mensual.arima.automatico,h = 12,level = c(.95)) |> print()
##          Point Forecast    Lo 95    Hi 95
## Jun 2022       123.8751 122.9552 124.7950
## Jul 2022       124.1337 122.6740 125.5933
## Aug 2022       124.3063 122.4583 126.1543
## Sep 2022       124.5086 122.3408 126.6765
## Oct 2022       124.7672 122.3210 127.2134
## Nov 2022       125.0435 122.3475 127.7396
## Dec 2022       125.1970 122.2724 128.1215
## Jan 2023       125.4378 122.3013 128.5742
## Feb 2023       125.7097 122.3747 129.0446
## Mar 2023       125.9875 122.4652 129.5097
## Apr 2023       126.2105 122.5104 129.9106
## May 2023       126.5076 122.6378 130.3774
#Gráfica
IPCg_mensual.arima.automatico |> forecast(h = 12,level = c(.95)) |> autoplot()

2. Generar pronósticos para 2022 completo (diciembre) usando un modelo de Holt-Winters:

Usando Stats y forecast

library(forecast)

Yt <- IPCg_mensual.ts #Serie original

#Estimar el modelo
Modelo_HW<-HoltWinters(x = Yt,
                      seasonal = "multiplicative",
                      optim.start = c(0.9,0.9,0.9))
Modelo_HW
## 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.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 pronostico

Pronosticos_HW<-forecast(object = Modelo_HW,h = 7,level = c(0.95))
Pronosticos_HW
##          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

Gráfico de la serie original y del pronóstico

Pronosticos_HW %>% autoplot()

Usando Forecast (Aproximación por Espacios de los Estados ETS )

library(forecast)
#Generar el pronóstico:
Pronosticos_HW2<-hw(y = Yt,
                   h = 7,
                   level = c(0.95),
                   seasonal = "multiplicative",
                   initial = "optimal")
Pronosticos_HW2
##          Point Forecast    Lo 95    Hi 95
## Jun 2022       123.8302 121.9984 125.6620
## Jul 2022       124.5773 122.6371 126.5174
## Aug 2022       125.1928 123.0924 127.2932
## Sep 2022       125.6452 123.3312 127.9591
## Oct 2022       126.4964 123.9071 129.0858
## Nov 2022       127.6980 124.7742 130.6219
## Dec 2022       127.9518 124.6694 131.2342

Gráfico de la serie original y del pronóstico.

Pronosticos_HW2 %>% autoplot()

3. Realice una validación cruzada para los modelos estimados. Interprete el MAPE

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 = 40,.step = 1)

TSCV1<-data.cross.validation %>% 
model(ARIMA(IVAE ~ pdq(0, 1, 1) + PDQ(1, 0, 0))) %>% 
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(TSCV1)
## # 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(IVAE ~ pdq(0… Test  0.00319 0.424 0.292 -0.00139 0.261 0.158 0.156 0.206

Interpretación del MAPE: 0.2610455 = en promedio, los pronósticos del modelo tienen un error absoluto del 0.26% en relación con los valores reales de la serie temporal.

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

SE HACE LA CARGA DE DATOS Y CONVERSIÓN DE OBJETO A SERIE TEMPORAL

library(readxl)

ImpHidro_m<- read_excel("C:/Users/iveth/OneDrive/Documents/Importacion de hidrocarburos.xlsx", 
    col_types = c("skip", "numeric"), skip = 5)
names(ImpHidro_m)<-c("ImpHidro_m")

# convertir en objeto de Serie Temporal (ts)
ImpHidro_m.ts<-ts(data = ImpHidro_m$ImpHidro_m,start = c(2017,1),frequency = 12) |> print()
##            Jan       Feb       Mar       Apr       May       Jun       Jul
## 2017 245281.10 173399.96 233320.45 220284.95 203655.71 275330.95 176511.44
## 2018 176860.00 208040.00 241980.00 267580.00 211080.00 202470.00 199160.00
## 2019 244220.00 172590.00 227640.00 232750.00 268910.00 214990.00 228650.00
## 2020 209088.61 192750.84 232738.11 212858.25  97575.83 175258.65 190317.69
## 2021 239938.21 190943.14 252235.78 175849.15 211409.80 275893.64 238569.80
## 2022 207025.43 258770.90 288743.53 284014.55 171331.42                    
##            Aug       Sep       Oct       Nov       Dec
## 2017 181637.97 181178.43 162342.83 273551.34 159754.63
## 2018 229540.00 219830.00 246640.00 206880.00 220810.00
## 2019 228580.00 198920.00 236590.00 233050.00 231570.00
## 2020 134017.33 187064.26 161027.78 221739.41 185403.23
## 2021 196986.17 207613.25 227286.37 225971.47 280510.69
## 2022

1. Generar pronósticos para 2022 completo (diciembre) usando un modelo SARIMA:

a) Usando la metodología paso a paso mostrada en Clase. (manual)

Calculo de las diferencias

library(forecast)
library(kableExtra)

d<-ndiffs(ImpHidro_m.ts)
D<-nsdiffs(ImpHidro_m.ts)

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 0
D 0

Encontramos p,q y P,Q

library(TSstudio)
ImpHidro_m.ts_diff<-ImpHidro_m.ts |>  
  diff(lag = 12,diffences=D) |> 
  diff(diffences=d)
ImpHidro_m.ts_diff |> 
  ts_cor(lag.max = 3*12) # si fueran datos mensuales sería 3X12

Respuestas: p=1 y q=1 P=0 y Q=1

SARIMA(p,d,q)(P,D,Q) = SARIMA(1,0,1)(0,0,1)[12]

library(forecast)
ImpHidro_m.arima.manual<-Arima(y = ImpHidro_m.ts,order = c(1,0,1), seasonal = c(0,0,1)) |> print()
## Series: ImpHidro_m.ts 
## ARIMA(1,0,1)(0,0,1)[12] with non-zero mean 
## 
## Coefficients:
##           ar1     ma1     sma1        mean
##       -0.8285  0.9116  -0.1688  214454.867
## s.e.   0.1718  0.1224   0.1267    4059.913
## 
## sigma^2 = 1.392e+09:  log likelihood = -774.69
## AIC=1559.37   AICc=1560.39   BIC=1570.24

Pronostico manual

library(forecast)

#Tabla
pronostico.arima.manual<-forecast(object = ImpHidro_m.arima.manual,h = 12,level = c(0.95)) |> print()
##          Point Forecast    Lo 95    Hi 95
## Jun 2022       200101.5 126987.6 273215.4
## Jul 2022       215246.3 141880.6 288612.0
## Aug 2022       216118.7 142580.6 289656.7
## Sep 2022       219364.9 145708.7 293021.0
## Oct 2022       211282.7 137545.6 285019.7
## Nov 2022       214203.3 140410.7 287995.9
## Dec 2022       202403.6 128572.9 276234.3
## Jan 2023       216355.5 142498.7 290212.3
## Feb 2023       206652.7 132778.0 280527.5
## Mar 2023       201675.1 127788.0 275562.1
## Apr 2023       203022.3 129126.8 276917.8
## May 2023       222982.6 149081.3 296883.9
#Gráfica
ImpHidro_m.arima.manual |> forecast(h = 12,level = c(.95)) |> autoplot()

b) Usando la generación automática disponible en forecast.

library(forecast)
library(TSstudio)
ImpHidro_m.arima.automatico<-auto.arima(y = ImpHidro_m.ts)
summary(ImpHidro_m.arima.automatico)
## Series: ImpHidro_m.ts 
## ARIMA(0,0,0)(1,0,0)[12] with non-zero mean 
## 
## Coefficients:
##          sar1        mean
##       -0.2204  214330.441
## s.e.   0.1437    3836.624
## 
## sigma^2 = 1.362e+09:  log likelihood = -775.07
## AIC=1556.14   AICc=1556.53   BIC=1562.66
## 
## Training set error measures:
##                     ME     RMSE      MAE     MPE     MAPE      MASE       ACF1
## Training set -258.2539 36337.82 29355.11 -3.5468 14.88772 0.6657831 0.08935384
ImpHidro_m.arima.automatico |> forecast(h = 12,level = c(.95)) |> autoplot()

2. Generar pronósticos para 2022 completo (diciembre) usando un modelo de Holt-Winters:

Usando Stats y forecast

library(forecast)

Yt <- ImpHidro_m.ts #Serie original

#Estimar el modelo
Modelo_HW<-HoltWinters(x = Yt,
                      seasonal = "multiplicative",
                      optim.start = c(0.9,0.9,0.9))
Modelo_HW
## 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.1429275
##  beta : 0
##  gamma: 0.4446265
## 
## Coefficients:
##             [,1]
## a   2.521978e+05
## b   1.239453e+03
## s1  1.031430e+00
## s2  9.655959e-01
## s3  8.373676e-01
## s4  8.997235e-01
## s5  9.154143e-01
## s6  1.034878e+00
## s7  1.012064e+00
## s8  9.427592e-01
## s9  9.466304e-01
## s10 1.117844e+00
## s11 1.021463e+00
## s12 7.895019e-01

Generar el pronostico

Pronosticos_HW<-forecast(object = Modelo_HW,h = 7,level = c(0.95))
Pronosticos_HW
##          Point Forecast    Lo 95    Hi 95
## Jun 2022       261402.7 215838.4 306967.1
## Jul 2022       245914.8 198888.4 292941.2
## Aug 2022       214295.9 166399.4 262192.4
## Sep 2022       231369.0 181308.4 281429.6
## Oct 2022       236538.6 184759.4 288317.8
## Nov 2022       268690.0 213493.2 323886.9
## Dec 2022       264021.2 207877.7 320164.8

Gráfico de la serie original y del pronóstico.

Pronosticos_HW %>% autoplot()

Usando Forecast (Aproximación por Espacios de los Estados ETS )

library(forecast)

#Generar el pronóstico:
Pronosticos_HW2<-hw(y = Yt,
                   h = 7,
                   level = c(0.95),
                   seasonal = "multiplicative",
                   initial = "optimal")
Pronosticos_HW2
##          Point Forecast    Lo 95    Hi 95
## Jun 2022       255125.1 164897.2 345353.1
## Jul 2022       235290.6 150766.7 319814.5
## Aug 2022       221427.7 140674.9 302180.4
## Sep 2022       226003.7 142372.1 309635.3
## Oct 2022       234133.6 146263.6 322003.5
## Nov 2022       263006.3 162944.2 363068.5
## Dec 2022       247224.6 151914.6 342534.6

Gráfico de la serie original y del pronóstico.

Pronosticos_HW2 %>% autoplot()

3. Realice una validación cruzada para los modelos estimados. Interprete el MAPE

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 = 40,.step = 1)

TSCV1<-data.cross.validation %>% 
model(ARIMA(IVAE ~ pdq(0, 0, 0) + PDQ(1, 0, 0))) %>% 
forecast(h=1) %>% accuracy(Yt)

print(TSCV1)
## # 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(IVAE ~ pdq(0, … Test  -2843. 46885. 39050. -7.46  21.1 0.886 0.843 0.231

Interpretación del MAPE: 21.11083 = en promedio, los pronósticos del modelo tienen un error absoluto del 21.11% en relación con los valores reales de la serie temporal.