#
# 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.
# 
#-----

save_plots = F

set.seed(0) 

# P 1: EPage 500
#
data(Tbrate,package="Ecdat")
library(tseries)
library(fGarch)
# r = the 91-day treasure bill rate
# y = the log of real GDP
# pi = the inflation rate
Tbill = Tbrate[,1]
Del.Tbill = diff(Tbill)

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter18/chap_18_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/Chapter18/chap_18_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 500
#
#garch.model.Tbill = garchFit(formula = ~arma(1,0) + garch(1,0), Tbill) # Tbill is not stationary 
garch.model.Tbill = garchFit(formula = ~arma(1,0) + garch(1,0), Del.Tbill) # Del.Tbill is almost stationary (has GARCH effects)
summary(garch.model.Tbill)
garch.model.Tbill@fit$matcoef


# P 3: EPage 500
#
res = residuals(garch.model.Tbill)
res_std = res / garch.model.Tbill@sigma.t # standardized residuals 
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter18/chap_18_prob_3_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 3: EPage 500
#
Del.Log.Tbill = diff(log(Tbill))
garch.model.Tbill = garchFit(formula = ~arma(1,0) + garch(1,0), data=Del.Log.Tbill)

res = residuals(garch.model.Tbill)
res_std = res / garch.model.Tbill@sigma.t # standardized residuals 
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter18/chap_18_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() }