#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

library(tseries)

if( !require('rugarch') ){
   install.packages('rugarch', dependencies=TRUE, repos='http://cran.rstudio.com/')
}
library(rugarch)

if( !require('fGarch') ){
   install.packages('fGarch', dependencies=TRUE, repos='http://cran.rstudio.com/')
}
library(fGarch)

save_plots = FALSE

set.seed(0) 

# P 1: EPage 500
#
TbGdpPi = read.csv('../../BookCode/Data/TbGdpPi.csv', header=TRUE)
# r = the 91-day Treasury bill rate
# y = the log of real GDP
# pi = the inflation rate
#
TbGdpPi = ts(TbGdpPi, start=1955, freq=4)
Tbrate = TbGdpPi
Tbill = Tbrate[,1]
Del.Tbill = diff(Tbill)

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter14/chap_14_prob_1_tbill.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(1,2))
plot(Tbill)
acf(Tbill)
par(mfrow=c(1,1))
if( save_plots ){ dev.off() }

adf.test(Tbill)
kpss.test(Tbill)


if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter14/chap_14_prob_1_diff_tbill.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(1,2))
plot(Del.Tbill)
acf(Del.Tbill)
par(mfrow=c(1,1))
if( save_plots ){ dev.off() }

adf.test(Del.Tbill)
kpss.test(Del.Tbill)


# P 2: EPage 444
#
arma.garch.norm = ugarchspec(mean.model=list(armaOrder=c(1, 0)), variance.model=list(garchOrder=c(1, 1)))
Del.Tbill.arma.garch.norm = ugarchfit(data=Del.Tbill, spec=arma.garch.norm) # Del.Tbill is almost stationary (has GARCH effects)
show(Del.Tbill.arma.garch.norm)


# P 3: EPage 444
#
res = ts(residuals(Del.Tbill.arma.garch.norm, standardize=FALSE), start=1955, freq=4)
res_std = ts(residuals(Del.Tbill.arma.garch.norm, standardize=TRUE), start=1955, freq=4)
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter14/chap_14_prob_3_residuals.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(2,3))
plot(res, main='residuals')
acf(res, main='acf(residuals)')
acf(res^2, main='acf(residuals^2)')
plot(res_std, main='standardized residuals')
acf(res_std, main='acf(standardized residuals)')
acf(res_std^2, main='acf(standardized residuals^2)')
if( save_plots ){ dev.off() }


# P 4: EPage 444
#
Del.Log.Tbill = diff(log(Tbill))
Del.Log.Tbill.arma.garch.norm = ugarchfit(data=Del.Log.Tbill, spec=arma.garch.norm) # diff(log(Tbi11)) is more stationary
show(Del.Log.Tbill.arma.garch.norm)

res = ts(residuals(Del.Log.Tbill.arma.garch.norm, standardize=FALSE), start=1955, freq=4)
res_std = ts(residuals(Del.Log.Tbill.arma.garch.norm, standardize=TRUE), start=1955, freq=4)
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter14/chap_14_prob_4_residuals.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(2,3))
plot(res)
acf(res)
acf(res^2)
plot(res_std)
acf(res_std)
acf(res_std^2)
if( save_plots ){ dev.off() }


# (P 5-6): EPage 445:
#
GPRO = read.table('../../BookCode/Data/GPRO.csv')
garchm = ugarchspec(mean.model=list(armaOrder=c(0, 0), archm=TRUE, archpow=1), variance.model=list(garchOrder=c(1, 1)))
GPR0.garchm = ugarchfit(garchm, data=GPRO)
show(GPR0.garchm)


# P 7: EPage 446:
#
# r = the 91-day Treasury bill rate
# y = the log of real GDP
# pi = the inflation rate
#
TbGdpPi = read.csv('../../BookCode/Data/TbGdpPi.csv', header=TRUE)
TbGdpPi.diff = ts(apply(TbGdpPi[, -2], 2, diff), start=c(1955, 2), freq=4)
plot(TbGdpPi.diff)
acf(TbGdpPi.diff^2)

source('../../BookCode/SDAFE2.R')
mLjungBox(TbGdpPi.diff^2, lag=8)


# P 8: EPage 446:
#
EWMA.param = est.ewma(lambda.0=0.95, innov=TbGdpPi.diff)
print(EWMA.param$lambda.hat)
EWMA.Sigma = sigma.ewma(lambda=EWMA.param$lambda.hat, innov=TbGdpPi.diff)
par(mfrow = c(2,2))
plot(ts(EWMA.Sigma[1,1,]*.95, start = c(1955, 2), frequency = 4), type='l', xlab = 'year', ylab = NULL, main = expression(paste("(a) ", hat(sigma)["1,t"])))
plot(ts(EWMA.Sigma[1,2,], start = c(1955, 2), frequency = 4), type='l', xlab = "year", ylab = NULL, main = expression(paste("(b) ", hat(sigma)["12,t"])))
plot(ts(EWMA.Sigma[1,2,]/(sqrt(EWMA.Sigma[1,1,]* EWMA.Sigma[2,2,])), start = c(1955, 2), frequency = 4), type = 'l', xlab = 'year', ylab = NULL, main = expression(paste("(c) ", hat(rho)["12,t"])))
plot(ts(EWMA.Sigma[2,2,]*.5, start = c(1955, 2), frequency = 4), type ='l', xlab = "year", ylab = NULL, main = expression(paste("(d) ", hat(sigma)["2,t"])))
par(mfrow = c(1,1))


# P 9: EPage 446:
#
n = dim(TbGdpPi.diff)[1]
d = dim(TbGdpPi.diff)[2]
stdResid.EWMA = matrix(0, n, d)
for(t in 11:n){
    stdResid.EWMA[t,] = TbGdpPi.diff[t, ] %*% matrix.sqrt.inv(EWMA.Sigma[,,t])
}
mLjungBox(stdResid.EWMA^2, lag=8)


# P 10: EPage 446:
#
DOC.test(TbGdpPi.diff^2, 8)