Taller Análisis de Series de tiempo multivariadas

El comportamiento económico social es un factor importante que permite medir la capacidad de endeudamiento existente en las personas. Normalmente, cuando una persona se encuentra en un estado de desempleo, hay probabilidades mas altas de un incumplimiento en sus obligaciones con entidades financieras. A continuación se hará un análisis donde se evalúa una serie de tiempo que posee información relacionada al nivel de desempleo mensual desde el mes de enero del 2008, hasta febrero del 2019. De la misma forma, se cuenta con información del nivel de morosidad durante el mismo tiempo.

Procedimiento del análisis

#Carga de librerías

library(tseries)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
library(vars)
## Warning: package 'vars' was built under R version 4.2.3
## Loading required package: MASS
## Warning: package 'MASS' was built under R version 4.2.3
## Loading required package: strucchange
## Warning: package 'strucchange' was built under R version 4.2.3
## Loading required package: zoo
## Warning: package 'zoo' was built under R version 4.2.3
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: sandwich
## Warning: package 'sandwich' was built under R version 4.2.3
## Loading required package: urca
## Loading required package: lmtest
library(highcharter)
## Highcharts (www.highcharts.com) is a Highsoft software product which is
## not free for commercial and Governmental use
## 
## Attaching package: 'highcharter'
## The following object is masked from 'package:lmtest':
## 
##     unemployment
library(dplyr)
## 
## Attaching package: 'dplyr'
## The following object is masked from 'package:MASS':
## 
##     select
## The following objects are masked from 'package:stats':
## 
##     filter, lag
## The following objects are masked from 'package:base':
## 
##     intersect, setdiff, setequal, union
library(ggfortify)
## Warning: package 'ggfortify' was built under R version 4.2.3
## Loading required package: ggplot2
## Warning: package 'ggplot2' was built under R version 4.2.3
library(plotly)
## 
## Attaching package: 'plotly'
## The following object is masked from 'package:ggplot2':
## 
##     last_plot
## The following object is masked from 'package:MASS':
## 
##     select
## The following object is masked from 'package:stats':
## 
##     filter
## The following object is masked from 'package:graphics':
## 
##     layout
library(urca)

library(readxl)

Cargamos el documento

analisis<- "series1.xlsx"
data<-read_xlsx(analisis)
morosidad <- ts(data$Morosidad, start = c(2008,1), end = c(2019,2), frequency = 12)
desempleo <- ts(data$Desempleo2, start = c(2008,1), end = c(2019,2), frequency = 12)
plot(cbind(morosidad, desempleo), main="MOROSIDAD y DESEMPLEO" )

hchart(cbind(morosidad, desempleo))
## Warning: Deprecated function. Use the `create_axis` function.

#Pruebas de raiz unitaria

adf.mor <- ur.df(morosidad, type = "trend", selectlags = "BIC")
summary(adf.mor)   
## 
## ############################################### 
## # 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 
## -0.0054254 -0.0010795 -0.0002387  0.0010075  0.0088360 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)  
## (Intercept)  5.651e-04  1.006e-03   0.562   0.5753  
## z.lag.1     -7.463e-03  1.223e-02  -0.610   0.5426  
## tt           3.131e-06  4.397e-06   0.712   0.4777  
## z.diff.lag   1.815e-01  8.732e-02   2.079   0.0397 *
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.001913 on 128 degrees of freedom
## Multiple R-squared:  0.03837,    Adjusted R-squared:  0.01583 
## F-statistic: 1.702 on 3 and 128 DF,  p-value: 0.1698
## 
## 
## Value of test-statistic is: -0.6105 0.6593 0.4047 
## 
## 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

Se observa que la serie de tiempo de morosidad no es estacionaria

adf.des <- ur.df(desempleo, type = "trend", selectlags = "BIC")
summary(adf.des) 
## 
## ############################################### 
## # 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 
## -0.0084458 -0.0025878 -0.0000842  0.0022558  0.0082096 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  8.944e-03  5.909e-03   1.514    0.133    
## z.lag.1     -7.759e-02  4.893e-02  -1.586    0.115    
## tt          -1.567e-05  1.501e-05  -1.044    0.298    
## z.diff.lag  -3.690e-01  8.431e-02  -4.377 2.48e-05 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.003521 on 128 degrees of freedom
## Multiple R-squared:  0.1868, Adjusted R-squared:  0.1677 
## F-statistic: 9.801 on 3 and 128 DF,  p-value: 7.25e-06
## 
## 
## Value of test-statistic is: -1.586 0.9597 1.4081 
## 
## 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

