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.3.6 ✔ 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
#Cargar datos Consumo y producción
series2<-uschange
autoplot(uschange[,c(2,4)])
#plot de serie de datos
ts.plot(series2[,c(2,4)], xlab="Tiempo",col=c(1,2))
#Búsqueda de parámetros
a <- VARselect(uschange[,c(2,4)], lag.max=15,type="const")
a$selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 5 1 1 5
#Creación de modelo
modelo1<-VAR(uschange[,c(2,4)],p=3,type=c("const")) #el P = se puede modificar. Usualmente, entre más mejor.
modelo_s<-summary(modelo1)
#si no pasan de 1, el modelo es estacionario
modelo_s$roots
## [1] 0.7534347 0.5434097 0.5434097 0.5314931 0.3796348 0.3796348
summary(modelo1,equation="Consumption") #Produccion no parece tan causal para consumo
## Warning in summary.varest(modelo1, equation = "Consumption"):
## Invalid variable name(s) supplied, using first variable.
##
## VAR Estimation Results:
## =========================
## Endogenous variables: Income, Savings
## Deterministic variables: const
## Sample size: 184
## Log Likelihood: -896.742
## Roots of the characteristic polynomial:
## 0.7534 0.5434 0.5434 0.5315 0.3796 0.3796
## Call:
## VAR(y = uschange[, c(2, 4)], p = 3, type = c("const"))
##
##
## Estimation results for equation Income:
## =======================================
## Income = Income.l1 + Savings.l1 + Income.l2 + Savings.l2 + Income.l3 + Savings.l3 + const
##
## Estimate Std. Error t value Pr(>|t|)
## Income.l1 0.076348 0.104471 0.731 0.46586
## Savings.l1 -0.022760 0.007189 -3.166 0.00182 **
## Income.l2 0.012334 0.106747 0.116 0.90814
## Savings.l2 -0.006541 0.007375 -0.887 0.37632
## Income.l3 0.294169 0.104897 2.804 0.00560 **
## Savings.l3 -0.023609 0.007099 -3.326 0.00107 **
## const 0.498738 0.114920 4.340 2.39e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.8766 on 177 degrees of freedom
## Multiple R-Squared: 0.1519, Adjusted R-squared: 0.1231
## F-statistic: 5.283 on 6 and 177 DF, p-value: 4.934e-05
##
##
##
## Covariance matrix of residuals:
## Income Savings
## Income 0.7684 8.264
## Savings 8.2637 171.343
##
## Correlation matrix of residuals:
## Income Savings
## Income 1.0000 0.7202
## Savings 0.7202 1.0000
summary(modelo1,equation="Production") # consumo falla
## Warning in summary.varest(modelo1, equation = "Production"):
## Invalid variable name(s) supplied, using first variable.
##
## VAR Estimation Results:
## =========================
## Endogenous variables: Income, Savings
## Deterministic variables: const
## Sample size: 184
## Log Likelihood: -896.742
## Roots of the characteristic polynomial:
## 0.7534 0.5434 0.5434 0.5315 0.3796 0.3796
## Call:
## VAR(y = uschange[, c(2, 4)], p = 3, type = c("const"))
##
##
## Estimation results for equation Income:
## =======================================
## Income = Income.l1 + Savings.l1 + Income.l2 + Savings.l2 + Income.l3 + Savings.l3 + const
##
## Estimate Std. Error t value Pr(>|t|)
## Income.l1 0.076348 0.104471 0.731 0.46586
## Savings.l1 -0.022760 0.007189 -3.166 0.00182 **
## Income.l2 0.012334 0.106747 0.116 0.90814
## Savings.l2 -0.006541 0.007375 -0.887 0.37632
## Income.l3 0.294169 0.104897 2.804 0.00560 **
## Savings.l3 -0.023609 0.007099 -3.326 0.00107 **
## const 0.498738 0.114920 4.340 2.39e-05 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.8766 on 177 degrees of freedom
## Multiple R-Squared: 0.1519, Adjusted R-squared: 0.1231
## F-statistic: 5.283 on 6 and 177 DF, p-value: 4.934e-05
##
##
##
## Covariance matrix of residuals:
## Income Savings
## Income 0.7684 8.264
## Savings 8.2637 171.343
##
## Correlation matrix of residuals:
## Income Savings
## Income 1.0000 0.7202
## Savings 0.7202 1.0000
#significancia solo si es menor a 0.05, ambos aprueban
#Validación del 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 = 34.837, df = 28, p-value = 0.1747
#Raíz unitaria < 1
roots(modelo1)
## [1] 0.7534347 0.5434097 0.5434097 0.5314931 0.3796348 0.3796348
#normalidad Jarque Bera < 0.05 validamos normalidad
normality.test(modelo1, multivariate.only=FALSE) #normalidad con 3 lags
## $Income
##
## JB-Test (univariate)
##
## data: Residual of Income equation
## Chi-squared = 242.03, df = 2, p-value < 2.2e-16
##
##
## $Savings
##
## JB-Test (univariate)
##
## data: Residual of Savings equation
## Chi-squared = 146.88, df = 2, p-value < 2.2e-16
##
##
## $JB
##
## JB-Test (multivariate)
##
## data: Residuals of VAR object modelo1
## Chi-squared = 288.93, df = 4, p-value < 2.2e-16
##
##
## $Skewness
##
## Skewness only (multivariate)
##
## data: Residuals of VAR object modelo1
## Chi-squared = 2.4701, df = 2, p-value = 0.2908
##
##
## $Kurtosis
##
## Kurtosis only (multivariate)
##
## data: Residuals of VAR object modelo1
## Chi-squared = 286.46, df = 2, p-value < 2.2e-16
#heteroscedasticity >0.05 NO HAY
arch<-arch.test(modelo1, lags.multi = 39, multivariate.only = FALSE) # si se pueden cmabiar los lags multi porque son para la prueba
arch
## $Income
##
## ARCH test (univariate)
##
## data: Residual of Income equation
## Chi-squared = 15.556, df = 16, p-value = 0.4843
##
##
## $Savings
##
## ARCH test (univariate)
##
## data: Residual of Savings equation
## Chi-squared = 25.85, df = 16, p-value = 0.05618
##
##
##
## ARCH (multivariate)
##
## data: Residuals of VAR object modelo1
## Chi-squared = 347.56, df = 351, p-value = 0.5418
#Structural breaks
stab<-stability(modelo1, type = "OLS-CUSUM")
par(mar=c(1,1,1,1))
plot(stab) #si las lineas estan estnre las lineas de protuberancia creo, las variables son estables.
#BoxCox.ar(abs(uschange[,2]))
#Causalidad de granger
#granger < 0.05 para que exista causalidad
GrangerConsumptions <-causality(modelo1, cause = 'Income')
GrangerConsumptions #aprueba
## $Granger
##
## Granger causality H0: Income do not Granger-cause Savings
##
## data: VAR object modelo1
## F-Test = 1.4203, df1 = 3, df2 = 354, p-value = 0.2365
##
##
## $Instant
##
## H0: No instantaneous causality between: Income and Savings
##
## data: VAR object modelo1
## Chi-squared = 62.841, df = 1, p-value = 2.22e-15
GrangerProduction <-causality(modelo1, cause = 'Savings')
GrangerProduction #Esta si parece tener causalidad por ser menor a 0.05.
## $Granger
##
## Granger causality H0: Savings do not Granger-cause Income
##
## data: VAR object modelo1
## F-Test = 8.5558, df1 = 3, df2 = 354, p-value = 1.688e-05
##
##
## $Instant
##
## H0: No instantaneous causality between: Savings and Income
##
## data: VAR object modelo1
## Chi-squared = 62.841, df = 1, p-value = 2.22e-15
#Respuesta de impulso
#Como se comporta una variable si la otra variable recibe un "shock"
IncomeIRF <- irf(modelo1, impulse = "Savings", response="Income", n.ahead = 20, boot = T )
plot(IncomeIRF, ylab = "Income", main = "Shock desde Consumptions") #le meti un shock al income para verificar la causalidad a largo plazo generada por consumo
ConsumptionIRF <- irf(modelo1, impulse = "Income", response="Savings", n.ahead = 20, boot = T )
plot(ConsumptionIRF, ylab = "Consumption", main = "Shock desde Income")
#Descomposición de la varianza
FEVD1 <- fevd(modelo1, n.ahead = 10)
plot(FEVD1) #parece como el consumo parece estar dando toda la varianza, es decir, quito el income para ver si es mejor modelo.
#lo podriamos volver un ARIMA debiso a esto
#prediccion
fore<-predict(modelo1, n.ahead = 10, ci=0.95)#nahead son los periodos a pronosticar, ci es nivel de confianza
fanchart(fore)
autoplot(forecast(modelo1))
