names = c("Durable consumer goods","Nondurable consumer goods","Business equivalent","Materials")
data = read.table("http://faculty.chicagobooth.edu/ruey.tsay/teaching/mts/sp2015/m-ip4comp.txt",header=TRUE)
n = nrow(data)
date = data[,1]+(data[,2]-1)/12
date1=date[2:n]


pdf(file="usipi-tsplots.pdf",width=10,height=8)
par(mfrow=c(2,2))
plot(date,log(data[,3]),main=names[1],type="l",ylab="Log",xlab="")
plot(date,log(data[,4]),main=names[2],type="l",ylab="Log",xlab="")
plot(date,log(data[,5]),main=names[3],type="l",ylab="Log",xlab="")
plot(date,log(data[,6]),main=names[4],type="l",ylab="Log",xlab="")
dev.off()

pdf(file="usipi-acfplots.pdf",width=10,height=8)
par(mfrow=c(2,2))
acf(log(data[,3]),main="",lag.max=50)
acf(log(data[,4]),main="",lag.max=50)
acf(log(data[,5]),main="",lag.max=50)
acf(log(data[,6]),main="",lag.max=50)
dev.off()


n = nrow(data)
t0=756
y   = log(data[,3:6])
y1 = log(data[1:t0,3:6])
y2 = log(data[(t0+1):n,3:6])
n1 = nrow(y1)
coefs = matrix(0,4,2)
sigs = rep(0,4)
for (i in 1:4){
  reg = lm(y1[2:n1,i]~y1[1:(n1-1),i])
  coefs[i,] = reg$coef
  sigs[i] = summary(reg)$sigma
}
cbind(coefs,sigs)

h = 36
hs = 0:(h-1)
forecast = matrix(0,h,4)
sd.fore = matrix(0,h,4)
L = matrix(0,h,4)
U = matrix(0,h,4)
for (i in 1:4){
  forecast[,i]=coefs[i,1]*cumsum(coefs[i,2]^hs)+coefs[i,2]^hs*y[n1,i]
  sd.fore[,i] = sigs[i]*sqrt(cumsum(coefs[i,2]^(2*hs)))
  U[,i] = forecast[,i]+2*sd.fore[,i]
  L[,i]  = forecast[,i]-2*sd.fore[,i]
}

pdf(file="usipi-forecast1.pdf",width=10,height=8)
par(mfrow=c(2,2))
for (i in 1:4){
  plot(date,y[,i],main=names[i],type="l",ylab="Log",xlab="",ylim=range(y[,i],U[,i],L[,i]))
  lines(date[n1]+(hs+1)/12,forecast[,i],col=4,lwd=2)
  lines(date[n1]+(hs+1)/12,L[,i],col=4,lty=3,lwd=2)
  lines(date[n1]+(hs+1)/12,U[,i],col=4,lty=3,lwd=2)
  abline(v=date[n1],lty=2)
}
dev.off()

pdf(file="usipi-forecast2.pdf",width=10,height=8)
par(mfrow=c(2,2))
for (i in 1:4){
  ind = (n1-120):n    
  plot(date[ind],y[ind,i],main=names[i],type="l",ylab="Log",xlab="",
          xlim=c(date[ind[1]],date[n]),ylim=range(y[ind,i],U[,i],L[,i]))
  lines(date[n1]+(hs+1)/12,forecast[,i],col=4,lwd=2)
  lines(date[n1]+(hs+1)/12,L[,i],col=4,lty=3,lwd=2)
  lines(date[n1]+(hs+1)/12,U[,i],col=4,lty=3,lwd=2)
  abline(v=date[n1],lty=2)
}
dev.off()


# Verifying possible cointegration between Durable consumer goods and Materials
y1 = log(data[,3])
y2 = log(data[,6])
y1 = (y1-mean(y1))/sqrt(var(y1))
y2 = (y2-mean(y2))/sqrt(var(y2))
reg = lm(y1~y2)

pdf(file="usipi-cointegration.pdf",width=10,height=6)
par(mfrow=c(1,2))
plot(date,y1,main="",type="l",ylab="Log (standardized)",xlab="",ylim=range(y1,y2))
lines(date,y2,col=2)
legend("topleft",legend=c(names[1],names[4]),col=1:2,lty=2)
acf(reg$res,main="ACF of residuals",lag=50)
dev.off()

install.packages("tseries")
library(tseries)
adf.test(reg$res)