Cargar las librerias

library(vars)
library(tseries)
library(forecast)
library(urca)
library(highcharter)
library(bvartools)

Cargar los datos y darles formato de series de tiempo

library(rio)
datos <- rio::import("https://github.com/Wilsonsr/Series-de-Tiempo/raw/main/bases/blanchQua.csv")
gdp <- ts(datos$GDP, start = c(1948, 2), freq = 4) # gdp
une <- ts(datos$U, start = c(1948, 2), freq = 4)   # unemployment
plot(cbind(gdp, une), main="PIB y DESEMPLEO" )

hchart(cbind(gdp, une))
## Warning: Deprecated function. Use the `create_axis` function.

Pruebas de raíz unitaria

Verificar si hay una raíz unitaria, donde notamos que podemos rechazar nulos de raíz unitaria al usar la prueba Dickey-Fuller

adf.une <- ur.df(une, type = "trend", selectlags = "BIC")
summary(adf.une)   # La series es estacionaria
## 
## ############################################### 
## # 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.70448 -0.21077 -0.04731  0.19844  1.03067 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  0.0059997  0.0536401   0.112    0.911    
## z.lag.1     -0.0972035  0.0204701  -4.749 4.67e-06 ***
## tt          -0.0001167  0.0005838  -0.200    0.842    
## z.diff.lag   0.6846610  0.0595395  11.499  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3311 on 153 degrees of freedom
## Multiple R-squared:  0.4776, Adjusted R-squared:  0.4674 
## F-statistic: 46.63 on 3 and 153 DF,  p-value: < 2.2e-16
## 
## 
## Value of test-statistic is: -4.7486 7.5583 11.3199 
## 
## 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

La serie de tiempo del desempleo es estacionaria.

adf.gdp <- ur.df(gdp, type = "trend", selectlags = "BIC")
summary(adf.gdp)   
## 
## ############################################### 
## # 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 
## -3.0512 -0.6762  0.0349  0.6835  3.7335 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.0240062  0.1677018  -0.143    0.886    
## z.lag.1     -0.5580038  0.0906818  -6.153 6.35e-09 ***
## tt           0.0002905  0.0018240   0.159    0.874    
## z.diff.lag  -0.1203398  0.0802804  -1.499    0.136    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.036 on 153 degrees of freedom
## Multiple R-squared:  0.326,  Adjusted R-squared:  0.3128 
## F-statistic: 24.67 on 3 and 153 DF,  p-value: 4.456e-13
## 
## 
## Value of test-statistic is: -6.1534 12.6359 18.9522 
## 
## 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
plot(adf.gdp)

La serie de tiempo del PIB es estacionaria.

summary(ur.pp(gdp, type = "Z-tau"))
## 
## ################################## 
## # Phillips-Perron Unit Root Test # 
## ################################## 
## 
## Test regression with intercept 
## 
## 
## Call:
## lm(formula = y ~ y.l1)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -3.1487 -0.6714  0.0956  0.6718  3.7773 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -0.003187   0.082224  -0.039    0.969    
## y.l1         0.364073   0.074614   4.879  2.6e-06 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.034 on 156 degrees of freedom
## Multiple R-squared:  0.1324, Adjusted R-squared:  0.1268 
## F-statistic: 23.81 on 1 and 156 DF,  p-value: 2.605e-06
## 
## 
## Value of test-statistic, type: Z-tau  is: -8.6003 
## 
##          aux. Z statistics
## Z-tau-mu           -0.0392
## 
## Critical values for Z statistics: 
##                     1pct      5pct     10pct
## critical values -3.47264 -2.879764 -2.576381

Selección y estimación del modelo

Para crear un objeto bivariado para las dos series de tiempo que modelaremos, podemos simplemente unir en columna los dos objetos existentes.

Selección

A partir de entonces, podemos usar criterios de información para decidir el número de retrasos que se incluirán.

