Series de Tiempo Multivariantes

Caso de Estudio: Datos macroeconómicos de Estados Unidos

Pasos:

  1. Análisis exploratorio de los datos.
  2. Dividir la serie en conjuntos de entrenamiento y prueba.
  3. Prueba de estacionariedad.
  4. Transformar la serie de entrenamiento si es necesario.
  5. Construir un modelo VAR sobre las series transformadas.
  6. Causalidad de Granger.
  7. Diagnóstico del modelo.
  8. Realizar pronósticos utilizando el modelo finalmente elegido.
  9. Transformación inversa del pronóstico a la escala original.
  10. Realizar una evaluación del pronóstico.

La siguiente data tiene información trimestral del PIB y el Desempleo en estados unidos.

data=read_csv("C:/Users/Usuario/Desktop/MACROECONOMETRIA/VAR/blanchQua.csv")
## Rows: 159 Columns: 3
## ── Column specification ────────────────────────────────────────────────────────
## Delimiter: ","
## dbl (3): dates, GDP, U
## 
## ℹ Use `spec()` to retrieve the full column specification for this data.
## ℹ Specify the column types or set `show_col_types = FALSE` to quiet this message.
head(data)
## # A tibble: 6 × 3
##   dates     GDP       U
##   <dbl>   <dbl>   <dbl>
## 1 1948.  0.854  -0.144 
## 2 1948.  0.0121 -0.0675
## 3 1948. -0.0851 -0.0245
## 4 1949. -2.08    0.785 
## 5 1949. -1.49    1.96  
## 6 1949. -0.149   2.77

El paso sigueinte es convertir a formato de series de tiempo la variables de la data.

gdp=ts(data$GDP, start= c(1948,2), frequency = 4)
une=ts(data$U, start= c(1948,2), frequency = 4)

Graficamos

dat.bv = cbind(gdp, une)
autoplot(dat.bv)

Dividamos la data en un set para entrenamiento y prubea

n_obs=10 ## estos 10 lo dejamos para contrastar los pronicticos con los ultimos 10 datos reales lo que es alrededor de 3 años
end=dim(dat.bv)[1] ## aqui me cuenta la cantidad de datos de la fila 1
X_train = dat.bv[1:(end-n_obs),] ## aqui me cuenta la cantidad de datos de la fila menos los 10 datos iniciales, es decir 159 - 10
X_test = dat.bv[(end-n_obs+1):end,] ## aqui esocjo los datos que vienen a partir del 149 + 10, es decir 159. y esto va ser mis datos de prueba
dim(X_test)
## [1] 10  2
dim(X_train)
## [1] 149   2

Prueba de estacionariedad

apply(X_train, 2, adf.test) ## se coloca el 2 para especificar que lo queremos aplicar por columnas
## Warning in FUN(newX[, i], ...): p-value smaller than printed p-value
## $gdp
## 
##  Augmented Dickey-Fuller Test
## 
## data:  newX[, i]
## Dickey-Fuller = -4.7992, Lag order = 5, p-value = 0.01
## alternative hypothesis: stationary
## 
## 
## $une
## 
##  Augmented Dickey-Fuller Test
## 
## data:  newX[, i]
## Dickey-Fuller = -3.4885, Lag order = 5, p-value = 0.0459
## alternative hypothesis: stationary

Bajo el anterior resultado los p-valor tanto para GDP y Desempleo son menores a 5%, por tal motivo no hay necesidad de diferenciar.

VAR MODELO

VARselect(X_train, type= "none", lag.max = 10) ## identificar el orden del modelo, el none hace referencia a ningun modelo en espesifico, y el 10 a la cantidad maxima de modleos
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      2      2      1      2 
## 
## $criteria
##                  1           2           3           4           5           6
## AIC(n) -2.82819854 -2.90791028 -2.89783899 -2.86623464 -2.82935387 -2.82623785
## HQ(n)  -2.79388222 -2.83927764 -2.79489003 -2.72896935 -2.65777227 -2.62033992
## SC(n)  -2.74375325 -2.73901970 -2.64450312 -2.52845346 -2.40712741 -2.31956609
## FPE(n)  0.05911949  0.05459142  0.05514817  0.05692732  0.05908041  0.05928636
##                  7           8          9          10
## AIC(n) -2.78060501 -2.73652219 -2.7052103 -2.72286531
## HQ(n)  -2.54039077 -2.46199163 -2.3963634 -2.37970211
## SC(n)  -2.18948796 -2.06095985 -1.9452027 -1.87841238
## FPE(n)  0.06208603  0.06492826  0.0670521  0.06595099
Selección del Orden del Modelo VAR

