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

# E 1 EPage 503
# 
library(Ecdat)
library(fGarch)
data(SP500,package="Ecdat")
returnBlMon = SP500$r500[1805] 
x = SP500$r500[(1804-2*253+1):1804]
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter18/chap_18_exercise_8_r_black_monday.eps", onefile=FALSE, horizontal=FALSE) }
plot(c(x,returnBlMon))
if( save_plots ){ dev.off() }

# Fit an AR(1)/ARCH(1,1) model:
# 
results = garchFit( ~arma(1,0) + garch(1,1), data=x, cond.dist="std" )
dfhat = as.numeric( results@fit$par[6] )
forecast = predict( results, n.ahead=1 )

# Part (a): 
z_score = ( returnBlMon - forecast$meanForecast ) / forecast$standardDeviation
pt( z_score, dfhat )

# Part (b): 
res = residuals(results)
res_std = res / results@sigma.t # standardized residuals
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter18/chap_18_exercise_8_standard_residuals.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(1,3))
plot(res_std)
acf(res_std)
acf(res_std^2)
par(mfrow=c(1,1))
if( save_plots ){ dev.off() }

# c: Look at fitting a AR(1)/ARCH(1) model:
#
results_c = garchFit( ~arma(1,0) + garch(1,0), data=x, cond.dist="std" )

res = residuals(results_c)
res_std = res / results_c@sigma.t # standardized residuals 
par(mfrow=c(1,3))
plot(res_std)
acf(res_std)
acf(res_std^2)
par(mfrow=c(1,1))

summary(results_c)

# d: Fit a AR(1) model and plot statistics:
#
results_d = arima( x, order=c(1,0,0) )

res = residuals(results_d)
res_std = res / sqrt(results_d$sigma2) # standardized residuals
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter18/chap_18_exercise_8_ar_1_plot.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(1,3))
plot(res_std)
acf(res_std)
acf(res_std^2)
par(mfrow=c(1,1))
if( save_plots ){ dev.off() }


# E 9 EPage 504
# 
library(fGarch)
library(Ecdat)
data(Irates)
r = as.numeric(log(Irates[,2]))
n = length(r)
lagr = r[1:(n-1)]
diffr = r[2:n] - lagr
model = garchFit( ~arma(1,0) + garch(1,1), data=diffr, cond.dist="std" )