library(fBasics)
library(timeDate)
library(timeSeries)
library(fGarch)
library(quantmod)
library(forecast)

getSymbols("IPDCONGD",src="FRED")
dim(IPDCONGD)
getSymbols("IPNCONGD",src="FRED")
dim(IPNCONGD)
getSymbols("IPBUSEQ",src="FRED")
dim(IPBUSEQ)
getSymbols("IPMAT",src="FRED")
dim(IPMAT)

IP = cbind(as.numeric(IPDCONGD), as.numeric(IPNCONGD), as.numeric(IPBUSEQ), as.numeric(IPMAT[-c(1:96)]))
dim(IP)
colnames(IP) <- c("IPD","IPN","IPB","IPM")

require(MTS)

#Question 1.1
IP=log(IP)
zt=diffM(IP)
tdx=c(1:891)/12+1947
MTSplot(zt,tdx)

#Question 1.2
VARorder(zt)
m1 <- VAR(zt, p=6)
m2 <- refVAR(m1, thres=1.645)
MTSdiag(m2)

#Question 1.3
detach("package:MTS", unload = TRUE)
require(vars)
varfit1 <- VAR(zt,p=6)
summary(varfit1)
varfit2 <- restrict(varfit1, thresh=1.645)
summary(varfit2)
impresp <- irf(varfit2)
plot(impresp)

#Question 1.4
fevd(varfit2, start = 879, n.ahead=6)

#Question 1.5
detach("package:vars", unload = TRUE)
require(MTS)
m1 <- VAR(zt, p=6)
m2 <- refVAR(m1, thres=1.645)
VARpred(m2, h = 6, orig = 879)
rm(list=ls())