Bajo lo anterior, el mejor orden del modelo de acuerdo a los 4 criterios sugiere una recomendación de AIC con 2 rezagos.

En un modelo de Vectores Autorregresivos (VAR), los rezagos representan los valores pasados de las variables endógenas (por ejemplo, el PIB y el desempleo), y se usan para capturar cómo el pasado influye en el presente.

Ejemplos:
  • gdp.l1 = PIB en el periodo anterior (t - 1)
  • une.l2 = Desempleo dos periodos atrás (t - 2)
¿Cuál es el objetivo?

El objetivo es analizar cómo los valores pasados de una variable afectan:

  • Su propio valor actual
  • El valor actual de otras variables 
var.a=vars::VAR(X_train, lag.max = 10, ic="AIC", type = "const") 
summary(var.a)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: gdp, une 
## Deterministic variables: const 
## Sample size: 147 
## Log Likelihood: -202.613 
## Roots of the characteristic polynomial:
## 0.8335 0.8335 0.1612 0.04178
## Call:
## vars::VAR(y = X_train, type = "const", lag.max = 10, ic = "AIC")
## 
## 
## Estimation results for equation gdp: 
## ==================================== 
## gdp = gdp.l1 + une.l1 + gdp.l2 + une.l2 + const 
## 
##        Estimate Std. Error t value Pr(>|t|)  
## gdp.l1  0.08324    0.11169   0.745   0.4573  
## une.l1 -0.54573    0.32008  -1.705   0.0904 .
## gdp.l2  0.05625    0.09471   0.594   0.5535  
## une.l2  0.84113    0.32662   2.575   0.0110 *
## const  -0.03378    0.08092  -0.418   0.6769  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.9794 on 142 degrees of freedom
## Multiple R-Squared: 0.2753,  Adjusted R-squared: 0.2549 
## F-statistic: 13.49 on 4 and 142 DF,  p-value: 2.417e-09 
## 
## 
## Estimation results for equation une: 
## ==================================== 
## une = gdp.l1 + une.l1 + gdp.l2 + une.l2 + const 
## 
##        Estimate Std. Error t value Pr(>|t|)    
## gdp.l1 -0.12717    0.03730  -3.409 0.000849 ***
## une.l1  1.30855    0.10691  12.240  < 2e-16 ***
## gdp.l2 -0.03194    0.03163  -1.010 0.314311    
## une.l2 -0.39444    0.10909  -3.616 0.000415 ***
## const   0.00129    0.02703   0.048 0.961997    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.3271 on 142 degrees of freedom
## Multiple R-Squared: 0.9427,  Adjusted R-squared: 0.9411 
## F-statistic: 583.8 on 4 and 142 DF,  p-value: < 2.2e-16 
## 
## 
## 
## Covariance matrix of residuals:
##         gdp     une
## gdp  0.9592 -0.2116
## une -0.2116  0.1070
## 
## Correlation matrix of residuals:
##         gdp     une
## gdp  1.0000 -0.6606
## une -0.6606  1.0000
Como interpretar VAR
  • ¿El desempleo del trimestre anterior (une.l1) tiene algún efecto sobre el PIB actual?
  • ¿El PIB de hace dos trimestres (gdp.l2) afecta al desempleo hoy?
IMPORTANTE: Limite es hasta el segundo rezago (lag 2) dado el AIC

Esto implica que se están analizando efectos con hasta dos periodos de retraso.

Significado de los términos:
Término Significado
gdp.l PIB del periodo anterior
gdp.l2 PIB de dos periodos atrás
une.l1 Desempleo del periodo anterior
une.l2 Desempleo de dos periodos atrás

Análisis de Resultados del Modelo VAR

PIB (GDP)

Para la primera ecuación, es decir, el PIB (GDP), se observa significancia en las variables une.l2.

El rezago 2 del desempleo (une.l2) afecta positivamente al PIB.

El modelo explica el 25.5% de la variabilidad del PIB. Sin embargo, en conclusión, el modelo VAR resulta débil para modelar adecuadamente el comportamiento del PIB.

