setwd("D:/uni year 3/EE435")
cat(rep("\n", 50)) #clear R console
#install.packages("quantmod")  
#install.packages("fBasics") 
#install.packages("sn")  
#install.packages("PerformanceAnalytics") 
#install.packages("car") 
#install.packages("tseries")  
#install.packages("forecast") 
#install.packages("fGarch")
#install.packages("vars")
#install.packages("fUnitRoots")
library(fGarch)
library(quantmod) 
library(fBasics)
library(sn)
library(PerformanceAnalytics)
library(car)
library(tseries)
library(forecast)
require(quantmod)
require(vars)
library(fUnitRoots)

##Question 1
#Load data from FRED
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)
#Combine data
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")
zt = diff(log(IP))
#1.1 Plot time series data
plot.ts(zt)
#1.2 VAR model
m1 = VAR(zt, p = 2)
summary(m1)

#1.3 Impulse response function
impresp_1 = irf(m1)
plot(impresp_1)
#1.4 Forecast error variance decomposition
fevd(m1, n.ahead = 20)
#1.5 Prediction
m1.prd <- predict(m1, n.ahead = 6, ci = 0.95)
m1.prd
fanchart(m1.prd)

#Question 2
getSymbols("DOGE-USD",src="yahoo",from = "2015-01-01", to = "2021-05-24")
da1 = na.omit(`DOGE-USD`)
doge_price = da1[,6]
plot(doge_price)
log_dogeprice = log(doge_price)
logrt_doge = diff(log_dogeprice)
nlogrt_doge = logrt_doge[2:nrow(logrt_doge),]
plot(nlogrt_doge)
chart.Histogram(nlogrt_doge, methods = c("add.normal"))
acf(nlogrt_doge)
pacf(nlogrt_doge)
pacf(nlogrt_doge^2)
Box.test(nlogrt_doge^2, lag = 10, type = "Ljung")
table.Stats(nlogrt_doge)
t.test(nlogrt_doge)
Box.test(nlogrt_doge, lag = 15, type = "Ljung")

m2 = garchFit(~arma(1,1)+garch(1,1),data = nlogrt_doge, trace = F)
summary(m2)
m3 = garchFit(~garch(1,1),data = nlogrt_doge, trace = F)
summary(m3)
m4 = garchFit(~garch(3,3), data = nlogrt_doge, trace = F)
summary(m4)
m5 = garchFit(~arma(6,6)+garch(1,3), data = nlogrt_doge, trace = F)
summary(m5)

#Question 3
getSymbols("CLVMNACSCAB1GQUK", src = "FRED")
uk_gdp = CLVMNACSCAB1GQUK[c(21:146)]
getSymbols("NAEXKP01CAQ189S", src = "FRED")
can_gdp = NAEXKP01CAQ189S[c(77:202)]
getSymbols("GDPC1", src = "FRED")
us_gdp = GDPC1[c(133:258)]
gdp = cbind(as.numeric(uk_gdp), as.numeric(can_gdp), as.numeric(us_gdp))
colnames(gdp) <- c("uk_gdp", "can_gdp", "us_gdp")
head(gdp)
zt2 = diff(log(gdp))
head(zt2)
#3.1
varfit = VAR(zt2, p = 1)
summary(varfit)
#3.2
impresp_2 = irf(varfit)
plot(impresp_2)
#3.3
fevd(varfit)
#3.4
m6 = adfTest(uk_gdp, lag = 3, type = c("ct"))
m6
m7= adfTest(can_gdp, lag = 3, type = c("ct"))
m7
m8 = adfTest(us_gdp, lag = 3, type = c("ct"))
m8
uk2 = diff(uk_gdp)
uk2 = uk2[c(2:nrow(uk_gdp))]
m9 = adfTest(uk2, lag = 3, type = c("ct"))
m9
can2 = diff(can_gdp)
can2 = can2[c(2:nrow(can_gdp))]
m10 = adfTest(can2, lag = 3, type = c("ct"))
m10
us2 = diff(us_gdp)
us2 = us2[c(2:nrow(us_gdp))]
m11 = adfTest(us2, lag = 3, type = c("ct"))
m11
#uk canada us I(1)
fit <- lm(uk_gdp~can_gdp+us_gdp)
summary(fit)
error = residuals(fit)
m12 = adfTest(error, lags = 3 , type = c("ct"))
m12

#3.5 VECM
uk.l1 = Lag(uk2, k=1)
can.l1 = Lag(can2, k=1)
us.l1 = Lag(us2, k=1)
error.l1 = Lag(error, k=1)
error.l1 = error.l1[2:126]
fit2 <- lm(uk2~uk.l1+can.l1+us.l1+error.l1)
summary(fit2)
m12 = lm(can2~uk2+us2)
fit3<- lm(can2~uk.l1+can.l1+us.l1+error.l1)
summary(fit3)
fit4 <- lm(us2~uk.l1+can.l1+us.l1+error.l1)
summary(fit4)