dat.bv <- cbind(gdp, une)
colnames(dat.bv) <- c("gdp", "une")
info.bv <- VARselect(dat.bv, lag.max = 12, type = "const")
info.bv
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      2      2      2      2 
## 
## $criteria
##                  1           2           3           4          5           6
## AIC(n) -2.87304917 -2.97227231 -2.97044481 -2.95619395 -2.9179975 -2.92377526
## HQ(n)  -2.82345549 -2.88961617 -2.85472622 -2.80741291 -2.7361540 -2.70886931
## SC(n)  -2.75099070 -2.76884152 -2.68564170 -2.59001853 -2.4704498 -2.39485521
## FPE(n)  0.05652695  0.05118955  0.05128789  0.05203247  0.0540721  0.05378026
##                  7          8           9          10          11          12
## AIC(n) -2.88043354 -2.8338250 -2.80225128 -2.81189842 -2.81480573 -2.82165222
## HQ(n)  -2.63246514 -2.5527942 -2.48815796 -2.46474265 -2.43458751 -2.40837154
## SC(n)  -2.27014118 -2.1421604 -2.02921428 -1.95748911 -1.87902411 -1.80449828
## FPE(n)  0.05619048  0.0589099  0.06085003  0.06032776  0.06022773  0.05990607
info.bv$selection   
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      2      2      2      2
library(MTS)
## 
## Attaching package: 'MTS'
## The following object is masked from 'package:vars':
## 
##     VAR
VARorder(dat.bv)
## selected order: aic =  3 
## selected order: bic =  2 
## selected order: hq =  2 
## Summary table:  
##        p     AIC     BIC      HQ     M(p) p-value
##  [1,]  0  0.5502  0.5502  0.5502   0.0000  0.0000
##  [2,]  1 -2.9115 -2.8343 -2.8802 500.4698  0.0000
##  [3,]  2 -3.0157 -2.8613 -2.9530  21.7110  0.0002
##  [4,]  3 -3.0256 -2.7940 -2.9316   8.3393  0.0799
##  [5,]  4 -3.0095 -2.7007 -2.8841   4.6658  0.3233
##  [6,]  5 -2.9712 -2.5851 -2.8144   1.6125  0.8065
##  [7,]  6 -2.9787 -2.5155 -2.7906   7.6672  0.1046
##  [8,]  7 -2.9408 -2.4003 -2.7213   1.6133  0.8064
##  [9,]  8 -2.8981 -2.2805 -2.6473   0.9842  0.9122
## [10,]  9 -2.8706 -2.1758 -2.5884   2.8840  0.5774
## [11,] 10 -2.8816 -2.1096 -2.5681   7.6368  0.1058
## [12,] 11 -2.8879 -2.0387 -2.5431   6.9351  0.1394
## [13,] 12 -2.8969 -1.9704 -2.5206   7.1415  0.1286
## [14,] 13 -2.8646 -1.8609 -2.4570   2.1373  0.7105

Todos los criterios de información sugieren el uso de 2 rezagos sería apropiado, \(p=2\)

Estimación del modelo VAR

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

summary(modelo)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: gdp, une 
## Deterministic variables: both 
## Sample size: 157 
## Log Likelihood: -211.351 
## Roots of the characteristic polynomial:
## 0.8359 0.8359 0.1096 0.1096
## Call:
## vars::VAR(y = dat.bv, p = 2, type = "both")
## 
## 
## Estimation results for equation gdp: 
## ==================================== 
## gdp = gdp.l1 + une.l1 + gdp.l2 + une.l2 + const + trend 
## 
##          Estimate Std. Error t value Pr(>|t|)   
## gdp.l1  0.0682370  0.1083678   0.630  0.52986   
## une.l1 -0.6185338  0.3093219  -2.000  0.04733 * 
## gdp.l2  0.0467653  0.0914260   0.512  0.60974   
## une.l2  0.9053891  0.3164947   2.861  0.00483 **
## const   0.0014743  0.1578306   0.009  0.99256   
## trend  -0.0002127  0.0017050  -0.125  0.90090   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.9644 on 151 degrees of freedom
## Multiple R-Squared: 0.2688,  Adjusted R-squared: 0.2446 
## F-statistic:  11.1 on 5 and 151 DF,  p-value: 4.025e-09 
## 
## 
## Estimation results for equation une: 
## ==================================== 
## une = gdp.l1 + une.l1 + gdp.l2 + une.l2 + const + trend 
## 
##          Estimate Std. Error t value Pr(>|t|)    
## gdp.l1 -0.1200753  0.0361378  -3.323 0.001118 ** 
## une.l1  1.3298154  0.1031506  12.892  < 2e-16 ***
## gdp.l2 -0.0306113  0.0304881  -1.004 0.316965    
## une.l2 -0.4130320  0.1055425  -3.913 0.000137 ***
## const   0.0119820  0.0526323   0.228 0.820223    
## trend  -0.0002412  0.0005686  -0.424 0.672003    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.3216 on 151 degrees of freedom
## Multiple R-Squared: 0.943,   Adjusted R-squared: 0.9411 
## F-statistic: 499.8 on 5 and 151 DF,  p-value: < 2.2e-16 
## 
## 
## 
## Covariance matrix of residuals:
##         gdp     une
## gdp  0.9300 -0.2036
## une -0.2036  0.1034
## 
## Correlation matrix of residuals:
##         gdp     une
## gdp  1.0000 -0.6566
## une -0.6566  1.0000

  • El sistema es estable (las raíces características pueden interpretarse como valores propios en este caso).

  • Hay variables no significativas en el modelo, donde gdp está influenciado por el desempleo pasado;

  • El desempleo está influenciado por medidas pasadas de gdp y desempleo.

