Econometría II: ARIMA

Departamento de Economía

Carlos A. Yanes G.

2024-03-19

Paquetes con que se trabaja la sesión

R libraries

Los paquetes que vamos a utilizar y que se deben instalar para ser usados en clase son:

install.packages(c("pacman","TSstudio", "urca", "forecast", "devtools"))

# Si tiene problemas con TSstudio podemos mirar entonces en Github
devtools::install_github("RamiKrispin/TSstudio")

Ejecución

No olvide que debe desde luego tener presente en su documento de R Markdown cargar estos paquetes.

library(pacman)
p_load(tidyverse, TSstudio, tseries, urca, forecast)

Preambulo

Modelar ARIMA

Debemos primero mirar la estacionariedad.
Si no es estacionaria la serie hay que diferenciar.
Después de testear procedemos a leer sus funciones de autocovarianzas denominadas correlogramas
Constatar el RUIDO BLANCO
Predecir

Datos

Cartera comercial de bancos

Enero de 2016 hasta enero de 2024
Web site de Superintendencia Financiera
Frecuencia Mensual

Agradecimientos

Esta sesión se construye a partir de los elementos de los paquetes de tseries y (Krispin 2019), quien desarrolló el paquete TSstudio respectivamente

A lo que vinimos

Cargar datos

library(pacman)
p_load(readxl, TSstudio, tidyverse, stats, urca, forecast, ggfortify, ggplot2, tseries)
bd <- read_excel("cartera.xlsx")
bd <- bd |> select(Cartera)
cartera <- ts(bd, frequency=12, start=c(2016,1))

Graficar normal

Code

ts_plot(cartera,
        title = "Cartera comercial de Bancos", 
        Ytitle = "En Miles de Millones de $",
        Xtitle = "Años")

Code

prim <- ur.df(y=cartera, lags=0, type='none')
summary(prim)


############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression none 


Call:
lm(formula = z.diff ~ z.lag.1 - 1)

Residuals:
    Min      1Q  Median      3Q     Max 
-8178.3 -2346.7   -61.7  1996.2 15252.0 

Coefficients:
         Estimate Std. Error t value Pr(>|t|)    
z.lag.1 0.0059658  0.0007548   7.904 4.79e-12 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3662 on 95 degrees of freedom
Multiple R-squared:  0.3967,    Adjusted R-squared:  0.3904 
F-statistic: 62.48 on 1 and 95 DF,  p-value: 4.79e-12


Value of test-statistic is: 7.9042 

Critical values for test statistics: 
     1pct  5pct 10pct
tau1 -2.6 -1.95 -1.61

Code

prim2 <- ur.df(y=cartera, lags=0, type='drift')
summary(prim2)


############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression drift 


Call:
lm(formula = z.diff ~ z.lag.1 + 1)

Residuals:
    Min      1Q  Median      3Q     Max 
-8177.7 -2364.0   -51.6  2005.6 15255.0 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept) -6.100e+01  2.117e+03  -0.029    0.977
z.lag.1      6.087e-03  4.275e-03   1.424    0.158

Residual standard error: 3681 on 94 degrees of freedom
Multiple R-squared:  0.02111,   Adjusted R-squared:  0.0107 
F-statistic: 2.028 on 1 and 94 DF,  p-value: 0.1578


Value of test-statistic is: 1.4239 30.9097 

Critical values for test statistics: 
      1pct  5pct 10pct
tau2 -3.51 -2.89 -2.58
phi1  6.70  4.71  3.86

Code

prim3 <- ur.df(y=cartera, lags=0, type='trend')
summary(prim3)


############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression trend 


Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt)

Residuals:
    Min      1Q  Median      3Q     Max 
-8473.6 -1890.1    75.1  1723.4 14870.7 

Coefficients:
              Estimate Std. Error t value Pr(>|t|)
(Intercept) 8918.21856 5878.88458   1.517    0.133
z.lag.1       -0.02124    0.01724  -1.232    0.221
tt            89.40563   54.68086   1.635    0.105