Por otro lado, se observa que la serie de desempleo si posee estacionareidad

Para solucionar la no estacionareidad de la serie de morosidad, se debe aplicar una diferencia:

adf.mor <- ur.df(diff(morosidad), type = "trend", selectlags = "BIC")
summary(adf.mor)  
## 
## ############################################### 
## # 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 
## -0.0052319 -0.0010810 -0.0001969  0.0008731  0.0083613 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -5.126e-05  3.311e-04  -0.155  0.87719    
## z.lag.1     -6.273e-01  1.103e-01  -5.686 8.47e-08 ***
## tt           2.668e-06  4.321e-06   0.617  0.53806    
## z.diff.lag  -2.426e-01  8.586e-02  -2.826  0.00547 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.001862 on 127 degrees of freedom
## Multiple R-squared:  0.4491, Adjusted R-squared:  0.4361 
## F-statistic: 34.51 on 3 and 127 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -5.6864 10.7853 16.1721 
## 
## 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

En este momento ya se asegura que está convertida en estacionaria.

A continuación se grafica la serie de desempleo

plot(adf.des)

y posteriormente la serie de morosidad

plot(adf.mor)

#Selección y estimación del modelo

dat.bv <- cbind(desempleo, morosidad)
colnames(dat.bv) <- c("desempleo", "morosidad")
info.bv <- VARselect(dat.bv, lag.max = 14, type = "const")
info.bv
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      3      3      3      3 
## 
## $criteria
##                    1             2             3             4             5
## AIC(n) -2.396597e+01 -2.408946e+01 -2.422351e+01 -2.422118e+01 -2.419424e+01
## HQ(n)  -2.390937e+01 -2.399512e+01 -2.409144e+01 -2.405138e+01 -2.398670e+01
## SC(n)  -2.382659e+01 -2.385717e+01 -2.389830e+01 -2.380305e+01 -2.368320e+01
## FPE(n)  3.905899e-11  3.452418e-11  3.019807e-11  3.027757e-11  3.111884e-11
##                    6             7             8             9            10
## AIC(n) -2.415728e+01 -2.410254e+01 -2.411959e+01 -2.411007e+01 -2.405412e+01
## HQ(n)  -2.391201e+01 -2.381953e+01 -2.379885e+01 -2.375160e+01 -2.365792e+01
## SC(n)  -2.355332e+01 -2.340567e+01 -2.332980e+01 -2.322736e+01 -2.307850e+01
## FPE(n)  3.231239e-11  3.416193e-11  3.362496e-11  3.399878e-11  3.602370e-11
##                   11            12            13            14
## AIC(n) -2.408038e+01 -2.403057e+01 -2.409741e+01 -2.405081e+01
## HQ(n)  -2.364644e+01 -2.355890e+01 -2.358801e+01 -2.350367e+01
## SC(n)  -2.301184e+01 -2.286911e+01 -2.284304e+01 -2.270352e+01
## FPE(n)  3.517180e-11  3.707091e-11  3.478805e-11  3.658803e-11
info.bv$selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      3      3      3      3
library(MTS)
## Warning: package 'MTS' was built under R version 4.2.3
## 
## Attaching package: 'MTS'
## The following object is masked from 'package:vars':
## 
##     VAR
VARorder(dat.bv)
## selected order: aic =  4 
## selected order: bic =  3 
## selected order: hq =  3 
## Summary table:  
##        p      AIC      BIC       HQ     M(p) p-value
##  [1,]  0 -17.5420 -17.5420 -17.5420   0.0000  0.0000
##  [2,]  1 -23.9916 -23.9051 -23.9565 764.8469  0.0000
##  [3,]  2 -24.1366 -23.9636 -24.0663  23.6361  0.0001
##  [4,]  3 -24.2546 -23.9951 -24.1491  20.1688  0.0005
##  [5,]  4 -24.2818 -23.9358 -24.1412   9.6937  0.0459
##  [6,]  5 -24.2660 -23.8335 -24.0902   4.8027  0.3082
##  [7,]  6 -24.2326 -23.7135 -24.0216   2.8258  0.5874
##  [8,]  7 -24.1830 -23.5775 -23.9370   1.0749  0.8982
##  [9,]  8 -24.2228 -23.5308 -23.9416  10.2968  0.0357
## [10,]  9 -24.2216 -23.4431 -23.9052   5.9337  0.2042
## [11,] 10 -24.1724 -23.3074 -23.8209   1.0468  0.9026
## [12,] 11 -24.1999 -23.2483 -23.8132   8.4992  0.0749
## [13,] 12 -24.1566 -23.1186 -23.7348   1.5698  0.8142
## [14,] 13 -24.2280 -23.1034 -23.7710  12.2551  0.0156