#roots(modelo)

Evaluación del Modelo

correlación serial

Para probar la correlación serial, se aplica la 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 = 46.791, df = 40, p-value = 0.2137

\(p\)-valor mayor 0.05 indica ausencia de correlación serial

Test heteroscedasticidad

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

\(p\)-valor mayor al 5% indicates ausencia de heterocedasticidad

Test de Normalidad de los residuales

bv.norm <- normality.test(modelo, multivariate.only = FALSE)
bv.norm
## $gdp
## 
##  JB-Test (univariate)
## 
## data:  Residual of gdp equation
## Chi-squared = 1.139, df = 2, p-value = 0.5658
## 
## 
## $une
## 
##  JB-Test (univariate)
## 
## data:  Residual of une equation
## Chi-squared = 16.102, df = 2, p-value = 0.0003188
## 
## 
## $JB
## 
##  JB-Test (multivariate)
## 
## data:  Residuals of VAR object modelo
## Chi-squared = 2.0915, df = 4, p-value = 0.7189
## 
## 
## $Skewness
## 
##  Skewness only (multivariate)
## 
## data:  Residuals of VAR object modelo
## Chi-squared = 1.8877, df = 2, p-value = 0.3891
## 
## 
## $Kurtosis
## 
##  Kurtosis only (multivariate)
## 
## data:  Residuals of VAR object modelo
## Chi-squared = 0.20383, df = 2, p-value = 0.9031

p-value indica que los residuales son normales

Ruptura estructural en los residuos,

aplicar una prueba CUSUM.

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

No parece haber una ruptura en los respectivos intervalos de confianza.

Causalidad de Granger, IRF y descomposiciones de varianza

Causalidad de Granger

bv.cause.gdp <- causality(modelo, cause = "gdp")
bv.cause.gdp
## $Granger
## 
##  Granger causality H0: gdp do not Granger-cause une
## 
## data:  VAR object modelo
## F-Test = 5.6111, df1 = 2, df2 = 302, p-value = 0.004049
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: gdp and une
## 
## data:  VAR object modelo
## Chi-squared = 47.292, df = 1, p-value = 6.115e-12

Segun los resultados el crecimiento trimestral de gdp causa Granger el crecimiento trimestral de une, al menos en nuestro periodo de estimación.

bv.cause.une <- causality(modelo, cause = "une")
bv.cause.une
## $Granger
## 
##  Granger causality H0: une do not Granger-cause gdp
## 
## data:  VAR object modelo
## F-Test = 12.731, df1 = 2, df2 = 302, p-value = 4.917e-06
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: une and gdp
## 
## data:  VAR object modelo
## Chi-squared = 47.292, df = 1, p-value = 6.115e-12

El crecimiento trimestral de desempleo causa Granger el crecimiento trimestral del gdp.

Funcion de Impulso - Respuesta

Las funciones de impulso respuesta muestran los efectos de los shock en la trayectoria de ajuste de las variables.

Para generar funciones de respuesta a impulsos para describir la respuesta del producto a un shock de desempleo, procedemos de la siguiente manera:

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

Un choque positivo al desempleo tiene un efecto negativo en la salida (es decir, un menor poder adquisitivo).

Para considerar la respuesta del desempleo al shock de salida,

irf.une <- vars::irf(modelo, impulse = "gdp", response = "une",n.ahead = 40, boot = TRUE)
plot(irf.une) # reponse of unemployment to gdp shock

# positive shock to output decreases unemployment by large and persistent amount

Aquí observamos que un shock positivo en GDP disminuye el desempleo en una cantidad relativamente grande y persistente.

  • Del mismo modo, podemos ver el efecto de un shock de desempleo en el desempleo.
irf.une_un <- vars::irf(modelo, impulse = "une", response = "une",n.ahead = 40, boot = TRUE)
#plot(irf.une, ylab = "unemployment", main = "Shock from GDP", xlab = "")
plot(irf.une_un, ylab = "unemployment", main = "Shock from unemployment")

  • Descomposición de La varianza

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

gdp está determinado en gran medida por las perturbaciones del pib, mientras que el desempleo está influenciado por ambas perturbaciones

Pronósticos

Para pronosticar hacia adelante podemos usar el comando \(\textit{predict}\), donde en este caso pronosticamos 8 pasos adelante. También buscamos utilizar intervalos de confianza del 95% para el pronóstico

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

plot(predictions, names = "une")

fanchart(predictions, names = "gdp")

fanchart(predictions, names = "une")