Residual standard error: 3649 on 93 degrees of freedom
Multiple R-squared:  0.04847,   Adjusted R-squared:  0.028 
F-statistic: 2.368 on 2 and 93 DF,  p-value: 0.09925


Value of test-statistic is: -1.2318 21.8644 2.3685 

Critical values for test statistics: 
      1pct  5pct 10pct
tau3 -4.04 -3.45 -3.15
phi2  6.50  4.88  4.16
phi3  8.73  6.49  5.47

Code

# adf.test(cartera,k=0)
adf.test(cartera) # automático


    Augmented Dickey-Fuller Test

data:  cartera
Dickey-Fuller = -2.4087, Lag order = 4, p-value = 0.4075
alternative hypothesis: stationary

Analisis de funciones de autocorrelación

Anotaciones
ACF
PAC

Permiten mirar o -sospechar- si la serie es estacionaria
La serie de cartera por ninguna parte parece ser estacionaria

Code

acf(cartera, lag.max = 20)

Code

pacf(cartera, lag.max = 20)

Comportamiento Lags

Code

ts_lags(cartera)

Comportamiento Lags

Code

ts_lags(cartera, lags=c(1,4,6,11))

Diferenciemos

Code

ts_plot(diff(cartera, lag = 1),
        title = "Transformación en Nivel",
        Xtitle = "Años",
        Ytitle = "Diferencia de cartera")

Diferenciación opcional

Code

ts_plot(diff(log(cartera), lag = 1),
        title = "Transformación en logaritmo",
        Xtitle = "Años",
        Ytitle = "Diferencia/logaritmica de cartera")

Nuevamente test e identificación

DF+
DF C+T
Recuerde

Code

dcartera<-diff(cartera, lag=1) # creo el objeto
adf.test(dcartera)


    Augmented Dickey-Fuller Test

data:  dcartera
Dickey-Fuller = -2.6412, Lag order = 4, p-value = 0.3115
alternative hypothesis: stationary

Code

prim4 <- ur.df(dcartera, lags=0, type='trend')
summary(prim4)


############################################### 
# Augmented Dickey-Fuller Test Unit Root Test # 
############################################### 

Test regression trend 


Call:
lm(formula = z.diff ~ z.lag.1 + 1 + tt)

Residuals:
    Min      1Q  Median      3Q     Max 
-7972.0 -1896.3    30.4  1817.9 13792.0 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 943.03465  701.55496   1.344    0.182    
z.lag.1      -0.55233    0.09511  -5.807 8.99e-08 ***
tt           12.96238   12.68604   1.022    0.310    
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 3312 on 92 degrees of freedom
Multiple R-squared:  0.2685,    Adjusted R-squared:  0.2526 
F-statistic: 16.89 on 2 and 92 DF,  p-value: 5.66e-07


Value of test-statistic is: -5.8072 11.2731 16.8877 

Critical values for test statistics: 
      1pct  5pct 10pct
tau3 -4.04 -3.45 -3.15
phi2  6.50  4.88  4.16
phi3  8.73  6.49  5.47

Identificación para pronostico

Code

par(mfrow=c(1,2)) 
acf(dcartera)
pacf(dcartera)

Modelos

Modelo 1
Modelo 2
Modelo 3

Code

arima_1 <- arima(cartera, order = c(1,1,0))
arima_1


Call:
arima(x = cartera, order = c(1, 1, 0))

Coefficients:
         ar1
      0.6682
s.e.  0.0750

sigma^2 estimated as 12084388:  log likelihood = -919.27,  aic = 1842.54

Code

arima_2 <- arima(cartera, order = c(2,1,0))
arima_2


Call:
arima(x = cartera, order = c(2, 1, 0))

Coefficients:
         ar1     ar2
      0.4917  0.2590
s.e.  0.0985  0.0984