De acuerdo a los diferentes criterios, se observa que los rezagos sugeridos a utilizar son 3.

#Estimación del modelo VAR

modelo <- vars::VAR(dat.bv, p = 3, type = "const")

summary(modelo)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: desempleo, morosidad 
## Deterministic variables: const 
## Sample size: 131 
## Log Likelihood: 1204.969 
## Roots of the characteristic polynomial:
## 0.9845 0.9845 0.5097 0.4398 0.4398 0.3497
## Call:
## vars::VAR(y = dat.bv, p = 3, type = "const")
## 
## 
## Estimation results for equation desempleo: 
## ========================================== 
## desempleo = desempleo.l1 + morosidad.l1 + desempleo.l2 + morosidad.l2 + desempleo.l3 + morosidad.l3 + const 
## 
##               Estimate Std. Error t value Pr(>|t|)    
## desempleo.l1  0.460676   0.087663   5.255 6.25e-07 ***
## morosidad.l1  0.064752   0.161170   0.402  0.68855    
## desempleo.l2  0.281637   0.093604   3.009  0.00318 ** 
## morosidad.l2 -0.228378   0.245106  -0.932  0.35328    
## desempleo.l3  0.204546   0.086937   2.353  0.02021 *  
## morosidad.l3  0.213942   0.162060   1.320  0.18922    
## const         0.001372   0.002831   0.485  0.62877    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.00344 on 124 degrees of freedom
## Multiple R-Squared: 0.9198,  Adjusted R-squared: 0.9159 
## F-statistic: 237.1 on 6 and 124 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation morosidad: 
## ========================================== 
## morosidad = desempleo.l1 + morosidad.l1 + desempleo.l2 + morosidad.l2 + desempleo.l3 + morosidad.l3 + const 
## 
##               Estimate Std. Error t value Pr(>|t|)    
## desempleo.l1  0.039246   0.047329   0.829  0.40858    
## morosidad.l1  1.106029   0.087016  12.711  < 2e-16 ***
## desempleo.l2  0.011585   0.050537   0.229  0.81906    
## morosidad.l2  0.121020   0.132334   0.915  0.36223    
## desempleo.l3 -0.067720   0.046937  -1.443  0.15160    
## morosidad.l3 -0.234190   0.087497  -2.677  0.00844 ** 
## const         0.002442   0.001528   1.598  0.11259    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.001857 on 124 degrees of freedom
## Multiple R-Squared: 0.983,   Adjusted R-squared: 0.9822 
## F-statistic:  1198 on 6 and 124 DF,  p-value: < 2.2e-16 
## 
## 
## 
## Covariance matrix of residuals:
##           desempleo morosidad
## desempleo 1.184e-05  1.28e-06
## morosidad 1.280e-06  3.45e-06
## 
## Correlation matrix of residuals:
##           desempleo morosidad
## desempleo    1.0000    0.2002
## morosidad    0.2002    1.0000

#Evaluación del modelo ##Evaluación de correlación serial

###Prueba de Portmanteau

bv.serial <- serial.test(modelo, lags.pt = 12, type = "PT.asymptotic")
bv.serial$serial
## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object modelo
## Chi-squared = 48.403, df = 36, p-value = 0.08109

a razón de que se obtiene un valor del 8,10%, se concluye que hay ausencia de correlación serial.

