Para empezar cargaremos las librerias a utilizar para el script y también los datos de la serie de tiempo. Asímismo, estaremos seleccionando el Ingreso y la Producción como las variables a utilizar dentro del script.

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,3)])

Seleccionamos dos variables con un grado alto de causalidad entre ellas. Se utilizó la serie de tiempo de “uschange” para la recopilación de las variables:

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

Luego se procedió con la búsqueda de los parámetros para el modelo

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

De acuerdo al resultado, el número de lags recomendable para el modelo fue de 3. El modelo muestra un R cuadrado aceptable y las variables son significativas por lo que procedimos a utilizar 3 lags para el modelo a pesar de que existen ciertos lags con baja o nula significancia.

modelo1<-VAR(uschange[,c(2,3)],p=3,type=c("const"))

aic2<-summary(modelo1)$logLik

summary(modelo1,equation="Production")
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: Income, Production 
## Deterministic variables: const 
## Sample size: 184 
## Log Likelihood: -523.213 
## Roots of the characteristic polynomial:
## 0.7242 0.5824 0.5824 0.5104 0.5065 0.5065
## Call:
## VAR(y = uschange[, c(2, 3)], p = 3, type = c("const"))
## 
## 
## Estimation results for equation Production: 
## =========================================== 
## Production = Income.l1 + Production.l1 + Income.l2 + Production.l2 + Income.l3 + Production.l3 + const 
## 
##               Estimate Std. Error t value Pr(>|t|)    
## Income.l1      0.05861    0.10598   0.553   0.5809    
## Production.l1  0.66060    0.07922   8.339 2.04e-14 ***
## Income.l2      0.06511    0.10526   0.619   0.5370    
## Production.l2 -0.21153    0.09338  -2.265   0.0247 *  
## Income.l3      0.17850    0.10313   1.731   0.0852 .  
## Production.l3  0.07336    0.07613   0.964   0.3366    
## const          0.04137    0.15720   0.263   0.7927    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 1.236 on 177 degrees of freedom
## Multiple R-Squared: 0.3759,  Adjusted R-squared: 0.3548 
## F-statistic: 17.77 on 6 and 177 DF,  p-value: 4.486e-16 
## 
## 
## 
## Covariance matrix of residuals:
##            Income Production
## Income     0.8129     0.3851
## Production 0.3851     1.5271
## 
## Correlation matrix of residuals:
##            Income Production
## Income     1.0000     0.3457
## Production 0.3457     1.0000
summary(modelo1,equation="Income")
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: Income, Production 
## Deterministic variables: const 
## Sample size: 184 
## Log Likelihood: -523.213 
## Roots of the characteristic polynomial:
## 0.7242 0.5824 0.5824 0.5104 0.5065 0.5065
## Call:
## VAR(y = uschange[, c(2, 3)], p = 3, type = c("const"))
## 
## 
## Estimation results for equation Income: 
## ======================================= 
## Income = Income.l1 + Production.l1 + Income.l2 + Production.l2 + Income.l3 + Production.l3 + const 
## 
##               Estimate Std. Error t value Pr(>|t|)    
## Income.l1     -0.15996    0.07732  -2.069 0.040015 *  
## Production.l1  0.02411    0.05780   0.417 0.677019    
## Income.l2      0.07870    0.07680   1.025 0.306839    
## Production.l2 -0.04228    0.06813  -0.621 0.535691    
## Income.l3      0.02028    0.07524   0.270 0.787794    
## Production.l3  0.18850    0.05555   3.394 0.000851 ***
## const          0.66458    0.11469   5.795 3.07e-08 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 0.9016 on 177 degrees of freedom
## Multiple R-Squared: 0.1028,  Adjusted R-squared: 0.07238 
## F-statistic:  3.38 on 6 and 177 DF,  p-value: 0.00352 
## 
## 
## 
## Covariance matrix of residuals:
##            Income Production
## Income     0.8129     0.3851
## Production 0.3851     1.5271
## 
## Correlation matrix of residuals:
##            Income Production
## Income     1.0000     0.3457
## Production 0.3457     1.0000

Posteriormente, utilizamos distintas pruebas para poder validar el modelo

#>PortManteu Test > 0.05 Autocorrelación
serial.test(modelo1, lags.pt=10, type="PT.asymptotic")
## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object modelo1
## Chi-squared = 42.461, df = 28, p-value = 0.03921

En este caso, debido a que no se superó la prueba, se pudo confirmar que no existe un problema de autocorrelación en el modelo.

#Raíz unitaria > 0.05
roots(modelo1)
## [1] 0.7241614 0.5823841 0.5823841 0.5104046 0.5064787 0.5064787

Debido a que todos los números son menores a 1, se pudo confirmar que la serie de tiempo es estacionaria ya que no cuenta con raíces unitarias.

#normalidad Jarque Bera < 0.05
normality.test(modelo1, multivariate.only=FALSE)
## $Income
## 
##  JB-Test (univariate)
## 
## data:  Residual of Income equation
## Chi-squared = 265.45, df = 2, p-value < 2.2e-16
## 
## 
## $Production
## 
##  JB-Test (univariate)
## 
## data:  Residual of Production equation
## Chi-squared = 51.196, df = 2, p-value = 7.637e-12
## 
## 
## $JB
## 
##  JB-Test (multivariate)
## 
## data:  Residuals of VAR object modelo1
## Chi-squared = 292, df = 4, p-value < 2.2e-16
## 
## 
## $Skewness
## 
##  Skewness only (multivariate)
## 
## data:  Residuals of VAR object modelo1
## Chi-squared = 8.543, df = 2, p-value = 0.01396
## 
## 
## $Kurtosis
## 
##  Kurtosis only (multivariate)
## 
## data:  Residuals of VAR object modelo1
## Chi-squared = 283.45, df = 2, p-value < 2.2e-16

Se realizó una prueba de Jarque Bera con la cual se pudo confirmar que los residuos presentan un comportamiento normal

#heteroscedasticity >0.05 NO HAY
arch<-arch.test(modelo1, lags.multi = 12, multivariate.only = TRUE)
arch
## 
##  ARCH (multivariate)
## 
## data:  Residuals of VAR object modelo1
## Chi-squared = 122.21, df = 108, p-value = 0.1654

De acuerdo al resultado de la prueba de heterocedasticidad, el modelo no tiene este problema ya que el p-value fue mayor a 0.05

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

Conforme a las gráficas anteriores, se pudo evidenciar que no existen quiebres estructurales dentro del modelo ya que las variables no superan los parámetros de tolerancia

Pruebas de Granger (Causalidad)

GrangerIncome <-causality(modelo1, cause = 'Income')
GrangerIncome
GrangerProduction <-causality(modelo1, cause = 'Production')
GrangerProduction

Ambas variables mostraron causalidad entre ellas.

A continuación se graficaron los tiempos en los cuales las variables se estabilizan luego de un shock con respecto a la otra variable

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

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

Descomposición de la varianza

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

Se pudo ver que existe una relación dentro del ingreso entre el ingreso y la producción, sin embargo, fue un poco mayor la relación dentro de la producción entre ambas.

Predicción

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