sigma^2 estimated as 11259283:  log likelihood = -915.94,  aic = 1837.88

Code

arima_3 <- arima(cartera, order = c(3,1,0))
arima_3


Call:
arima(x = cartera, order = c(3, 1, 0))

Coefficients:
         ar1     ar2     ar3
      0.4342  0.1483  0.2187
s.e.  0.0997  0.1085  0.0993

sigma^2 estimated as 10704947:  log likelihood = -913.58,  aic = 1835.16

Tenemos modelo!!

Test de Ruido Blanco

Opción Modelo Arima (3 1 0)

Code

forecast::checkresiduals(arima(cartera, order = c(3, 1, 0), include.mean = FALSE))


    Ljung-Box test

data:  Residuals from ARIMA(3,1,0)
Q* = 33.151, df = 16, p-value = 0.007056

Model df: 3.   Total lags used: 19

Test de Ruido Blanco

Opción Modelo Arima (1 1 0)

Code

forecast::checkresiduals(arima(cartera, order = c(1, 1, 0), include.mean = FALSE))


    Ljung-Box test

data:  Residuals from ARIMA(1,1,0)
Q* = 42.759, df = 18, p-value = 0.0008655

Model df: 1.   Total lags used: 19

Test de Ruido Blanco

Lectura

La prueba de correlación serial de nuestros modelos no esta pasando el test. Se requiere que esta siempre no rechace la hipotesis nula.
Este paso es vital. Si nuestro modelo no aprueba el ruido blanco vamos a tener problemas en el futuro.
Debe existir un modelo que nos ayude a establecer que realmente es ruido blanco

Pronostico de cartera

Code

arima_fc <- forecast(arima_3, h = 12)
plot_forecast(arima_fc,
 title = "Pronostico modelo AR(3)", 
 Ytitle = "Valores",
 Xtitle = "Años")

Detras de la serie

Code

# Crear el gráfico de las dos series de tiempo
plot.ts(cartera, col = "blue", xlab = "Tiempo", ylab = "Cartera Comercial", main = "Predicción de cartera comercial")
lines(arima_fc$fitted, col = "red")
legend("topleft", legend = c("Serie Original", "Predicción"), col = c("blue", "red"), lty = 1)

Valores

Code

arima_fc

         Point Forecast    Lo 80    Hi 80    Lo 95    Hi 95
Feb 2024       645793.1 641600.1 649986.1 639380.4 652205.8
Mar 2024       644824.7 637493.7 652155.8 633612.9 656036.6
Apr 2024       643703.1 633268.3 654138.0 627744.4 659661.9
May 2024       642869.4 628940.9 656797.9 621567.6 664171.2
Jun 2024       642129.3 624616.9 659641.7 615346.4 668912.2
Jul 2024       641439.0 620325.4 662552.6 609148.5 673729.5
Aug 2024       640847.2 616101.6 665592.8 603002.1 678692.3
Sep 2024       640326.0 611959.8 668692.2 596943.7 683708.3
Oct 2024       639861.0 607908.8 671813.1 590994.4 688727.6
Nov 2024       639452.4 603957.6 674947.1 585167.9 693736.8
Dec 2024       639092.0 600109.0 678075.0 579472.7 698711.3
Jan 2025       638773.2 596363.7 681182.7 573913.5 703632.9

Selección de modelos

Existen varias formas de buscar elegir un modelo entre varios que se han estimado por pensar que es nuestro mejor modelo generador de datos. Estos vienen a ser los criterios de Akaike (AIC) y Bayes (BIC)

\(AIC=-2\times log(L)+ 2\times p\)
P es el número de parámetros del modelo
L la razón de máxima verosimilitud del modelo

\(BIC=-2\times log(L)+ p \times log(n)\)
n es es tamaño de la muestra

Gracias por su atención!!

Referencias

Krispin, Rami. 2019. Hands-on Time Series Analysis with r: Perform Time Series Analysis and Forecasting Using r. Packt Publishing Ltd.