Primero se carga las librerías a utilizar y la serie de tiempo.

library(vars)
## Loading required package: MASS
## Loading required package: strucchange
## Loading required package: zoo
## 
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric
## Loading required package: sandwich
## Loading required package: urca
## Loading required package: lmtest
library(fpp2)
## Registered S3 method overwritten by 'quantmod':
##   method            from
##   as.zoo.data.frame zoo
## ── Attaching packages ────────────────────────────────────────────── fpp2 2.5 ──
## ✔ ggplot2   3.4.0     ✔ fma       2.5  
## ✔ forecast  8.20      ✔ expsmooth 2.3
## 
library(TSA)
## Registered S3 methods overwritten by 'TSA':
##   method       from    
##   fitted.Arima forecast
##   plot.Arima   forecast
## 
## Attaching package: 'TSA'
## The following objects are masked from 'package:stats':
## 
##     acf, arima
## The following object is masked from 'package:utils':
## 
##     tar
series<-uschange
autoplot(uschange[,c(2,5)])

ts.plot(series[,c(2,5)], xlab="Tiempo",col=c(1,3))

Se escogen las variables con mayor causalidad así como Ingreso y Desempleo y se extraen de la serie de tiempo de uschange.

Búsqueda de parámetros y elección del modelo candidato

a <- VARselect(uschange[,c(2,5)], lag.max=15,type="const")
a$selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##      8      1      1      8

La función Varselect sufiere que el mejor modelo candidato es de 8 lags, pero al utilizar 8 lags el modelo presenta autocorrelación e incluye variables que no son significativas al modelo. por esa razón usaremos 5 lags.

modelo1<-VAR(uschange[,c(2,5)],p=5,type=c("const"))
modelo_s<-summary(modelo1)
summary(modelo1,equation="Income")
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: Income, Unemployment 
## Deterministic variables: const 
## Sample size: 182 
## Log Likelihood: -262.292 
## Roots of the characteristic polynomial:
## 0.7776 0.7776 0.7661 0.7661 0.7488 0.7488 0.672 0.672 0.6364 0.5646
## Call:
## VAR(y = uschange[, c(2, 5)], p = 5, type = c("const"))
## 
## 
## Estimation results for equation Income: 
## ======================================= 
## Income = Income.l1 + Unemployment.l1 + Income.l2 + Unemployment.l2 + Income.l3 + Unemployment.l3 + Income.l4 + Unemployment.l4 + Income.l5 + Unemployment.l5 + const 
## 
##                 Estimate Std. Error t value Pr(>|t|)    
## Income.l1       -0.13069    0.07628  -1.713  0.08848 .  
## Unemployment.l1 -0.47812    0.22736  -2.103  0.03693 *  
## Income.l2        0.08289    0.07790   1.064  0.28883    
## Unemployment.l2  0.39251    0.23868   1.644  0.10191    
## Income.l3        0.03524    0.07602   0.464  0.64351    
## Unemployment.l3 -0.65814    0.24414  -2.696  0.00773 ** 
## Income.l4       -0.09359    0.07556  -1.239  0.21717    
## Unemployment.l4 -0.14117    0.24531  -0.575  0.56571    
## Income.l5       -0.12761    0.07471  -1.708  0.08944 .  
## Unemployment.l5  0.39764    0.22214   1.790  0.07522 .  
## const            0.87038    0.14889   5.846 2.49e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.8936 on 171 degrees of freedom
## Multiple R-Squared: 0.1347,  Adjusted R-squared: 0.08413 
## F-statistic: 2.663 on 10 and 171 DF,  p-value: 0.004781 
## 
## 
## 
## Covariance matrix of residuals:
##               Income Unemployment
## Income        0.7985     -0.05280
## Unemployment -0.0528      0.09033
## 
## Correlation matrix of residuals:
##               Income Unemployment
## Income        1.0000      -0.1966
## Unemployment -0.1966       1.0000
summary(modelo1,equation="Unemployment")
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: Income, Unemployment 
## Deterministic variables: const 
## Sample size: 182 
## Log Likelihood: -262.292 
## Roots of the characteristic polynomial:
## 0.7776 0.7776 0.7661 0.7661 0.7488 0.7488 0.672 0.672 0.6364 0.5646
## Call:
## VAR(y = uschange[, c(2, 5)], p = 5, type = c("const"))
## 
## 
## Estimation results for equation Unemployment: 
## ============================================= 
## Unemployment = Income.l1 + Unemployment.l1 + Income.l2 + Unemployment.l2 + Income.l3 + Unemployment.l3 + Income.l4 + Unemployment.l4 + Income.l5 + Unemployment.l5 + const 
## 
##                  Estimate Std. Error t value Pr(>|t|)    
## Income.l1       -0.050174   0.025657  -1.956  0.05215 .  
## Unemployment.l1  0.436087   0.076472   5.703 5.08e-08 ***
## Income.l2       -0.004766   0.026202  -0.182  0.85589    
## Unemployment.l2  0.156143   0.080280   1.945  0.05342 .  
## Income.l3       -0.023081   0.025569  -0.903  0.36796    
## Unemployment.l3  0.067925   0.082118   0.827  0.40930    
## Income.l4       -0.012685   0.025414  -0.499  0.61833    
## Unemployment.l4 -0.183030   0.082510  -2.218  0.02785 *  
## Income.l5        0.068785   0.025128   2.737  0.00685 ** 
## Unemployment.l5  0.038701   0.074718   0.518  0.60516    
## const            0.012796   0.050080   0.256  0.79863    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.3006 on 171 degrees of freedom
## Multiple R-Squared: 0.3665,  Adjusted R-squared: 0.3295 
## F-statistic: 9.893 on 10 and 171 DF,  p-value: 5.479e-13 
## 
## 
## 
## Covariance matrix of residuals:
##               Income Unemployment
## Income        0.7985     -0.05280
## Unemployment -0.0528      0.09033
## 
## Correlation matrix of residuals:
##               Income Unemployment
## Income        1.0000      -0.1966
## Unemployment -0.1966       1.0000

