GEM Perú
Diplomado en Econometría Aplicada
Econometria Financiera



Trabajo Final
Modelo SARIMA

Introduccion

Según se asevera en un reporte preliminar del Instituto Nacional de Estadística (INE) en Bolivia, la economía del país andino ha acumulado un crecimiento interanual del 9,67% en el período correspondiente a enero-agosto del 2021, en franca recuperación en medio de la pandemia global del coronavirus. El INE señala que el Índice Global de la Actividad Económica (IGAE), de enero a agosto de 2021, registró una tasa de variación acumulada positiva de 9,67% con relación a similar período de 2020», resaltando, que el desempeño económico reflejado en el IGAE en este período estuvo impulsado principalmente por los sectores de minería, que se expandió 55,57%, transporte 40,06% y construcción 33,72%.

Desarrollo

Definición de la serie de tiempo.

El Índice Global de Actividad Económica (IGAE) y sus variaciones son estimaciones anticipadas y aproximadas de la evolución de la actividad económica, resultado del agregado de índices de producción sectorial, se diferencia del Producto Interno Bruto (PIB) que es resultante de un balance macroeconómico contable entre la oferta y demanda a precios constantes

Importación de los datos

library(readxl)

## Warning: package 'readxl' was built under R version 4.2.1

ruta <- "F:/Econometria Aplicada/Econometria Financiera"

IGAE <- read_excel(paste0(ruta, "/IGAE.xlsx"))

Visualización de los datos

Librerias a utilizar

• library(tidyverse)
• library(zoo)
• library(xts)
• library(lubridate)

IGAE <- IGAE %>% 
  mutate(Fecha = seq(as.Date("2008-01-01", format="%Y-%m-%d"),
                        as.Date("2022-05-01", format="%Y-%m-%d"),
                        by= "1 month"),
         Mes = as.yearmon(Fecha, "%m"))

library(scales)

## 
## Attaching package: 'scales'

## The following object is masked from 'package:purrr':
## 
##     discard

## The following object is masked from 'package:readr':
## 
##     col_factor

## Gráfico con escala real

IGAE %>%
  ggplot(aes(x=Mes, y=Act_Eco))+
  geom_line()+
  scale_x_yearmon(format = "%Y-%m", n=50)+
  scale_y_continuous(labels = comma,
                     breaks = seq())+
  labs(title = 'Indice Global de Actividad Economica de 2008-2022',
       subtitle = 'Serie Mensual',
       caption = 'Elaboración propia con datos de INE',
       x='Mes',
       y=expression(IGAE[t]))+
  theme(axis.text.x = element_text(angle=90,
                                   hjust=1),
        axis.title.y = element_text(hjust=1,
                                    vjust=0.5,
                                    angle=0))

Identificación de estacionariedad

Declaración de serie de tiempo

IGAE_ts <- ts(IGAE$Act_Eco, start=c(2008,1), frequency = 12)
IGAE_ts

##         Jan    Feb    Mar    Apr    May    Jun    Jul    Aug    Sep    Oct
## 2008 173.16 167.19 190.98 207.59 206.50 204.06 201.22 193.40 202.01 207.98
## 2009 178.47 173.29 195.24 211.58 210.36 209.81 206.19 199.90 212.18 216.77
## 2010 183.89 177.81 202.92 216.94 219.62 219.08 211.36 206.23 223.63 228.56
## 2011 197.77 190.96 210.78 228.85 229.30 225.42 222.23 217.94 234.53 238.54
## 2012 205.62 198.87 225.06 240.07 239.26 234.23 231.20 228.37 246.13 254.05
## 2013 221.95 212.45 238.29 257.09 254.08 252.94 250.11 247.46 259.60 271.24
## 2014 235.85 226.97 249.61 269.96 266.92 263.88 263.40 261.89 277.50 285.48
## 2015 247.48 234.04 265.15 283.67 279.66 279.48 270.50 274.23 289.61 300.88
## 2016 262.24 247.74 276.66 291.97 288.49 289.60 282.72 287.93 304.93 312.26
## 2017 268.33 253.88 290.71 305.32 299.15 298.55 293.95 298.13 320.97 327.51
## 2018 283.48 266.50 303.00 318.34 311.48 316.98 303.28 313.25 332.88 340.12
## 2019 293.12 275.73 310.47 325.19 319.06 327.56 312.07 318.08 340.39 347.43
## 2020 305.48 277.78 297.85 232.74 231.14 267.51 272.20 278.48 314.40 345.37
## 2021 291.29 268.82 315.60 306.50 292.62 301.17 297.39 296.19 318.92 343.76
## 2022 302.29 284.19 324.01 317.24 304.53                                   
##         Nov    Dec
## 2008 199.81 198.82
## 2009 209.70 208.20
## 2010 220.02 222.00
## 2011 231.96 235.53
## 2012 247.10 250.32
## 2013 257.89 267.50
## 2014 272.27 280.15
## 2015 285.29 297.09
## 2016 299.02 304.53
## 2017 315.06 321.18
## 2018 324.80 330.36
## 2019 315.43 342.94
## 2020 332.71 337.39
## 2021 333.68 340.38
## 2022

Pruebas Dickey-Fuller aumentada

library(tseries)

## Warning: package 'tseries' was built under R version 4.2.1

## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo

ADF_test <- function(x, alpha=0.05){
  # Requiere tseries y tidyverse
  adf<-adf.test(x)
  df<-data.frame(adf[1],adf[2],adf[4])
  names(df)<-c("Estadístico", "Rezagos", "Probabilidad")
  df<-df %>% 
    mutate(Decisión=ifelse(Probabilidad>alpha, "No rechazar H0",
                           "Rechazar H0"),
           Conclusión =ifelse(Probabilidad>alpha, "No estacionariedad",
                              "Estacionariedad"))
  df
}

adf<-ADF_test(IGAE_ts, alpha=0.05)
print(adf)

##               Estadístico Rezagos Probabilidad    Decisión      Conclusión
## Dickey-Fuller   -3.476656       5   0.04664796 Rechazar H0 Estacionariedad

Pueba URCA para de Dickey-Fuller para estimacion de estadisticos

