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

Metodología Box-Jenkins

library(readxl)
library(forecast)
IVAE <- read_excel("C:/Users/HP/Desktop/REBE UNIVERSIDAD/ECONOMETRIA/ECONOMETRIA/ejericios/IVAE.xlsx")
             col_types = c("skip", "numeric")
             skip = 5
# Serie
IVAE.ts<- ts(data = IVAE, start = c(2009,1),
                frequency = 12)
Yt<-IVAE.ts

Generacion de Grafica general

library(TSstudio)
library(forecast)

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

Identificación

Orden de Integración

library(kableExtra)
library(magrittr)
d<-ndiffs(Yt)
D<-nsdiffs(Yt)
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 
# Gráfico de la serie diferenciada
Yt %>% 
  diff(lag = 12,diffences=D) %>% 
  diff(diffences=d) %>% 
  ts_plot(title = "Yt estacionaria")

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

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

Estimación del modelo.

Usando forecast

library(forecast)
library(ggthemes)
modelo_estimado <- Yt %>% 
  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.2040  -0.6422
## s.e.   0.1586   0.1554
## 
## sigma^2 = 6.897:  log likelihood = -323.43
## AIC=652.85   AICc=653.04   BIC=661.55
## 
## Training set error measures:
##                      ME     RMSE      MAE        MPE     MAPE      MASE
## Training set 0.04712358 2.488576 1.666188 0.02203353 1.658867 0.4981404
##                    ACF1
## Training set 0.08970896
modelo_estimado %>% autoplot(type="both")+theme_solarized()

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

Importacion de PronosticosHW de la practica anterior

Estimación del modelo

library(forecast)


modeloHW<-HoltWinters(x = Yt,
                      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 = Yt, seasonal = "multiplicative", optim.start = c(0.9,     0.9, 0.9))
## 
## Smoothing parameters:
##  alpha: 0.8408163
##  beta : 0
##  gamma: 1
## 
## Coefficients:
##            [,1]
## a   117.0799442
## b     0.1600306
## s1    0.9502255
## s2    1.0233274
## s3    1.0518154
## s4    0.9900276
## s5    1.0007537
## s6    0.9807088
## s7    0.9600206
## s8    1.0149628
## s9    1.0915423
## s10   0.9796752
## s11   0.9584676
## s12   0.9994880

Generar pronostico

# Generar el pronostico:
pronosticosHW<-forecast(object = modeloHW,h=12, level=c(0.95,0.99))
pronosticosHW
##          Point Forecast     Lo 95    Hi 95     Lo 99    Hi 99
## Apr 2021       111.4044 105.94952 116.8593 104.23546 118.5734
## May 2021       120.1386 112.77972 127.4976 110.46737 129.8099
## Jun 2021       123.6515 114.83356 132.4694 112.06278 135.2401
## Jul 2021       116.5461 107.01096 126.0813 104.01480 129.0774
## Aug 2021       117.9689 107.30373 128.6341 103.95248 131.9854
## Sep 2021       115.7630 104.33417 127.1918 100.74297 130.7830
## Oct 2021       113.4746 101.36814 125.5810  97.56402 129.3852
## Nov 2021       120.1312 106.57272 133.6897 102.31234 137.9500
## Dec 2021       129.3698 114.12877 144.6109 109.33967 149.4000
## Jan 2022       116.2681 101.78191 130.7543  97.23002 135.3062
## Feb 2022       113.9046  98.99575 128.8134  94.31107 133.4981
## Mar 2022       118.9394  81.61739 156.2614  69.88997 167.9888

Grafico

# Gráfico de la serie original y del pronóstico.
 pronosticosHW%>% autoplot()

Yt_Sarima<-modelo_estimado$fitted
Yt_HW<-pronosticosHW$fitted
grafico_comparativo<-cbind(Yt,Yt_Sarima,Yt_HW)
ts_plot(grafico_comparativo)

Verificación de sobre ajuste/sub ajuste

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 %>% 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 × 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.30   -320.  649.  649.  660. <cpl [1]>  <cpl [12]>
## 2 arima_010_011    Orders   6.88   -324.  652.  653.  658. <cpl [0]>  <cpl [12]>
## 3 arima_original   Orders   6.90   -323.  653.  653.  662. <cpl [12]> <cpl [12]>
## 4 arima_010_110    Orders   7.60   -328.  660.  660.  666. <cpl [12]> <cpl [0]>

Validación Cruzada (Cross Validated)

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

TSCV<-data.cross.validation %>% 
model(ARIMA(IVAE ~ 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 2021 abr.
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(IVAE ~ pdq(0, … Test  0.0234  2.98  2.08 -0.0260  2.01 0.623 0.612 0.184