#************************************************************ 
# 1. ARIMA for 3 Month T-Bill  
#install.packages("nortest")
# install.packages("TTR")
# install.packages("forecast")
#------------------------------------------------------------

library("TTR")
library("forecast")
library(quantmod)

#############################################################################
# Here is a funciton to compare manual vs. auto.arima estimation
# It also perform out sample forecast
# The data is divided by train set and test set
#----------------------------------------------------------------------#
test_arima <- function(n,arma,h){
  par(mfrow = c(2, 2))
  ts.plot(y, main = "ARMA Process ", col = "red")
  acf(y)
  pacf(y)
  
  # ARIMA Estimation
  p1 <- arima(y, order = c(1, 0, 0))
  summary(p1)
  plot(forecast(p1, h = h), include = 2*h)

  # Residual CHeck
  par(mfrow = c(2, 2))
  plot.ts(p1$residuals)
  Box.test(p1$residuals, lag = 20, type = "Ljung-Box")
  jarque.bera.test(p1$residuals)
  acf(p1$residuals)
  pacf(p1$residuals)
  
  # Train and Check out of sample forecast
  par(mfrow = c(1, 1))
  ts.plot(y, y.train, y.test, col = c("black", "blue", "red"))
  
  par(mfrow = c(2, 1))

  p1 <- arima(y.train, order = arma)
  p2 <- auto.arima(y.train, ic = "bic", trace = TRUE)
  summary(p1)
  summary(p2)
  
  p1.forc = forecast(p1, h = h)
  p2.forc = forecast(p2, h = h)
  plot(p1.forc)
  plot(p2.forc)
  accuracy(p1.forc, y.test)
  accuracy(p2.forc, y.test)
  
  plotarimapred(y.test, p1, xlim = c(n-3*h,n))
  plotarimapred(y.test, p2, xlim = c(n-3*h,n))
  
}
##########################################################################

title="3-Month Treasury Bill: Secondary Market Rate (TB3MS)"
getSymbols("TB3MS", src = "FRED")

startyear = 1934
endyear   = 2017
yraw <- ts(TB3MS,start=c(1934,1),end=c(2007,12),frequency = 12)
plot.ts(yraw)
adf.test(yraw)

y <- diff(ts(yraw, start=1, frequency=1))
n <- nrow(y) 
h <- 36
arma <- c(1, 0, 1)  # Find your ARMA terms based on the ACF and PACF 
y.train <- window(y, frequency = 1, start = 1  , end=n-h )
y.test  <- window(y, frequency = 1, start = n-h+1, end=n ) 

test_arima(n,arma,h)


# ARFIMA
fit1 <- arfima(y, order = c(0, 0, 1), cpus = 2, back=T)
fit1
AIC(fit1)