###Prueba de heterocedasticidad

bv.arch <- arch.test(modelo, lags.multi = 12, multivariate.only = TRUE)
bv.arch
## 
##  ARCH (multivariate)
## 
## data:  Residuals of VAR object modelo
## Chi-squared = 129.97, df = 108, p-value = 0.07368

el valor en esta prueba es del 7%, indicando que hay ausencia de heterocedasticidad

###Prueba de normalidad de los residuales

bv.norm <- normality.test(modelo, multivariate.only = FALSE)
bv.norm
## $desempleo
## 
##  JB-Test (univariate)
## 
## data:  Residual of desempleo equation
## Chi-squared = 1.1956, df = 2, p-value = 0.55
## 
## 
## $morosidad
## 
##  JB-Test (univariate)
## 
## data:  Residual of morosidad equation
## Chi-squared = 110.44, df = 2, p-value < 2.2e-16
## 
## 
## $JB
## 
##  JB-Test (multivariate)
## 
## data:  Residuals of VAR object modelo
## Chi-squared = 77.76, df = 4, p-value = 5.551e-16
## 
## 
## $Skewness
## 
##  Skewness only (multivariate)
## 
## data:  Residuals of VAR object modelo
## Chi-squared = 25.426, df = 2, p-value = 3.011e-06
## 
## 
## $Kurtosis
## 
##  Kurtosis only (multivariate)
## 
## data:  Residuals of VAR object modelo
## Chi-squared = 52.334, df = 2, p-value = 4.324e-12

###Ruptura estructural de los residuos

bv.cusum <- stability(modelo, type = "OLS-CUSUM")
plot(bv.cusum)

En la prueba anterior, se observa que no hay una ruptura que supere los intervalos de confianza

###Causalidad de Granger

Prueba con la variable morosidad

bv.cause.mor <- causality(modelo, cause = "morosidad")
bv.cause.mor
## $Granger
## 
##  Granger causality H0: morosidad do not Granger-cause desempleo
## 
## data:  VAR object modelo
## F-Test = 1.8124, df1 = 3, df2 = 248, p-value = 0.1454
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: morosidad and desempleo
## 
## data:  VAR object modelo
## Chi-squared = 5.0499, df = 1, p-value = 0.02463

Prueba con la variable Desempleo

bv.cause.des <- causality(modelo, cause = "desempleo")
bv.cause.des
## $Granger
## 
##  Granger causality H0: desempleo do not Granger-cause morosidad
## 
## data:  VAR object modelo
## F-Test = 1.0948, df1 = 3, df2 = 248, p-value = 0.3519
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: desempleo and morosidad
## 
## data:  VAR object modelo
## Chi-squared = 5.0499, df = 1, p-value = 0.02463
bv.cause.mor <- causality(modelo, cause = "morosidad")
bv.cause.mor
## $Granger
## 
##  Granger causality H0: morosidad do not Granger-cause desempleo
## 
## data:  VAR object modelo
## F-Test = 1.8124, df1 = 3, df2 = 248, p-value = 0.1454
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: morosidad and desempleo
## 
## data:  VAR object modelo
## Chi-squared = 5.0499, df = 1, p-value = 0.02463

##Función de impulso y respuesta

irf.des <- vars::irf(modelo, impulse = "morosidad", response = "desempleo",n.ahead = 40, boot = TRUE)
plot(irf.des) 

irf.mor <- vars::irf(modelo, impulse = "desempleo", response = "morosidad",n.ahead = 40, boot = TRUE)
plot(irf.mor)

irf.une_un <- vars::irf(modelo, impulse = "morosidad", response = "morosidad",n.ahead = 40, boot = TRUE)
#plot(irf.une, ylab = "unemployment", main = "Shock from GDP", xlab = "")
plot(irf.une_un, ylab = "Desepmleo", main = "Shock from unemployment")

###Descomposición de la varianza

bv.vardec <- vars::fevd(modelo, n.ahead = 10)
plot(bv.vardec)

predictions <- predict(modelo, n.ahead = 8, ci = 0.95)
plot(predictions, names = "desempleo")

fanchart(predictions, names = "desempleo")

plot(predictions, names = "morosidad")

fanchart(predictions, names = "morosidad")