setwd("D:/")
library(quantmod) 
library(fBasics)
library(sn)
library(PerformanceAnalytics)
library(car)
library(tseries)
library(forecast)
library(fGarch)
getSymbols("IPDCONGD",src="FRED")
dim(IPDCONGD)
tail(IPDCONGD)

getSymbols("IPNCONGD", src="FRED")
dim(IPNCONGD)

getSymbols("IPBUSEQ",src="FRED")
dim(IPBUSEQ)

getSymbols("IPMAT", src="FRED")
dim(IPMAT)

par(mfrow=c(2,1))
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")

par(mar=c(1,1,1,1))

require(MTS)
#Question 1
#QUestion 1.1 
IPA <- IP
IPA = log(IP,base=exp(10))
zt=diffM(IPA)*100
colnames(zt) <- c("IPD","IPN","IPB","IPM")
MTSplot(zt)

#Comment
#We could see that n_t (Non-durable goods) is more volatile compared to others.
#Also, we could see that d_t and b_t seems like to be correlated to each other,
#However, we are not sure whether they are correlated or not, we need to conduct
#statistical test.
#There are some extreme volatile in the series, especially the
#volatility during the pandemic

#Question 1.2
require(MTS)
#detach("package:vars")
VARorder(zt)
#BIC suggests that we should use lag-2, AIC suggests lag-6
m1 = VAR(zt,p=2)
MTSdiag(m1)
varfit=VAR(zt,p=2)
summary(varfit)
m2 = refVAR(m1,thres = 1.645)
#ModelChecking
MTSdiag(m2)

#Question 1.3

#Impulse response function of the fitted model
detach("package:MTS")
require(vars)
VARselect(zt,lag.max=15)
#SC suggests 2-lag model the same as BIC
m3=VAR(zt,p=2)
summary(m3)
impls=irf(m3)

plot(impls)
#Interpretation


#Question 1.4
#Forecast error variance decomposition
fevd(m3,n.ahead=1)
fevd(m3,n.ahead=6)
#Interpretation

#Question 1.5
#Compute 1-step to 6-step ahead origin h=
fm3 <- predict(m3, n.ahead=6, ci=.9)
fanchart(fm3)