#### Urca ----

library(urca)

ADF <- ur.df(y=IGAE_ts,
             type='trend',
             lags=5,
             selectlags = 'AIC')

summary(ADF)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression trend 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 + 1 + tt + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -66.068  -6.195   0.013  11.797  37.808 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 55.71570   15.69984   3.549 0.000508 ***
## z.lag.1     -0.28035    0.08127  -3.449 0.000718 ***
## tt           0.22027    0.07090   3.107 0.002237 ** 
## z.diff.lag1  0.18252    0.08400   2.173 0.031258 *  
## z.diff.lag2 -0.29591    0.08408  -3.520 0.000563 ***
## z.diff.lag3  0.14158    0.07511   1.885 0.061227 .  
## z.diff.lag4 -0.37052    0.07521  -4.927 2.07e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 14.78 on 160 degrees of freedom
## Multiple R-squared:  0.3536, Adjusted R-squared:  0.3293 
## F-statistic: 14.59 on 6 and 160 DF,  p-value: 3.016e-13
## 
## 
## Value of test-statistic is: -3.4495 4.451 6.041 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau3 -3.99 -3.43 -3.13
## phi2  6.22  4.75  4.07
## phi3  8.43  6.49  5.47

Interpretacion

Podemos evidenciar que los estadisticos de urca de la prueba de Dickey-Fuller estan por debajo del valor critico que este caso trabajaresmos al 5%, por lo que rechazamos la hipotesis nula de no estacionariedad y concluimos que el modelo es estacionario, siendo esta una conclucion muy cobtradictoria, ya que podemos evidenciar el grafico que los datos no siguen un patron estacionario.

Pruebas Phillips-Perron

PP_test <- function(x, alpha=0.05){
  # Requiere tseries y tidyverse
  pp<-pp.test(x)
  df<-data.frame(pp[1],pp[2],pp[4])
  names(df)<-c("Estadístico", "Rezagos", "Probabilidad")
  df<-df %>% 
    mutate(Decisión=ifelse(Probabilidad>alpha, "No rechazar H0",
                           "Rechazar H0"),
           Conclusión =ifelse(Probabilidad>alpha, "No estacionariedad",
                              "Estacionariedad"))
  df
}

PP <-PP_test(IGAE_ts)

## Warning in pp.test(x): p-value smaller than printed p-value

print(PP)

##                        Estadístico Rezagos Probabilidad    Decisión
## Dickey-Fuller Z(alpha)    -58.1666       4         0.01 Rechazar H0
##                             Conclusión
## Dickey-Fuller Z(alpha) Estacionariedad

Pueba URCA para de Phillips-Perron para estimacion de estadisticos

PP <- ur.pp(x=IGAE_ts,
            type='Z-tau',
            model='trend',
            lags='short')

summary(PP)

## 
## ################################## 
## # Phillips-Perron Unit Root Test # 
## ################################## 
## 
## Test regression with intercept and trend 
## 
## 
## Call:
## lm(formula = y ~ y.l1 + trend)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -70.573  -8.280   0.892  11.918  29.713 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 94.37739   15.59195   6.053 8.91e-09 ***
## y.l1         0.64380    0.05914  10.886  < 2e-16 ***
## trend        0.28162    0.05454   5.163 6.77e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 16.35 on 169 degrees of freedom
## Multiple R-squared:  0.8725, Adjusted R-squared:  0.871 
## F-statistic: 578.2 on 2 and 169 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic, type: Z-tau  is: -5.8991 
## 
##            aux. Z statistics
## Z-tau-mu              6.5617
## Z-tau-beta            5.0296
## 
## Critical values for Z statistics: 
##                      1pct      5pct     10pct
## critical values -4.013968 -3.436685 -3.142214

Interpretacion

Podemos evidenciar que el estadistico de urca de la prueba de Phillips-Perron esta por debajo del valor critico que este caso trabajaresmos al 5%, por lo que rechazamos la hipotesis nula de no estacionariedad y concluimos que el modelo es estacionario, siendo esta una conclucion cobtradictoria, ya que podemos evidenciar el grafico que los datos no siguen un patron estacionario.

Pruebas KPSS

KPSS_test <- function(x, alpha=0.05, type='Trend'){
  # Requiere tseries y tidyverse
  kpss<-kpss.test(x, null=type)
  df<-data.frame(kpss[1],kpss[2],kpss[3])
  names(df)<-c("Estadístico", "Rezagos", "Probabilidad")
  df<-df %>% 
    mutate(Decisión=ifelse(Probabilidad>alpha, "No rechazar H0",
                           "Rechazar H0"),
           Conclusión =ifelse(Probabilidad>alpha, "Estacionariedad",
                              "No estacionariedad"))
  df
}

KPSS<-KPSS_test(IGAE_ts)

## Warning in kpss.test(x, null = type): p-value smaller than printed p-value

print(KPSS)

##            Estadístico Rezagos Probabilidad    Decisión         Conclusión
## KPSS Trend   0.3251277       4         0.01 Rechazar H0 No estacionariedad

Pueba URCA para KPSS para estimacion de estadisticos

KPSS <- ur.kpss(y=IGAE_ts,
            type='mu',
            lags='short')

summary(KPSS)

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: mu with 4 lags. 
## 
## Value of test-statistic is: 3.2563 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.347 0.463  0.574 0.739

Interpretacion

Podemos evidenciar que el estadistico de urca de la prueba de KPSS esta por encima del valor critico que este caso trabajaresmos al 5%, por lo que aceptamos la hipotesis nula de no estacionariedad y concluimos que el modelo es no estacionario, siendo esta una conclucion acertada, ya que podemos evidenciar en el grafico que los datos no siguen un patron estacionario.

Gráfico de Autocorrelacion Simple

library(forecast)

## Warning: package 'forecast' was built under R version 4.2.1

### ACF ----

Acf(IGAE_ts, plot=F) %>% 
  autoplot()+
  labs(title = 'Función de autocorrelación simple',
       subtitle = 'Serie: IGAE en I(0)',
       x='Rezagos',
       y=expression(rho[t]))

Declaracion en primera diferencia de la serie de tiempo

