#EE435 final take-home exam by jaruvit chantana 6104641482 seat 16
setwd("/Users/Phumj/Documents/EE435/6104641482_exam")
cat(rep("\n",50))
library(quantmod)
library(fBasics)
library(sn)
library(PerformanceAnalytics)
library(car)
library(tseries)
library(forecast)


#question 1
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)
IP = cbind(as.numeric(IPDCONGD),as.numeric(IPNCONGD),as.numeric(IPBUSEQ),as.numeric(IPMAT[-c(1:96)]))
IP = log(IP)           
zt = diff(IP)*100
colnames(zt) <- c("IPDCONGD","IPNCONGD","IPBUSEQ","IPMAT")
library(MTS)
require(MTS)

par('mar')
par(mar=c(1,1,1,1))
#question 1.1
MTSplot(zt)


#question 1.2
ccm(zt)
VARorder(zt)
# I choose BIC
m1.MTS=VAR(zt,2)
MTSdiag(m1.MTS)
#looking for p-value the model is adequate
#refine the model by removing insignificant parameters with treshold of t-ratio 1.645
m2.MTS=refVAR(m1.MTS, thres = 1.645)
MTSdiag(m2.MTS)


#question 1.3
detach_package <- function(pkg, character.only = FALSE)
{
  if(!character.only)
  {
    pkg <- deparse(substitute(pkg))
  }
  search_item <- paste("package", pkg, sep = ":")
  while(search_item %in% search())
  {
    detach(search_item, unload = TRUE, character.only = TRUE)
  }
}
detach_package("MTS", TRUE)
library(vars)
require(vars)
VARselect(zt,lag.max = 10)
m1.VAR=VAR(zt,p=6)
summary(m1.VAR)
imp1=irf(m1.VAR)
plot(imp1)


#question 1.4
fevd(m1.VAR, n.ahead = 20)


#question 1.5
m1.VAR.prd <- predict(m1.VAR, n.ahead = 6, ci= 0.95)
fanchart(m1.VAR.prd)
m1.VAR.prd



#question 2 using dogecoin-USD series
getSymbols("DOGE-USD", src = "yahoo",from="2015-01-01",to="2021-05-24")
da=`DOGE-USD`
plot(da)
dcrtn <- diff(log(da[,6]))*100
dcrtn=na.omit(dcrtn)
ts.plot(dcrtn)
library(fUnitRoots)
adfTest(dcrtn)
require(quantmod)
#trying to find the optimal lag
m1 <- ar(dcrtn,order.max=10) 
m1$order
pacf(dcrtn)
m2 <- arima(dcrtn, order=c(6,0,0))
m2
#doing refinement process by dividing ar-value by its s.e. it is 0 if the consequences is less than 1.96, if not, it is considered to be NA.
c1 <- c(0,NA,NA,0,0,NA,0)
m3 <- arima(dcrtn,order=c(6,0,0),fixed=c1)
m3
require(forecast)
auto.arima(dcrtn)
#when we know the lag given by autoarima then we run the code by ourselves
m4 <- arima(dcrtn,order=c(3,0,3))
m4
#doing back test to check whether m3 or m4 is better
######################################################################################################
"backtest" <- function(m1,rt,orig,h=1,xre=NULL,fixed=NULL,include.mean=TRUE){
  # m1: is a time-series model object
  # orig: is the starting forecast origin
  # rt: the time series
  # xre: the independent variables
  # h: forecast horizon
  # fixed: parameter constraint
  # inc.mean: flag for constant term of the model.
  #
  regor=c(m1$arma[1],m1$arma[6],m1$arma[2])
  seaor=list(order=c(m1$arma[3],m1$arma[7],m1$arma[4]),period=m1$arma[5])
  nT=length(rt)
  if(orig > nT)orig=nT
  if(h < 1) h=1
  rmse=rep(0,h)
  mabso=rep(0,h)
  nori=nT-orig
  err=matrix(0,nori,h)
  jlast=nT-1
  if(!is.null(xre))xre <- matrix(xre)
  for (n in orig:jlast){
    jcnt=n-orig+1
    x=rt[1:n]
    if (is.null(xre))
      pretor=NULL else pretor=xre[1:n,]
    mm=arima(x,order=regor,seasonal=seaor,xreg=pretor,fixed=fixed,include.mean=include.mean)
    if (is.null(xre)){nx=NULL} 
    else {nx=matrix(xre[(n+1):(n+h),],h,ncol(xre))}
    fore=predict(mm,h,newxreg=nx)
    kk=min(nT,(n+h))
    # nof is the effective number of forecast at the forecast origin n.
    nof=kk-n
    pred=fore$pred[1:nof]
    obsd=rt[(n+1):kk]
    err[jcnt,1:nof]=obsd-pred
  }
  #
  for (i in 1:h){
    iend=nori-i+1
    tmp=err[1:iend,i]
    mabso[i]=sum(abs(tmp))/iend
    rmse[i]=sqrt(sum(tmp^2)/iend)
  }
  print("RMSE of out-of-sample forecasts")
  print(rmse)
  print("Mean absolute error of out-of-sample forecasts")
  print(mabso)
  backtest <- list(origin=orig,error=err,rmse=rmse,mabso=mabso)
}
######################################################################################################
backtest(m3,dcrtn,2300,fixed=c1)
backtest(m4,dcrtn,2300)
#observing in RMSE and MAE each has one that higher and lower comparing to another
#so selecting from information criteria, m4 or ARMA(2,2) provides less value in AIC
#so I picked m4 -> ARMA(3,3) as my mean equation


