library(quantmod)
library(vars)
library(tseries)

# Download data from FED's FRED database

# FED's total assets, weekly
getSymbols('WALCL',src='FRED')
plot(WALCL)
start = start(WALCL)
end = end(WALCL)

dfed <- diff(log(WALCL))*100
plot(dfed)
head(dfed)
dfed <- as.double(dfed[-1])
class(dfed)

# Download S&P 500 from Yahoo Finance

getSymbols.yahoo("^GSPC", env=globalenv(), periodicity = "weekly", from = start, to = end)
head(GSPC)
plot(GSPC[,6])


rsp500 <- diff(log(GSPC[,6]))*100
plot(rsp500)
head(rsp500)
rsp500 <- as.double(rsp500[-1])

m <- lm(rsp500~dfed)
summary(m)

dat <- cbind(rsp500,dfed)
names(dat) <- c("rsp500","dfed")
head(dat)

# Select lag length

VARselect(dat, lag.max = 20)

var <- VAR(dat,p=4)
summary(var)

# Granger causality

causality(var,cause='dfed')

causality(var,cause='rsp500')

graphics.off()

# Impulse-response functions

ir = irf(var, n.ahead = 20,ci = 0.95, cumulative=T)
plot(ir)

# Forecast error variance decomposition

vd = fevd(var,n.ahead = 20)
plot(vd)
vd