IGAE_d <- diff(IGAE_ts, differences = 1)

IGAE_diff <- data.frame(Fecha=seq(as.Date("2008-02-01", format='%Y-%m-%d'),
                                 as.Date("2022-05-01", format='%Y-%m-%d'),
                                 by='1 month'),
                       IGAE_d) %>% 
  mutate(Mes = as.yearmon(Fecha, "%m"))

Visualización de la diferencia estacional

IGAE_diff %>%
  ggplot(aes(x=Mes, y=IGAE_d))+
  geom_line()+
  scale_x_yearmon(format = "%Y-%m", n=50)+
  scale_y_continuous(labels = comma,
                     breaks = seq())+
  labs(title = 'Indice Global de Actividad Economica de 2008-2022',
       subtitle = 'Serie Mensual',
       caption = 'Elaboración propia con datos de INE',
       x='Mes',
       y=expression(PIB[t]))+
  theme(axis.text.x = element_text(angle=90,
                                   hjust=1),
        axis.title.y = element_text(hjust=1,
                                    vjust=0.5,
                                    angle=0))

Pruebas Dickey-Fuller aumentada

adf<-ADF_test(IGAE_d, alpha=0.05)

## Warning in adf.test(x): p-value smaller than printed p-value

print(adf)

##               Estadístico Rezagos Probabilidad    Decisión      Conclusión
## Dickey-Fuller   -6.116418       5         0.01 Rechazar H0 Estacionariedad

Pueba URCA para de Dickey-Fuller para estimacion de estadisticos

#### Urca ----
adf <- ur.df(y=IGAE_d,
             type='none',
             lags=5,
             selectlags = 'AIC')

summary(adf)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -65.905  -4.255   2.648  10.506  59.159 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## z.lag.1     -1.59882    0.26332  -6.072 8.87e-09 ***
## z.diff.lag1  0.60538    0.22771   2.659  0.00864 ** 
## z.diff.lag2  0.22169    0.18641   1.189  0.23609    
## z.diff.lag3  0.24473    0.15383   1.591  0.11360    
## z.diff.lag4 -0.15638    0.11142  -1.404  0.16238    
## z.diff.lag5 -0.17768    0.07939  -2.238  0.02660 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.13 on 160 degrees of freedom
## Multiple R-squared:  0.6623, Adjusted R-squared:  0.6497 
## F-statistic: 52.31 on 6 and 160 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -6.0718 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

Interpretacion

Podemos evidenciar que el estadistico de urca de la prueba de Dickey-Fuller esta por debajo del valor critico que este caso trabajaresmos al 5%, por lo que rechazamos la hipotesis nula de no estacionariedad y concluimos que el modelo es estacionario.

Pruebas Phillips-Perron

PP <-PP_test(IGAE_d)

## Warning in pp.test(x): p-value smaller than printed p-value

print(PP)

##                        Estadístico Rezagos Probabilidad    Decisión
## Dickey-Fuller Z(alpha)    -129.635       4         0.01 Rechazar H0
##                             Conclusión
## Dickey-Fuller Z(alpha) Estacionariedad

Pueba URCA para de Phillips-Perron para estimacion de estadisticos

PP <- ur.pp(x=IGAE_d,
            type='Z-tau',
            model='constant',
            lags='short')
summary(PP)

## 
## ################################## 
## # Phillips-Perron Unit Root Test # 
## ################################## 
## 
## Test regression with intercept 
## 
## 
## Call:
## lm(formula = y ~ y.l1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -65.877  -8.578  -1.416  11.397  45.933 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)
## (Intercept)  0.80473    1.37928   0.583    0.560
## y.l1        -0.00187    0.07702  -0.024    0.981
## 
## Residual standard error: 18.02 on 169 degrees of freedom
## Multiple R-squared:  3.487e-06,  Adjusted R-squared:  -0.005914 
## F-statistic: 0.0005893 on 1 and 169 DF,  p-value: 0.9807
## 
## 
## Value of test-statistic, type: Z-tau  is: -13.7138 
## 
##          aux. Z statistics
## Z-tau-mu            0.6153
## 
## Critical values for Z statistics: 
##                      1pct      5pct     10pct
## critical values -3.469582 -2.878398 -2.575663

Interpretacion

Podemos evidenciar que el estadistico de urca de la prueba de Phillips-Perron esta por debajo del valor critico que este caso trabajaresmos al 5%, por lo que rechazamos la hipotesis nula de no estacionariedad y concluimos que el modelo es estacionario.

Pruebas KPSS

KPSS<-KPSS_test(IGAE_d,type='Level')

## Warning in kpss.test(x, null = type): p-value greater than printed p-value

print(KPSS)

##            Estadístico Rezagos Probabilidad       Decisión      Conclusión
## KPSS Level  0.03828675       4          0.1 No rechazar H0 Estacionariedad

Pueba URCA para de KPSS para estimacion de estadisticos

kpss <- ur.kpss(y=IGAE_d,
                type='tau',
                lags='short')
summary(kpss)

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: tau with 4 lags. 
## 
## Value of test-statistic is: 0.0148 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.119 0.146  0.176 0.216

Interpretacion

Podemos evidenciar que el estadistico de urca de la prueba de KPSS esta por debajo del valor critico que este caso trabajaresmos al 5%, por lo que rechazamos la hipotesis nula de no estacionariedad y concluimos que el modelo es estacionario.

Gráfico de Autocorrelacion Simple

Acf(IGAE_d, plot=F) %>% 
  autoplot()+
  labs(title = 'Función de autocorrelación simple',
       subtitle = 'Serie: IGAE en I(1)',
       x='Rezagos',
       y=expression(rho[t]))

Pacf(IGAE_d, plot=F) %>% 
  autoplot()+
  labs(title = 'Función de autocorrelación parcial',
       subtitle = 'Serie: PIB en I(1)',
       x='Rezagos',
       y=expression(rho[t]))

Selección de ordenes

p, d=1 ,q

Para el orden p (se utiliza PACF): 2 y 4 Para el orden q (se utiliza ACF): 1-5