#next volatility model
library(fGarch)
acf(m4$residuals^2)
Box.test(m4$residuals^2,lag=10,type = 'Ljung')
#reject H0 there is the ARCH effect with 95% CI
#finding vola equation
g1 <- garchFit(~arma(3,3)+garch(1,1),data=dcrtn,trace=F)
summary(g1)
plot(g1,which=13)

g2 <- garchFit(~arma(3,3)+garch(1,1),data=dcrtn,trace=F,cond.dist = "std")
summary(g2)
plot(g2,which=13)
#model is not adequate and not distributed as normal
#finding the new way

g3= garchFit(~arma(2,2)+garch(1,1),data=dcrtn,cond.dist="std",trace=F)
summary(g3)
plot(g3,which=13)
#model is not adequate and there are many insignificant values
#finding the new way by cut mean equation out due to insignificance
g4=garchFit(~garch(1,1),data = dcrtn,trace=F,cond.dist = "std")
summary(g4)
plot(g4,which=13)
#adequate at 99% CI

#question 3
GDPUS=read.table("US.txt", header = T)
GDPUK=read.table("UK.txt", header = T)
GDPCD=read.table("CD.txt", header = T)
head(GDPUS)
head(GDPUK)
head(GDPCD)
x=cbind(as.numeric(GDPUK$CLVMNACSCAB1GQUK),as.numeric(GDPCD$NAEXKP01CAQ189S),as.numeric(GDPUS$GDPC1))
x=log(x)
detach_package("vars", TRUE)
require(MTS)
zt=diffM(x)*100
colnames(zt) <- c('gdpukg' , 'gdpusg' , 'gdpcdg')


#question 3.1
VARorder(zt)
m3.1=VAR(zt,1)
MTSdiag(m3.1)
#seems fine looking at Ljung-box Statistics
m3.1ref=refVAR(m3.1, thres = 1.645)
#final fitted model


#question 3.2
detach_package("MTS", TRUE)
require(vars)
varfit=VAR(zt,p=1)
summary(varfit)
#roots are less than 1 meaning by using this model, it satisfies the weak stationarity
imp3=irf(varfit)
plot(imp3)


#question 3.3
fevd(varfit, n.ahead = 20)
#their own shocks impact themselves the most


#question 3.4
#install.packages("fUnitRoots")
da341=GDPUK$CLVMNACSCAB1GQUK
da342=GDPCD$NAEXKP01CAQ189S
da343=GDPUS$GDPC1
m341=adfTest(da341, lags = 1, type =c("ct"),title = NULL,description = NULL)
m341@test$p.value
m341@test$parameter
m341@test$lm
# 0.5466456 > 0.05 cannot reject that there is unit root in the series
m342=adfTest(da342, lags = 1, type =c("ct"),title = NULL,description = NULL)
m342@test$p.value
m342@test$parameter
m342@test$lm
#there is unit root
m343=adfTest(da343, lags = 1, type =c("ct"),title = NULL,description = NULL)
m343@test$p.value
m343@test$parameter
m343@test$lm
#there is unit root
#so take log into it
m341diff=adfTest(diff(da341) , lags = 1, type=c("ct"),title = NULL, description = NULL)
m341diff@test$p.value
#no unit root now => stationarity => change in GDPUK --> GDPUK~I(1)
m342diff=adfTest(diff(da342) , lags = 1, type=c("ct"),title = NULL, description = NULL)
m342diff@test$p.value
#no unit root now => stationarity => change in GDPCD --> GDPCD~I(1)
m343diff=adfTest(diff(da343) , lags = 1, type=c("ct"),title = NULL, description = NULL)
m343diff@test$p.value
#no unit root now => stationarity => change in GDPUS --> GDPUS~I(1)
fit <- lm(da341 ~ da342 + da343)
summary(fit)
error=residuals(fit)
merror=adfTest(error, lag = 1,type = c("nc"),title = NULL, description = NULL)
merror@test$p.value
#reject H0 -> no unit root


#question 3.5
diff.UK= diff(da341)
diff.CD= diff(da342)
diff.US= diff(da343)
diff.UK.L.1=Lag(diff(da341),k=1)
diff.CD.L.1=Lag(diff(da342),k=1)
diff.Us.L.1=Lag(diff(da343),k=1)
error.L.1= Lag(error, k=1)
error.L.1= error.L.1[2:126]
#[2:126] due to the first difference
fit1 <- lm(diff.UK ~ diff.UK.L.1+diff.CD.L.1+diff.Us.L.1+error.L.1)
summary(fit1)
fit2 <- lm(diff.CD ~ diff.UK.L.1+diff.CD.L.1+diff.Us.L.1+error.L.1)
summary(fit2)
fit3 <- lm(diff.US ~ diff.UK.L.1+diff.CD.L.1+diff.Us.L.1+error.L.1)
summary(fit3)