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))

Utilizamos las variables ingreso y desempleo.

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

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

El modelo sugiere utilizar 8 lags, sin embargo al utilizar 8 el modelo presenta autocorrelación e incluye variables que no son significativas al modelo.

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

modelo_s<-summary(modelo1)
modelo_s$roots
##  [1] 0.7775660 0.7775660 0.7660914 0.7660914 0.7488061 0.7488061 0.6720236
##  [8] 0.6720236 0.6364359 0.5646411
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 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.

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.

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

En la prueba de normalidad el p.value es < 0.05 por lo que nuestro modelo es normal.

arch<-arch.test(modelo1, lags.multi = 12, multivariate.only = FALSE)
arch
## $Income
## 
##  ARCH test (univariate)
## 
## data:  Residual of Income equation
## Chi-squared = 11.621, df = 16, p-value = 0.7697
## 
## 
## $Unemployment
## 
##  ARCH test (univariate)
## 
## data:  Residual of Unemployment equation
## Chi-squared = 21.391, df = 16, p-value = 0.164
## 
## 
## 
##  ARCH (multivariate)
## 
## data:  Residuals of VAR object modelo1
## Chi-squared = 121.69, df = 108, p-value = 0.1737

El modelo no presenta heteroscedasticidad ya que el p.value es > 0.05

arch<-arch.test(modelo1, lags.multi = 12, multivariate.only = FALSE)
arch
## $Income
## 
##  ARCH test (univariate)
## 
## data:  Residual of Income equation
## Chi-squared = 11.621, df = 16, p-value = 0.7697
## 
## 
## $Unemployment
## 
##  ARCH test (univariate)
## 
## data:  Residual of Unemployment equation
## Chi-squared = 21.391, df = 16, p-value = 0.164
## 
## 
## 
##  ARCH (multivariate)
## 
## data:  Residuals of VAR object modelo1
## Chi-squared = 121.69, df = 108, p-value = 0.1737

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

Respuesta de impulso

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

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

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)

Fórmula

modelo1$varresult$Income$coefficients
##       Income.l1 Unemployment.l1       Income.l2 Unemployment.l2       Income.l3 
##     -0.13068647     -0.47812026      0.08288524      0.39250556      0.03524393 
## Unemployment.l3       Income.l4 Unemployment.l4       Income.l5 Unemployment.l5 
##     -0.65813631     -0.09359142     -0.14117449     -0.12760511      0.39764021 
##           const 
##      0.87038486
modelo1$varresult$Unemployment$coefficients
##       Income.l1 Unemployment.l1       Income.l2 Unemployment.l2       Income.l3 
##    -0.050173536     0.436086896    -0.004765779     0.156142861    -0.023081051 
## Unemployment.l3       Income.l4 Unemployment.l4       Income.l5 Unemployment.l5 
##     0.067925032    -0.012685006    -0.183030356     0.068785222     0.038700766 
##           const 
##     0.012796215
autoplot(forecast(modelo1))