ARIMA(2,1,1) ARIMA(2,1,2) ARIMA(2,1,3) ARIMA(2,1,4) ARIMA(2,1,5) ARIMA(4,1,1) ARIMA(4,1,2) ARIMA(4,1,3) ARIMA(4,1,4) ARIMA(4,1,5)

Identificación estacional

Gráfico estacional

library(forecast)
### Se requiere la libreria forecast

ggseasonplot(IGAE_ts, year.labels = TRUE,
             year.labels.left = TRUE)+
  labs(title = 'Gráfico estacional del IGAE',
       subtitle = 'Serie mensual por cortes anuales de 2008-2022',
       x='Mes',
       y= 'PIB')

ggseasonplot(IGAE_d, year.labels = TRUE,
             year.labels.left = TRUE)+
  labs(title = 'Gráfico estacional de la primera diferencia del IGAE',
       subtitle = 'Serie mensual por cortes anuales del 2008-2022',
       x='Mes',
       y= expression(paste(Delta, PIB[t])))

Diferencia estacional

Formato ts

IGAE_s <- diff(IGAE_ts, lag=12)

Formato en data.frame

IGAE_ds <- data.frame(Fecha=seq(as.Date("2009-01-01", format='%Y-%m-%d'),
                               as.Date("2022-05-01", format='%Y-%m-%d'),
                               by='1 month'),
                       IGAE_s) %>% 
  mutate(Mes = as.yearmon(Fecha, "%m"))

Visualización de la diferencia estacional

IGAE_ds %>% 
  ggplot(aes(x=Mes, y=IGAE_s))+
  geom_line()+
  scale_x_yearmon(format = "%Y-%m", n=40)+
  scale_y_continuous(labels = comma)+
  labs(title = 'Diferencia estacional: Indice Global de Actividad Economica 2009-2022',
       subtitle = 'Serie Mensual',
       caption = 'Elaboración propia con datos de INE',
       x='Mes',
       y=expression(paste(Delta[s], IGAE[t])))+
  theme(axis.text.x = element_text(angle=90,
                                   hjust=1),
        axis.title.y = element_text(hjust=1,
                                    vjust=0.5,
                                    angle=0)) 

Prueba DFS

Rezago estacional

IGAE_lag12 <- stats::lag(IGAE_ts, k=-12)

Regresion lineal

Orden 1 estacional

DFS <- lm(IGAE_s ~ IGAE_lag12[1:161])
summary(DFS)

## 
## Call:
## lm(formula = IGAE_s ~ IGAE_lag12[1:161])
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -94.458  -3.603   2.533   6.686  62.214 
## 
## Coefficients:
##                   Estimate Std. Error t value Pr(>|t|)    
## (Intercept)        35.5571     7.3016   4.870 2.68e-06 ***
## IGAE_lag12[1:161]  -0.1032     0.0277  -3.724 0.000272 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 15.81 on 159 degrees of freedom
## Multiple R-squared:  0.08023,    Adjusted R-squared:  0.07444 
## F-statistic: 13.87 on 1 and 159 DF,  p-value: 0.0002717

HEGY

library(uroot)

hegy<-hegy.test(IGAE_ts)
summary(hegy)

## 
##  HEGY test for unit roots
## 
## data:  IGAE_ts
## 
## Fitted model
## ------------
## 
## Call:
## lm(formula = dx ~ 0 + ypi + xreg)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -71.858  -3.819  -0.049   7.318  26.057 
## 
## Coefficients:
##           Estimate Std. Error t value Pr(>|t|)    
## ypiYpi1  -0.002798   0.001894  -1.477 0.141882    
## ypiYpi2  -0.264323   0.058685  -4.504 1.34e-05 ***
## ypiYpi3  -0.012019   0.012843  -0.936 0.350867    
## ypiYpi4  -0.046564   0.012331  -3.776 0.000230 ***
## ypiYpi5  -0.015206   0.014875  -1.022 0.308335    
## ypiYpi6  -0.029070   0.014823  -1.961 0.051748 .  
## ypiYpi7  -0.088498   0.038429  -2.303 0.022678 *  
## ypiYpi8  -0.161353   0.038510  -4.190 4.78e-05 ***
## ypiYpi9  -0.034129   0.024395  -1.399 0.163898    
## ypiYpi10  0.020527   0.024569   0.836 0.404777    
## ypiYpi11 -0.200188   0.053366  -3.751 0.000252 ***
## ypiYpi12  0.077353   0.057336   1.349 0.179357    
## xreg     13.045368   6.105242   2.137 0.034262 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 11.55 on 148 degrees of freedom
## Multiple R-squared:  0.6451, Adjusted R-squared:  0.6139 
## F-statistic: 20.69 on 13 and 148 DF,  p-value: < 2.2e-16
## 
## Test statistic
## --------------
##         statistic p-value    
## t_1       -1.4767  0.4971    
## t_2       -4.5041       0 ***
## F_3:4      7.6734   6e-04 ***
## F_5:6      2.4761  0.0889 .  
## F_7:8     12.2074       0 ***
## F_9:10      1.335  0.2712    
## F_11:12    8.1889   4e-04 ***
## F_2:12    15.4311       0 ***
## F_1:12    14.6868       0 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1 
## 
## Deterministic terms: constant 
## Lag selection criterion and order: fixed, 0
## P-values: based on response surface regressions

Análisis de estacionariedad en la diferencia estacional

Prueba de Dickey-Fuller aumentada

adf<-ADF_test(IGAE_s, alpha=0.05)

## Warning in adf.test(x): p-value smaller than printed p-value

print(adf)

##               Estadístico Rezagos Probabilidad    Decisión      Conclusión
## Dickey-Fuller   -4.206491       5         0.01 Rechazar H0 Estacionariedad

Pueba URCA para de Dickey-Fuller para estimacion de estadisticos

ADF <- ur.df(y=IGAE_s,
             type='none',
             lags=5,
             selectlags = 'AIC')

summary(ADF)