Desempleo (UNE)

Para la segunda ecuación, es decir, el desempleo (UNE), se encuentra significancia en las variables gdp.l1, une.l1 y une.l2.

El desempleo actual depende positiva y fuertemente de su rezago 1 (une.l1).

Por otro lado, el PIB rezagado (gdp.l1) reduce significativamente el desempleo, lo cual es coherente desde el punto de vista económico.

Finalmente, el rezago 2 del desempleo (une.l2) también presenta un efecto negativo y significativo.

En este caso, el modelo explica más del 94% de la variabilidad del desempleo, lo que indica un ajuste más adecuado en comparación con la ecuación del PIB.

Causalidad de GRANGER

causality(var.a, cause = c("gdp"))
## $Granger
## 
##  Granger causality H0: gdp do not Granger-cause une
## 
## data:  VAR object var.a
## F-Test = 5.8894, df1 = 2, df2 = 284, p-value = 0.003118
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: gdp and une
## 
## data:  VAR object var.a
## Chi-squared = 44.656, df = 1, p-value = 2.349e-11

Bajo lo anterior, dado que el p-valor es menor al 5% (0,3%), rechazamos la hipótesis nula y consideramos con 95% que los valores pasados del PIB ayudan a predecir el desempleo, es decir, hay una relación causal con rezago en el tiempo. Por otro lado, el p-valor de Instant es viblimente menor que 5% indicando que hay causalidad instantánea entre PIB y desempleo. Es decir, ambas variables están correlacionadas en el mismo periodo.

causality(var.a, cause = c("une"))
## $Granger
## 
##  Granger causality H0: une do not Granger-cause gdp
## 
## data:  VAR object var.a
## F-Test = 12.411, df1 = 2, df2 = 284, p-value = 6.803e-06
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: une and gdp
## 
## data:  VAR object var.a
## Chi-squared = 44.656, df = 1, p-value = 2.349e-11

Al igaul que el anterior, rechazamos la hipótesis nula y por ende aceptamos que el desempleo pasado ayuda a predecir el PIB. Es decir hay causalidad. Asi mismo, hay evidencia suficiente de una correlacion de las variables.

Evaluar el modelo

bv.serial=serial.test(var.a)
bv.serial
## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object var.a
## Chi-squared = 57.677, df = 56, p-value = 0.413

Bajo lo anterior el p-valor es mayor al 5%lo cual es ideal porque eso es lo que bsucamos. Si esto no fuese asi tendriamos que cambiar el orden del modelo, el tipo de modelo o aplicar diferenciacion o utilizar logaritmos. En este caso no hay que hacer nada.

plot(bv.serial, names = "gdp")

plot(bv.serial, names = "une")

bv.norm <- normality.test(var.a, multivariate.only = FALSE)
bv.norm
## $gdp
## 
##  JB-Test (univariate)
## 
## data:  Residual of gdp equation
## Chi-squared = 0.63924, df = 2, p-value = 0.7264
## 
## 
## $une
## 
##  JB-Test (univariate)
## 
## data:  Residual of une equation
## Chi-squared = 12.941, df = 2, p-value = 0.001549
## 
## 
## $JB
## 
##  JB-Test (multivariate)
## 
## data:  Residuals of VAR object var.a
## Chi-squared = 1.6378, df = 4, p-value = 0.802
## 
## 
## $Skewness
## 
##  Skewness only (multivariate)
## 
## data:  Residuals of VAR object var.a
## Chi-squared = 1.5251, df = 2, p-value = 0.4665
## 
## 
## $Kurtosis
## 
##  Kurtosis only (multivariate)
## 
## data:  Residuals of VAR object var.a
## Chi-squared = 0.11261, df = 2, p-value = 0.9452

En cuanto a la normalidad, aunque los residuos de la ecuación del desempleo no son normales (0.001549 - El p-valor es muy bajo) según el test de Jarque-Bera univariado, el test multivariado aplicado al sistema completo confirma que los residuos del modelo VAR (72,6%), en conjunto, sí siguen una distribución normal, sin problemas de asimetría ni curtosis. Estos resultados respaldan la validez del modelo y permiten avanzar con confianza hacia el análisis estructural, como funciones impulso-respuesta y descomposición de varianza. IMPORTANTE: El VAR no exige normalidad estricta para la estimación de los coeficientes.

