#install.packages('quantmod')
#install.packages('MTs')
#install.packages('vars')
#install.packages('tseries')
#install.packages('forecast')
library(quantmod)
library(MTS)
library(tseries)
library(forecast)
library(fUnitRoots)
#1.1
getSymbols("IPDCONGD",src="FRED")
dim(IPDCONGD)
tail(IPDCONGD)
getSymbols("IPNCONGD",src="FRED")
getSymbols("IPBUSEQ",src="FRED")
getSymbols("IPMAT",src="FRED")
dim(IPNCONGD)
dim(IPBUSEQ)
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")
logIP=log(IP)
Zt=diff(logIP)*100
Zt
par('mar')
par(mar=c(1,1,1,1))
MTSplot(Zt)

#1.2
VARorder(Zt) #with lowest possible bic criteria
m.MTS=VAR(Zt,2)
MTSdiag(m.MTS)

m2.MTS=refVAR(m.MTS, thres=1.645)
MTSdiag(m2.MTS)

#1.3
detach(package:MTS, unload = TRUE)
require(vars)
VARselect(Zt, lag.max = 10)
m2.vars=VAR(Zt, p = 2)
summary(m2.vars)
IRF1=irf(m2.vars)
plot(IRF1)

#1.4
fevd(m2.vars, n.ahead=20)
m2vars.forecast <- predict(m2.vars, n.ahead = 6, ci=0.95)
m2vars.forecast

#2
da <- read.table("DOGE-USD.txt", fill=TRUE, header=TRUE)
D <- (da)
MOON=as.numeric(D[,6])
DD=log(MOON)
DDD=diff(DD)*100
Dt=na.omit(DDD)
D1 <- ar(Dt,order.max=10, na.action = na.pass)
D1$order
pacf(Dt, na.action = na.pass)

auto.arima(Dt)
D2 <- arima(Dt, order=c(3,0,3))
summary(D2)
#checking for ARCH effect in the model
#install.packages('fGarch')
library(fGarch)
acf(Dt)
t.test(Dt)
Box.test(Dt,lag=10,type='Ljung')
Box.test(Dt^2,lag=10,type='Ljung')
#Rejecting null hull hypothesis at 95% confidence interval, the series contains ARCH effect
G1=garchFit(~arma(3,3)+garch(1,1),data=Dt,trace=F)
summary(G1)# mean model is insignificance
G2=garchFit(~arma(3,3)+garch(1,1),data=Dt,cond.dist='std',trace=F)
summary(G2)
G3=garchFit(~arma(2,2)+garch(1,1),data=Dt,cond.dist='std',trace=F)
summary(G3)
G4=garchFit(~arma(1,1)+garch(1,1),data=Dt,cond.dist='std',trace=F)
summary(G4)
G5=garchFit(~garch(1,1),data=Dt,cond.dist='std',trace=F)
summary(G5)
G6=garchFit(~garch(2,2),data=Dt,cond.dist='std',trace=F)
summary(G6)

#3.1
require(MTS)
US=read.table("USGDP.txt", header=T)
CA=read.table("CAGDP.txt", header=T)
UK=read.table("UKGDP.txt", header=T)
X=cbind(as.numeric(US$GDPC1), as.numeric(CA$NAEXKP01CAQ189S), as.numeric(UK$CLVMNACSCAB1GQUK))
Z=log(X)
Zt=diffM(Z)*100
colnames(Zt)<-c('US' , 'CA', 'UK')

#finding optimal lag for the model with lowest possible BIC criteria
#simple VAR model
VARorder(Zt)
G.MTS=VAR(Zt,1)
MTSdiag(G.MTS)
#refind VAR model with threshold equals to 1.645
G2.MTS=refVAR(G.MTS, thres=1.645)
MTSdiag(G2.MTS)

detach(package:MTS, unload=TRUE)
library(vars)
fittedmodel=VAR(Zt,p=1)
summary(fittedmodel)
#The fitted model is appropriate, since the value of the roots are lower than 1, the model construct the property of weak stationarity

#3.2
IRF2=irf(fittedmodel)
plot(IRF2)


#3.3
fevd(fittedmodel, n.ahead = 20)

#3.4
#install.packages('fUnitRoots')
library(quantmod)
USTEST=cbind(as.numeric(US$GDPC1))
UKTEST=cbind(as.numeric(UK$CLVMNACSCAB1GQUK))
CATEST=cbind(as.numeric(CA$NAEXKP01CAQ189S))
fitted <- lm(USTEST ~ UKTEST + CATEST)
summary(fitted)

#engle-granger cointegration test

adf.test(UKTEST)
adf.test(diff(UKTEST))

adf.test(CATEST)
adf.test(diff(CATEST))

adf.test(USTEST)
adf.test(diff(USTEST))

error = residuals(fitted)
adf.test(error)
#reject null hypothesis if the p-value is less than critical value, indicate no unit root problem

#3.5
DUK=diff(UKTEST)
DUS=diff(USTEST)
DCA=diff(CATEST)
DUKL1=Lag(DUK, k=1)
DUSL1=Lag(DUS, k=1)
DCAL1=Lag(DCA, k=1)
EL1=Lag(error, k=1)
EL1=EL1[2:126]
VECM1 <- lm(DUK ~ DUKL1+DCAL1+DUSL1+EL1)
VECM2 <- lm(DUS ~ DUKL1+DCAL1+DUSL1+EL1)
VECM3 <- lm(DCA ~ DUKL1+DCAL1+DUSL1+EL1)
summary(VECM1)
summary(VECM2)
summary(VECM3)