Pick 5 years of monthly adjusted closing data for two stocks or funds that you might feel are reasonably correlated. Build a VAR model to investigate the relationships. Discuss.

#install.packages("urca")
#install.packages("vars")
library(forecast)
library(fpp2)
library(urca)
library(vars)
library(quantmod)
library(tseries)
library(ggplot2)
library(PerformanceAnalytics)
library(xts)
maxDate<- "2015-01-01"

amzn<- Ad(getSymbols("AMZN", auto.assign=F, from=maxDate))

sp<- Ad(getSymbols("^GSPC", auto.assign=F, from=maxDate))
mydata<- cbind(amzn,sp)

autoplot(mydata) + xlab("Time in days")

Autocorrelation: ACF & PACF

acf_a<- acf(amzn, main = "ACF for AMZN")

pacf_a<- pacf(amzn, main = "PACF for AMZN")

acf_s<- acf(sp, main = "S&P 500")

pacf_s<- pacf(sp, main = "S&P 500")

All of the lags of the ACF are significant, hence there is a lot of persistence.

Find optimal lags

lagselect<- VARselect(mydata, lag.max=10, type = "both")

lagselect$selection
## AIC(n)  HQ(n)  SC(n) FPE(n) 
##     10     10      2     10
# Lag of BIC = 2

Although AIC and HQ align on 10 lags, for VAR models the BIC is preferred.

Modeling

fit.mod<- VAR(mydata, p=2, type="both", season= NULL, exog= NULL)

summary(fit.mod)
## 
## VAR Estimation Results:
## ========================= 
## Endogenous variables: AMZN.Adjusted, GSPC.Adjusted 
## Deterministic variables: both 
## Sample size: 1408 
## Log Likelihood: -13165.701 
## Roots of the characteristic polynomial:
## 0.9977 0.9853 0.2098 0.004825
## Call:
## VAR(y = mydata, p = 2, type = "both", exogen = NULL)
## 
## 
## Estimation results for equation AMZN.Adjusted: 
## ============================================== 
## AMZN.Adjusted = AMZN.Adjusted.l1 + GSPC.Adjusted.l1 + AMZN.Adjusted.l2 + GSPC.Adjusted.l2 + const + trend 
## 
##                   Estimate Std. Error t value Pr(>|t|)    
## AMZN.Adjusted.l1  0.996417   0.033054  30.145  < 2e-16 ***
## GSPC.Adjusted.l1 -0.115621   0.030463  -3.795 0.000154 ***
## AMZN.Adjusted.l2  0.003376   0.033280   0.101 0.919223    
## GSPC.Adjusted.l2  0.104055   0.030504   3.411 0.000665 ***
## const            21.714709  10.483109   2.071 0.038504 *  
## trend             0.013829   0.006970   1.984 0.047423 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 27.99 on 1402 degrees of freedom
## Multiple R-Squared: 0.9981,  Adjusted R-squared: 0.9981 
## F-statistic: 1.491e+05 on 5 and 1402 DF,  p-value: < 2.2e-16 
## 
## 
## Estimation results for equation GSPC.Adjusted: 
## ============================================== 
## GSPC.Adjusted = AMZN.Adjusted.l1 + GSPC.Adjusted.l1 + AMZN.Adjusted.l2 + GSPC.Adjusted.l2 + const + trend 
## 
##                   Estimate Std. Error t value Pr(>|t|)    
## AMZN.Adjusted.l1 -0.012860   0.035250  -0.365 0.715293    
## GSPC.Adjusted.l1  0.781646   0.032487  24.061  < 2e-16 ***
## AMZN.Adjusted.l2  0.015996   0.035491   0.451 0.652275    
## GSPC.Adjusted.l2  0.198344   0.032530   6.097 1.39e-09 ***
## const            37.244270  11.179418   3.332 0.000886 ***
## trend             0.014629   0.007433   1.968 0.049232 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## 
## Residual standard error: 29.85 on 1402 degrees of freedom
## Multiple R-Squared: 0.9942,  Adjusted R-squared: 0.9942 
## F-statistic: 4.848e+04 on 5 and 1402 DF,  p-value: < 2.2e-16 
## 
## 
## 
## Covariance matrix of residuals:
##               AMZN.Adjusted GSPC.Adjusted
## AMZN.Adjusted         783.4         490.1
## GSPC.Adjusted         490.1         890.9
## 
## Correlation matrix of residuals:
##               AMZN.Adjusted GSPC.Adjusted
## AMZN.Adjusted        1.0000        0.5866
## GSPC.Adjusted        0.5866        1.0000

Model Diagnostics

Test for serial autocorrelation

serial1<- serial.test(fit.mod, lags.pt = 12)

serial1
## 
##  Portmanteau Test (asymptotic)
## 
## data:  Residuals of VAR object fit.mod
## Chi-squared = 243.03, df = 40, p-value < 2.2e-16

Portmanteau test is significant which means there is autocorrelation.

Test for heteroskedasticity

arch<- arch.test(fit.mod, lags.multi=12, multivariate.only = TRUE)

arch
## 
##  ARCH (multivariate)
## 
## data:  Residuals of VAR object fit.mod
## Chi-squared = 1213.3, df = 108, p-value < 2.2e-16

Normal distribution of Residuals

normality.test(fit.mod, multivariate.only = TRUE)
## $JB
## 
##  JB-Test (multivariate)
## 
## data:  Residuals of VAR object fit.mod
## Chi-squared = 24261, df = 4, p-value < 2.2e-16
## 
## 
## $Skewness
## 
##  Skewness only (multivariate)
## 
## data:  Residuals of VAR object fit.mod
## Chi-squared = 609.42, df = 2, p-value < 2.2e-16
## 
## 
## $Kurtosis
## 
##  Kurtosis only (multivariate)
## 
## data:  Residuals of VAR object fit.mod
## Chi-squared = 23652, df = 2, p-value < 2.2e-16

Forecasting

fcast<- predict(fit.mod, n.ahead = 30, ci=0.95)

fanchart(fcast, names = "AMZN.Adjusted")

fanchart(fcast, names = "GSPC.Adjusted")