## 
## ############################################### 
## # Augmented Dickey-Fuller Test Unit Root Test # 
## ############################################### 
## 
## Test regression none 
## 
## 
## Call:
## lm(formula = z.diff ~ z.lag.1 - 1 + z.diff.lag)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -79.063  -0.666   1.808   4.141  51.994 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## z.lag.1     -0.16929    0.04556  -3.716 0.000284 ***
## z.diff.lag1  0.28796    0.07764   3.709 0.000291 ***
## z.diff.lag2 -0.12813    0.08052  -1.591 0.113612    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.662 on 152 degrees of freedom
## Multiple R-squared:  0.1604, Adjusted R-squared:  0.1439 
## F-statistic: 9.682 on 3 and 152 DF,  p-value: 6.912e-06
## 
## 
## Value of test-statistic is: -3.7157 
## 
## Critical values for test statistics: 
##       1pct  5pct 10pct
## tau1 -2.58 -1.95 -1.62

Interpretacion

Podemos evidenciar que el estadistico de urca de la prueba de Dickey-Fuller esta por debajo del valor critico que este caso trabajaresmos al 5%, por lo que rechazamos la hipotesis nula de no estacionariedad y concluimos que el modelo es estacionario.

Prueba de Phillips-Perron

PP <-PP_test(IGAE_s)

## Warning in pp.test(x): p-value smaller than printed p-value

print(PP)

##                        Estadístico Rezagos Probabilidad    Decisión
## Dickey-Fuller Z(alpha)   -34.83503       4         0.01 Rechazar H0
##                             Conclusión
## Dickey-Fuller Z(alpha) Estacionariedad

Pueba URCA para de Phillips-Perron para estimacion de estadisticos

PP <- ur.pp(x=IGAE_s,
            type='Z-tau',
            model='constant',
            lags='short')

summary(PP)

## 
## ################################## 
## # Phillips-Perron Unit Root Test # 
## ################################## 
## 
## Test regression with intercept 
## 
## 
## Call:
## lm(formula = y ~ y.l1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -84.050  -1.605   0.353   2.237  57.730 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  1.75178    0.87991   1.991   0.0482 *  
## y.l1         0.80439    0.04726  17.020   <2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 9.826 on 158 degrees of freedom
## Multiple R-squared:  0.6471, Adjusted R-squared:  0.6448 
## F-statistic: 289.7 on 1 and 158 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic, type: Z-tau  is: -4.2979 
## 
##          aux. Z statistics
## Z-tau-mu             2.064
## 
## Critical values for Z statistics: 
##                      1pct      5pct     10pct
## critical values -3.472136 -2.879539 -2.576262

Interpretacion

Podemos evidenciar que el estadistico de urca de la prueba de Phillips-Perron esta por debajo del valor critico que este caso trabajaresmos al 5%, por lo que rechazamos la hipotesis nula de no estacionariedad y concluimos que el modelo es estacionario.

Prueba KPSS

KPSS<-KPSS_test(IGAE_s, type='Level')

## Warning in kpss.test(x, null = type): p-value greater than printed p-value

print(KPSS)

##            Estadístico Rezagos Probabilidad       Decisión      Conclusión
## KPSS Level   0.2341268       4          0.1 No rechazar H0 Estacionariedad

Pueba URCA para de KPSS para estimacion de estadisticos

KPSS <- ur.kpss(y=IGAE_s,
                type='tau',
                lags='short')
summary(KPSS)

## 
## ####################### 
## # KPSS Unit Root Test # 
## ####################### 
## 
## Test is of type: tau with 4 lags. 
## 
## Value of test-statistic is: 0.1044 
## 
## Critical value for a significance level of: 
##                 10pct  5pct 2.5pct  1pct
## critical values 0.119 0.146  0.176 0.216

Interpretacion

Podemos evidenciar que el estadistico de urca de la prueba de KPSS esta por debajo del valor critico que este caso trabajaresmos al 5%, por lo que rechazamos la hipotesis nula de no estacionariedad y concluimos que el modelo es estacionario.

Funciones ACF y PACF

ACF_s<-Acf(IGAE_s, plot=F) %>% 
  autoplot()+
  labs(title = 'Función de autocorrelación simple',
       subtitle = 'Serie: IGAE en diferencia estacional \n Selección de orden Q',
       x='Rezagos',
       y=expression(rho[t]))

PACF_s<-Pacf(IGAE_s, plot=F) %>% 
  autoplot()+
  labs(title = 'Función de autocorrelación parcial',
       subtitle = 'Serie: IGAE en diferencia estacional \n Selección de orden P',
       x='Rezagos',
       y=expression(rho[t]))

library(gridExtra)

## 
## Attaching package: 'gridExtra'

## The following object is masked from 'package:dplyr':
## 
##     combine

grid.arrange(ACF_s, PACF_s, ncol=2)

Selecion de Ordenes P, D=0, Q

ARIMA(2,1,1) ARIMA(2,1,2) ARIMA(2,1,3) ARIMA(2,1,4) ARIMA(2,1,5) ARIMA(4,1,1) ARIMA(4,1,2) ARIMA(4,1,3) ARIMA(4,1,4) ARIMA(4,1,5)

P = 1 y 2 ; Q = 1-3

Estimación del modelo SARIMA

library(stats)

# SARIMA(2,1,1)(1,0,1)[4]
mod1 <- stats::arima(IGAE_ts, order=c(1,1,1),
                     seasonal=list(order=c(2,0,1), period=12))
# SARIMA(2,1,1)(1,0,2)[4]
mod2 <- stats::arima(IGAE_ts, order=c(2,1,1),
                     seasonal=list(order=c(1,0,2), period=12))
# SARIMA(2,1,1)(1,0,3)[4]
mod3 <- stats::arima(IGAE_ts, order=c(2,1,1),
                     seasonal=list(order=c(1,0,3), period=12))
# SARIMA(2,1,1)(2,0,1)[4]
mod4 <- stats::arima(IGAE_ts, order=c(2,1,1),
                     seasonal=list(order=c(2,0,1), period=12))
# SARIMA(2,1,1)(2,0,2)[4]
mod5 <- stats::arima(IGAE_ts, order=c(2,1,1),
                     seasonal=list(order=c(2,0,2), period=12))
