# S&P 500 data from Robert Shiller's web page

shiller <- read.csv("D:/data/shiller.csv")
head(shiller)
# P = stock (S&P 500 index) price
# D = stock dividend
# E = stock earnings

dat <- ts(shiller[,2:4], start=c(1871,1), frequency=12)
dim(dat)
dim(dat)[1]/12 # Number of years
dum <- rep(c(0,0,1,0,0,1,0,0,1,0,0,1),150)
length(dum)

dat <- cbind(dum,shiller)
head(dat)
colnames(dat) <- c("dum","month","P","D","E")
head(dat)

qdata <- dat[dum==1,c(2:5)]
head(qdata)

# Log differences
dp <- diff(log(qdata$P))*100 # Log return
dd <- diff(log(qdata$D))*100 # Log dividend growth
de <- diff(log(qdata$E))*100 # Log earnings growth

library(vars) # VAR model

dat0 <- cbind(dp,dd,de)

VARselect(dat0,lag.max=8,type="const")

var <- VAR(dat0, p=6, type=c("const"))
summary(var)

causality(var, cause=c("dd","de"))

causality(var, cause=c("dp","de"))

causality(var, cause=c("dp","dd"))

ir = irf(var, n.ahead = 24, cumulative=FALSE)
plot(ir)

ir = irf(var, n.ahead = 24, cumulative=TRUE)
plot(ir)

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

### Cointegration analysis

lp <- log(qdata$P) # log price
ld <- log(qdata$D) # log divideds
le <- log(qdata$E) # log earnings

library(fUnitRoots)

adfTest(lp, lags=4, type="c")
adfTest(ld, lags=4, type="c")
adfTest(le, lags=4, type="c")

dat1 <- cbind(lp,ld,le)

VARselect(dat1,lag.max=24,type="c")

summary(ca.jo(dat1, type="eigen", K=3, ecdet="const"))

summary(ca.jo(dat1, type="trace", K=3, ecdet="const"))

# VECM
vecm <- ca.jo(dat1, type="eigen", K=3, ecdet="const")

cajorls(vec, r=2)