const.cusum = stability(var.a, type = "OLS-CUSUM")
plot(const.cusum)

bajo lo naterior podemos evidenciar que los coeficientes del modelo se mantienen estables en todo momento, esto es asi porquer se mantiene dentro de las bandas. Es importante que el VAR sea estable dado que asi permite confiar en las predicciones, análisis impulso-respuesta y descomposición de varianza.

Pronosticar

predicciones = predict (var.a, n.ahead = 15, ci=0.95) ## prediccion de 10 periodos con 95% de intervalo confianza
plot(predicciones, names="gdp")

predicciones = predict (var.a, n.ahead = 15, ci=0.95) ## prediccion de 10 periodos con 95% de intervalo confianza
plot(predicciones, names="une")

Evaluemos el modelo

pred=predicciones$fcst
rmsegdp=sqrt(mean(X_test[,1]-pred$gdp)^2) ##error cuadratico medio de gdp
rmseune=sqrt(mean(X_test[,2]-pred$une)^2) ##error cuadratico medio de une

rmsegdp
## [1] 0.2788738
rmseune
## [1] 1.36005

reentrenemos el modelo con todo

VARselect(dat.bv, type = "none", lag.max = 10)
## $selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      2      2      2      2 
## 
## $criteria
##                  1           2           3           4           5           6
## AIC(n) -2.88055345 -2.97522838 -2.96572378 -2.93483896 -2.89661887 -2.89965712
## HQ(n)  -2.84778959 -2.90970066 -2.86743221 -2.80378352 -2.73279957 -2.70307397
## SC(n)  -2.79991059 -2.81394267 -2.72379522 -2.61226754 -2.49340460 -2.41580000
## FPE(n)  0.05610388  0.05103709  0.05152765  0.05315026  0.05523189  0.05508059
##                  7          8           9          10
## AIC(n) -2.85864313 -2.8191397 -2.78690253 -2.80496021
## HQ(n)  -2.62929613 -2.5570289 -2.49202780 -2.47732162
## SC(n)  -2.29414316 -2.1739969 -2.06111684 -1.99853167
## FPE(n)  0.05741032  0.0597565  0.06175805  0.06070675
var.a = vars::VAR(dat.bv, lag.max = 10, ic = "AIC", type = "const") ## escogimos todos los datos "dat.bv" a diferencia del anterior
summary(var.a)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: gdp, une 
## Deterministic variables: const 
## Sample size: 157 
## Log Likelihood: -211.593 
## Roots of the characteristic polynomial:
## 0.8361 0.8361 0.1014 0.1014
## Call:
## vars::VAR(y = dat.bv, type = "const", lag.max = 10, ic = "AIC")
## 
## 
## Estimation results for equation gdp: 
## ==================================== 
## gdp = gdp.l1 + une.l1 + gdp.l2 + une.l2 + const 
## 
##        Estimate Std. Error t value Pr(>|t|)   
## gdp.l1  0.06905    0.10782   0.640  0.52289   
## une.l1 -0.61530    0.30723  -2.003  0.04699 * 
## gdp.l2  0.04741    0.09098   0.521  0.60309   
## une.l2  0.90193    0.31426   2.870  0.00469 **
## const  -0.01570    0.07684  -0.204  0.83833   
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.9613 on 152 degrees of freedom
## Multiple R-Squared: 0.2687,  Adjusted R-squared: 0.2495 
## F-statistic: 13.97 on 4 and 152 DF,  p-value: 1e-09 
## 
## 
## Estimation results for equation une: 
## ==================================== 
## une = gdp.l1 + une.l1 + gdp.l2 + une.l2 + const 
## 
##         Estimate Std. Error t value Pr(>|t|)    
## gdp.l1 -0.119156   0.035975  -3.312 0.001157 ** 
## une.l1  1.333482   0.102510  13.008  < 2e-16 ***
## gdp.l2 -0.029883   0.030358  -0.984 0.326497    
## une.l2 -0.416951   0.104853  -3.977 0.000108 ***
## const  -0.007501   0.025639  -0.293 0.770243    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.3207 on 152 degrees of freedom
## Multiple R-Squared: 0.9429,  Adjusted R-squared: 0.9414 
## F-statistic: 628.1 on 4 and 152 DF,  p-value: < 2.2e-16 
## 
## 
## 
## Covariance matrix of residuals:
##         gdp     une
## gdp  0.9240 -0.2022
## une -0.2022  0.1029
## 
## Correlation matrix of residuals:
##         gdp     une
## gdp  1.0000 -0.6558
## une -0.6558  1.0000