# SARIMA(2,1,1)(2,0,3)[4]
mod6 <- stats::arima(IGAE_ts, order=c(2,1,1),
                     seasonal=list(order=c(2,0,3), period=12))
# SARIMA(2,1,2)(1,0,1)[4]
mod7 <- stats::arima(IGAE_ts, order=c(2,1,2),
                     seasonal=list(order=c(1,0,1), period=12))
# SARIMA(2,1,2)(1,0,2)[4]
mod8 <- stats::arima(IGAE_ts, order=c(2,1,2),
                     seasonal=list(order=c(1,0,2), period=12))
# SARIMA(2,1,2)(1,0,3)[4]
mod9 <- stats::arima(IGAE_ts, order=c(2,1,2),
                     seasonal=list(order=c(1,0,3), period=12))
# SARIMA(2,1,2)(2,0,1)[4]
mod10 <- stats::arima(IGAE_ts, order=c(2,1,2),
                     seasonal=list(order=c(2,0,1), period=12))
# SARIMA(2,1,2)(2,0,2)[4]
mod11 <- stats::arima(IGAE_ts, order=c(2,1,2),
                     seasonal=list(order=c(2,0,2), period=12))
# SARIMA(2,1,2)(2,0,3)[4]
mod12 <- stats::arima(IGAE_ts, order=c(2,1,2),
                     seasonal=list(order=c(2,0,3), period=12))
# SARIMA(2,1,3)(1,0,1)[4]
mod13 <- stats::arima(IGAE_ts, order=c(2,1,3),
                     seasonal=list(order=c(1,0,1), period=12))
# SARIMA(2,1,3)(1,0,2)[4]
mod14 <- stats::arima(IGAE_ts, order=c(2,1,3),
                     seasonal=list(order=c(1,0,2), period=12))
# SARIMA(2,1,3)(1,0,3)[4]
mod15 <- stats::arima(IGAE_ts, order=c(2,1,3),
                     seasonal=list(order=c(1,0,3), period=12))
# SARIMA(2,1,3)(2,0,1)[4]
mod16 <- stats::arima(IGAE_ts, order=c(2,1,3),
                     seasonal=list(order=c(2,0,1), period=12))
# SARIMA(2,1,3)(2,0,2)[4]
mod17 <- stats::arima(IGAE_ts, order=c(2,1,3),
                     seasonal=list(order=c(2,0,2), period=12))
# SARIMA(2,1,3)(2,0,3)[4]
mod18 <- stats::arima(IGAE_ts, order=c(2,1,3),
                     seasonal=list(order=c(2,0,3), period=12))
# SARIMA(2,1,4)(1,0,1)[4]
mod19 <- stats::arima(IGAE_ts, order=c(2,1,4),
                     seasonal=list(order=c(1,0,1), period=12))
# SARIMA(2,1,4)(1,0,2)[4]
mod20 <- stats::arima(IGAE_ts, order=c(2,1,4),
                     seasonal=list(order=c(1,0,2), period=12))
# SARIMA(2,1,4)(1,0,3)[4]
mod21 <- stats::arima(IGAE_ts, order=c(2,1,4),
                     seasonal=list(order=c(1,0,3), period=12))
# SARIMA(2,1,4)(2,0,1)[4]
mod22 <- stats::arima(IGAE_ts, order=c(2,1,4),
                     seasonal=list(order=c(2,0,1), period=12))
# SARIMA(2,1,4)(2,0,2)[4]
mod23 <- stats::arima(IGAE_ts, order=c(2,1,4),
                     seasonal=list(order=c(2,0,2), period=12))
# SARIMA(2,1,4)(2,0,3)[4]
mod24 <- stats::arima(IGAE_ts, order=c(2,1,4),
                     seasonal=list(order=c(2,0,3), period=12))
# SARIMA(2,1,5)(1,0,1)[4]
mod25 <- stats::arima(IGAE_ts, order=c(2,1,5),
                     seasonal=list(order=c(1,0,1), period=12))
# SARIMA(2,1,5)(1,0,2)[4]
mod26 <- stats::arima(IGAE_ts, order=c(2,1,5),
                     seasonal=list(order=c(1,0,2), period=12))
# SARIMA(2,1,5)(1,0,3)[4]
mod27 <- stats::arima(IGAE_ts, order=c(2,1,5),
                     seasonal=list(order=c(1,0,3), period=12))
# SARIMA(2,1,5)(2,0,1)[4]
mod28 <- stats::arima(IGAE_ts, order=c(2,1,5),
                     seasonal=list(order=c(2,0,1), period=12))
# SARIMA(2,1,5)(2,0,2)[4]
mod29 <- stats::arima(IGAE_ts, order=c(2,1,5),
                     seasonal=list(order=c(2,0,2), period=12))
# SARIMA(2,1,5)(2,0,3)[4]
mod30 <- stats::arima(IGAE_ts, order=c(2,1,5),
                     seasonal=list(order=c(2,0,3), period=12))
# SARIMA(4,1,1)(1,0,1)[4]
mod31 <- stats::arima(IGAE_ts, order=c(4,1,1),
                     seasonal=list(order=c(1,0,1), period=12))
# SARIMA(4,1,1)(1,0,2)[4]
mod32 <- stats::arima(IGAE_ts, order=c(4,1,1),
                     seasonal=list(order=c(1,0,2), period=12))
# SARIMA(4,1,1)(1,0,3)[4]
mod33 <- stats::arima(IGAE_ts, order=c(4,1,1),
                     seasonal=list(order=c(1,0,3), period=12))
# SARIMA(4,1,1)(2,0,1)[4]
mod34 <- stats::arima(IGAE_ts, order=c(4,1,1),
                     seasonal=list(order=c(2,0,1), period=12))
# SARIMA(4,1,1)(2,0,2)[4]
mod35 <- stats::arima(IGAE_ts, order=c(4,1,1),
                     seasonal=list(order=c(2,0,2), period=12))
# SARIMA(4,1,1)(2,0,3)[4]
mod36 <- stats::arima(IGAE_ts, order=c(4,1,1),
                     seasonal=list(order=c(2,0,3), period=12))