Hay ciertos lags de las variables que no son significativos al modelo, por lo que se podria considerar en eliminarlos del modelo.

Pruebas del Modelo

serial.test(modelo1, lags.pt=10, type="PT.asymptotic")
## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object modelo1
## Chi-squared = 31.773, df = 20, p-value = 0.04577

El modelo no presenta autocorrelación porque según el resultado de la prueba de PortManteu Test > 0.05 Autocorrelación. Básicamente el p-value tiene que ser mayor para pasar la prueba.

roots(modelo1)
##  [1] 0.7775660 0.7775660 0.7660914 0.7660914 0.7488061 0.7488061 0.6720236
##  [8] 0.6720236 0.6364359 0.5646411

El modelo no tiene raices unitarias por lo que es estacionaria. Para pasar esta prueba todos los números tienen que ser < 1.

normality.test(modelo1, multivariate.only=FALSE)
## $Income
## 
##  JB-Test (univariate)
## 
## data:  Residual of Income equation
## Chi-squared = 195.51, df = 2, p-value < 2.2e-16
## 
## 
## $Unemployment
## 
##  JB-Test (univariate)
## 
## data:  Residual of Unemployment equation
## Chi-squared = 33.197, df = 2, p-value = 6.184e-08
## 
## 
## $JB
## 
##  JB-Test (multivariate)
## 
## data:  Residuals of VAR object modelo1
## Chi-squared = 216.79, df = 4, p-value < 2.2e-16
## 
## 
## $Skewness
## 
##  Skewness only (multivariate)
## 
## data:  Residuals of VAR object modelo1
## Chi-squared = 7.8916, df = 2, p-value = 0.01934
## 
## 
## $Kurtosis
## 
##  Kurtosis only (multivariate)
## 
## data:  Residuals of VAR object modelo1
## Chi-squared = 208.9, df = 2, p-value < 2.2e-16

Si en la prueba de Jarque-Bera el p-value < 0.05 presenta normalidad en los residuos. Para este modelo el resultado presenta normalidad.

stab <-stability(modelo1, type = "OLS-CUSUM")
par(mar=c(1,1,1,1))
plot(stab)

El modelo no presenta quiebres estructurales, ambas gráficas no salen de los parámetros tolerados.

Causalidad de granger

granger < 0.05 para que exista causalidad o que tiene relación

GrangerIncome <-causality(modelo1, cause = 'Income')
GrangerIncome
## $Granger
## 
##  Granger causality H0: Income do not Granger-cause Unemployment
## 
## data:  VAR object modelo1
## F-Test = 2.8057, df1 = 5, df2 = 342, p-value = 0.01687
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: Income and Unemployment
## 
## data:  VAR object modelo1
## Chi-squared = 6.7734, df = 1, p-value = 0.009253
GrangerUnemployment <-causality(modelo1, cause = 'Unemployment')
GrangerUnemployment
## $Granger
## 
##  Granger causality H0: Unemployment do not Granger-cause Income
## 
## data:  VAR object modelo1
## F-Test = 3.7264, df1 = 5, df2 = 342, p-value = 0.002681
## 
## 
## $Instant
## 
##  H0: No instantaneous causality between: Unemployment and Income
## 
## data:  VAR object modelo1
## Chi-squared = 6.7734, df = 1, p-value = 0.009253

Ambas variables presentan causalidad.

Respuesta de impulso

UnemploymentIRF <- irf(modelo1,  impulse = "Income", response="Unemployment", n.ahead = 20, boot = T )
plot(UnemploymentIRF, ylab = "Unemployment", main = "Shock desde Income")

IncomeIRF <- irf(modelo1,  impulse = "Unemployment", response="Income", n.ahead = 20, boot = T )
plot(IncomeIRF, ylab = "Income", main = "Shock desde Unemployment")

Muestra como ambas variables en sus gráficas se estabilizan y vuelven a la normalidad después del shock.

Descomposición de la varianza

FEVD1 <- fevd(modelo1, n.ahead = 10)
plot(FEVD1)

En el ingreso se observa que aunque en menor medida tiene relación con el desempleo. De igual forma para el desempleo.

Predicción

fore<-predict(modelo1, n.ahead = 10, ci=0.95)
fanchart(fore)