Causalidad de Granger

causality(var.a, cause = c("gdp"))
## $Granger
## 
##  Granger causality H0: gdp do not Granger-cause une
## 
## data:  VAR object var.a
## F-Test = 5.5723, df1 = 2, df2 = 304, p-value = 0.0042
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: gdp and une
## 
## data:  VAR object var.a
## Chi-squared = 47.214, df = 1, p-value = 6.364e-12
causality(var.a, cause = c("une"))
## $Granger
## 
##  Granger causality H0: une do not Granger-cause gdp
## 
## data:  VAR object var.a
## F-Test = 12.821, df1 = 2, df2 = 304, p-value = 4.511e-06
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: une and gdp
## 
## data:  VAR object var.a
## Chi-squared = 47.214, df = 1, p-value = 6.364e-12

Diagnostico del modelo

bv.serial= serial.test(var.a)
bv.serial
## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object var.a
## Chi-squared = 59.717, df = 56, p-value = 0.3422

Graficamos

plot(bv.serial, names = "gdp")

plot(bv.serial, names = "une")

### GDP

predictions = predict(var.a, n.ahead = 10, ci = 0.95)
plot(predictions, names = "gdp")

#### Ver el pronóstico de gdp

predictions$fcst$gdp
##             fcst     lower    upper       CI
##  [1,] -0.2534468 -2.137477 1.630583 1.884030
##  [2,] -0.4389495 -2.383709 1.505810 1.944759
##  [3,] -0.6424544 -2.599747 1.314838 1.957293
##  [4,] -0.7096998 -2.673288 1.253888 1.963588
##  [5,] -0.6933338 -2.685233 1.298566 1.991900
##  [6,] -0.6178004 -2.656915 1.421315 2.039115
##  [7,] -0.5084496 -2.600303 1.583404 2.091854
##  [8,] -0.3859498 -2.524352 1.752452 2.138402
##  [9,] -0.2660481 -2.438882 1.906786 2.172834
## [10,] -0.1594979 -2.353994 2.034998 2.194496

Desempleo

predictions = predict(var.a, n.ahead = 10, ci = 0.95)
plot(predictions, names = "une")

#### Ver el pronóstico de une (desempleo)

predictions$fcst$une
##              fcst     lower      upper        CI
##  [1,] -1.79276328 -2.421379 -1.1641480 0.6286153
##  [2,] -1.67808411 -2.859196 -0.4969718 1.1811123
##  [3,] -1.43782351 -3.106645  0.2309985 1.6688220
##  [4,] -1.13546379 -3.186067  0.9151394 2.0506032
##  [5,] -0.81835556 -3.139340  1.5026286 2.3209842
##  [6,] -0.52150745 -3.014598  1.9715831 2.4930906
##  [7,] -0.26737394 -2.856749  2.3220011 2.5893751
##  [8,] -0.06754976 -2.702063  2.5669636 2.6345134
##  [9,]  0.07508649 -2.575370  2.7255431 2.6504566
## [10,]  0.16402467 -2.489540  2.8175896 2.6535649

Evaluamos

pred2=predictions$fcst
rmsegdp=sqrt(mean(X_test[,1]-pred2$gdp)^2) ##error cuadratico medio de gdp
rmseune=sqrt(mean(X_test[,2]-pred2$une)^2) ##error cuadratico medio de une

rmsegdp
## [1] 0.07175598
rmseune
## [1] 0.7176352

Impulso respuesta

impulso_gdp_une = irf(var.a, impulse = "gdp", response = "une", n.ahead = 10, boot = TRUE, ci = 0.95)
plot(impulso_gdp_une)

### Impuso Respuesta

impulso_une_gdp = irf(var.a, impulse = "une", response = "gdp", n.ahead = 10, boot = TRUE, ci = 0.95)
plot(impulso_une_gdp)

Hay una caída significativa en el PIB tras un aumento inesperado en el desempleo (más desempleo = menor producción). De forma posterior el PIB se recupera e incluso aumenta por encima de cero, aunque el efecto va disminuyendo.