# SARIMA(4,1,2)(1,0,1)[4]
mod37 <- stats::arima(IGAE_ts, order=c(4,1,2),
                     seasonal=list(order=c(1,0,1), period=12))
# SARIMA(4,1,2)(1,0,2)[4]
mod38 <- stats::arima(IGAE_ts, order=c(4,1,2),
                     seasonal=list(order=c(1,0,2), period=12))
# SARIMA(4,1,2)(1,0,3)[4]
mod39 <- stats::arima(IGAE_ts, order=c(4,1,2),
                     seasonal=list(order=c(1,0,3), period=12))
# SARIMA(4,1,2)(2,0,1)[4]
mod40 <- stats::arima(IGAE_ts, order=c(4,1,2),
                     seasonal=list(order=c(2,0,1), period=12))
# SARIMA(4,1,2)(2,0,2)[4]
mod41 <- stats::arima(IGAE_ts, order=c(4,1,2),
                     seasonal=list(order=c(2,0,2), period=12))
# SARIMA(4,1,2)(2,0,3)[4]
mod42 <- stats::arima(IGAE_ts, order=c(4,1,2),
                     seasonal=list(order=c(2,0,3), period=12))
# SARIMA(4,1,3)(1,0,1)[4]
mod43 <- stats::arima(IGAE_ts, order=c(4,1,3),
                     seasonal=list(order=c(1,0,1), period=12))
# SARIMA(4,1,3)(1,0,2)[4]
mod44 <- stats::arima(IGAE_ts, order=c(4,1,3),
                     seasonal=list(order=c(1,0,2), period=12))
# SARIMA(4,1,3)(1,0,3)[4]
mod45 <- stats::arima(IGAE_ts, order=c(4,1,3),
                     seasonal=list(order=c(1,0,3), period=12))
# SARIMA(4,1,3)(2,0,1)[4]
mod46 <- stats::arima(IGAE_ts, order=c(4,1,3),
                     seasonal=list(order=c(2,0,1), period=12))
# SARIMA(4,1,3)(2,0,2)[4]
mod47 <- stats::arima(IGAE_ts, order=c(4,1,3),
                     seasonal=list(order=c(2,0,2), period=12))
# SARIMA(4,1,3)(2,0,3)[4]
mod48 <- stats::arima(IGAE_ts, order=c(4,1,3),
                     seasonal=list(order=c(2,0,3), period=12))
# SARIMA(4,1,4)(1,0,1)[4]
mod49 <- stats::arima(IGAE_ts, order=c(4,1,4),
                     seasonal=list(order=c(1,0,1), period=12))
# SARIMA(4,1,4)(1,0,2)[4]
mod50 <- stats::arima(IGAE_ts, order=c(4,1,4),
                     seasonal=list(order=c(1,0,2), period=12))
# SARIMA(4,1,4)(1,0,3)[4]
mod51 <- stats::arima(IGAE_ts, order=c(4,1,4),
                     seasonal=list(order=c(1,0,3), period=12))
# SARIMA(4,1,4)(2,0,1)[4]
mod52 <- stats::arima(IGAE_ts, order=c(4,1,4),
                     seasonal=list(order=c(2,0,1), period=12))
# SARIMA(4,1,4)(2,0,2)[4]
mod53 <- stats::arima(IGAE_ts, order=c(4,1,4),
                     seasonal=list(order=c(2,0,2), period=12))
# SARIMA(4,1,4)(2,0,3)[4]
mod54 <- stats::arima(IGAE_ts, order=c(4,1,4),
                     seasonal=list(order=c(2,0,3), period=12))
# SARIMA(4,1,5)(1,0,1)[4]
mod55 <- stats::arima(IGAE_ts, order=c(4,1,5),
                     seasonal=list(order=c(1,0,1), period=12))
# SARIMA(4,1,5)(1,0,2)[4]
mod56 <- stats::arima(IGAE_ts, order=c(4,1,5),
                     seasonal=list(order=c(1,0,2), period=12))
# SARIMA(4,1,5)(1,0,3)[4]
mod57 <- stats::arima(IGAE_ts, order=c(4,1,5),
                     seasonal=list(order=c(1,0,3), period=12))
# SARIMA(4,1,5)(2,0,1)[4]
mod58 <- stats::arima(IGAE_ts, order=c(4,1,5),
                     seasonal=list(order=c(2,0,1), period=12))
# SARIMA(4,1,5)(2,0,2)[4]
mod59 <- stats::arima(IGAE_ts, order=c(4,1,5),
                     seasonal=list(order=c(2,0,2), period=12))
# SARIMA(4,1,5)(2,0,3)[4]
mod60 <- stats::arima(IGAE_ts, order=c(4,1,5),
                     seasonal=list(order=c(2,0,3), period=12))

Modelo AUTO

auto <- auto.arima(IGAE_ts, seasonal = TRUE,
                   trace = TRUE)

## 
##  Fitting models using approximations to speed things up...
## 
##  ARIMA(2,0,2)(1,1,1)[12] with drift         : 1092.324
##  ARIMA(0,0,0)(0,1,0)[12] with drift         : 1291.317
##  ARIMA(1,0,0)(1,1,0)[12] with drift         : 1108.029
##  ARIMA(0,0,1)(0,1,1)[12] with drift         : 1143.964
##  ARIMA(0,0,0)(0,1,0)[12]                    : 1329.767
##  ARIMA(2,0,2)(0,1,1)[12] with drift         : 1078.62
##  ARIMA(2,0,2)(0,1,0)[12] with drift         : 1117.004
##  ARIMA(2,0,2)(0,1,2)[12] with drift         : 1080.827
##  ARIMA(2,0,2)(1,1,0)[12] with drift         : 1103.408
##  ARIMA(2,0,2)(1,1,2)[12] with drift         : 1094.473
##  ARIMA(1,0,2)(0,1,1)[12] with drift         : 1075.932
##  ARIMA(1,0,2)(0,1,0)[12] with drift         : 1113.868
##  ARIMA(1,0,2)(1,1,1)[12] with drift         : 1089.655
##  ARIMA(1,0,2)(0,1,2)[12] with drift         : 1078.118
##  ARIMA(1,0,2)(1,1,0)[12] with drift         : 1101.16
##  ARIMA(1,0,2)(1,1,2)[12] with drift         : Inf
##  ARIMA(0,0,2)(0,1,1)[12] with drift         : 1107.424
##  ARIMA(1,0,1)(0,1,1)[12] with drift         : 1074.167
##  ARIMA(1,0,1)(0,1,0)[12] with drift         : 1112.064
##  ARIMA(1,0,1)(1,1,1)[12] with drift         : 1087.843
##  ARIMA(1,0,1)(0,1,2)[12] with drift         : 1076.319
##  ARIMA(1,0,1)(1,1,0)[12] with drift         : 1099.002
##  ARIMA(1,0,1)(1,1,2)[12] with drift         : Inf
##  ARIMA(1,0,0)(0,1,1)[12] with drift         : 1081.562
##  ARIMA(2,0,1)(0,1,1)[12] with drift         : 1076.836
##  ARIMA(0,0,0)(0,1,1)[12] with drift         : 1273.702
##  ARIMA(2,0,0)(0,1,1)[12] with drift         : 1076.915
##  ARIMA(1,0,1)(0,1,1)[12]                    : 1083.964
## 
##  Now re-fitting the best model(s) without approximations...
## 
##  ARIMA(1,0,1)(0,1,1)[12] with drift         : 1148.29
## 
##  Best model: ARIMA(1,0,1)(0,1,1)[12] with drift

Evaluación del modeo

Revisión de la ecuación

summary(auto)

## Series: IGAE_ts 
## ARIMA(1,0,1)(0,1,1)[12] with drift 
## 
## Coefficients:
##          ar1     ma1     sma1   drift
##       0.7466  0.3148  -0.5197  0.7589
## s.e.  0.0620  0.0932   0.0687  0.1427
## 
## sigma^2 = 68.27:  log likelihood = -568.95
## AIC=1147.9   AICc=1148.29   BIC=1163.31
## 
## Training set error measures:
##                     ME    RMSE      MAE       MPE     MAPE      MASE
## Training set 0.0146836 7.87111 3.780458 -0.048506 1.396329 0.2752209
##                      ACF1
## Training set -0.002557715

Revisar la significancia del modelo

library(lmtest)

## Warning: package 'lmtest' was built under R version 4.2.1

coeftest(auto)

## 
## z test of coefficients:
## 
##        Estimate Std. Error z value  Pr(>|z|)    
## ar1    0.746639   0.062001 12.0424 < 2.2e-16 ***
## ma1    0.314834   0.093221  3.3773  0.000732 ***
## sma1  -0.519653   0.068679 -7.5664 3.837e-14 ***
## drift  0.758870   0.142689  5.3184 1.047e-07 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Raíces invertidas

autoplot(auto)

Pruebas de autocorrelación

res <- auto$residuals

Función de autocorrelación simple

Acf(res, plot=F) %>% 
  autoplot()+
  labs(title = 'Función de autocorrelación simple',
       subtitle = paste0("Residuales del modelo auto: SARIMA(1,0,1)(0,1,1)[12]"),
       x='Rezagos',
       y=expression(rho[t]))

Prueba de Ljung-Box

# 4 lags son suficientes para validar el modelo

Box.test(res, lag=4, type='Ljung-Box')

## 
##  Box-Ljung test
## 
## data:  res
## X-squared = 0.19141, df = 4, p-value = 0.9957

Box.test(res, lag=4, type='Box-Pierce')

## 
##  Box-Pierce test
## 
## data:  res
## X-squared = 0.18548, df = 4, p-value = 0.996

Homocedasticidad

ACF con residuales al cuadrado

res2<-res^2
Acf(res2, plot=F) %>% 
  autoplot()+
  labs(title = 'Función de autocorrelación simple',
       subtitle = paste0("Residuales al cuadrado del modelo auto: SARIMA(1,0,1)(0,1,1)[12]"),
       x='Rezagos',
       y=expression(rho[t]))

Box.test(res2, lag=4, type='Ljung-Box')

## 
##  Box-Ljung test
## 
## data:  res2
## X-squared = 1.5104, df = 4, p-value = 0.8248

Box.test(res2, lag=4, type='Box-Pierce')

## 
##  Box-Pierce test
## 
## data:  res2
## X-squared = 1.4601, df = 4, p-value = 0.8337

Normalidad

jarque.bera.test(res)

## 
##  Jarque Bera Test
## 
## data:  res
## X-squared = 30348, df = 2, p-value < 2.2e-16

shapiro.test(res)

## 
##  Shapiro-Wilk normality test
## 
## data:  res
## W = 0.56939, p-value < 2.2e-16

Estacionariedad de los residuales

detach(unload = TRUE)
ADF_test(res)

##               Estadístico Rezagos Probabilidad    Decisión      Conclusión
## Dickey-Fuller   -5.062112       5         0.01 Rechazar H0 Estacionariedad

PP_test(res)

##                        Estadístico Rezagos Probabilidad    Decisión
## Dickey-Fuller Z(alpha)   -169.6463       4         0.01 Rechazar H0
##                             Conclusión
## Dickey-Fuller Z(alpha) Estacionariedad

KPSS_test(res, type='Level')

##            Estadístico Rezagos Probabilidad       Decisión      Conclusión
## KPSS Level   0.2742535       4          0.1 No rechazar H0 Estacionariedad

Pronostico

fcast <- forecast(auto, h=7, level=95)
fcast

##          Point Forecast    Lo 95    Hi 95
## Jun 2022       316.6424 300.4483 332.8366
## Jul 2022       312.3059 288.6896 335.9223
## Aug 2022       314.7254 287.8469 341.6040
## Sep 2022       340.0303 311.4947 368.5660
## Oct 2022       361.4960 332.0771 390.9150
## Nov 2022       347.2530 317.3530 377.1531
## Dec 2022       355.9941 325.8292 386.1590

autoplot(fcast, include = 50)

Fechas Pronosticadas

Fecha = (“2022-06-01”, “2022-07-01”, “2022-08-01”, “2022-09-01”, “2022-10-01”, “2022-11-01”, “2022-12-01”)