#
# 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 = FALSE

set.seed(0)

# E 1 EPage 77
# 
ford.s = read.csv('../../BookCode/Data/ford.csv', header=TRUE)
ford.s = ford.s[, c(2, 3)]

summary(ford.s[,2])
sd(ford.s[,2])

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter4/chap_4_ex_1_normal_qqplot.eps", onefile=FALSE, horizontal=FALSE) }
qqnorm(ford.s[,2],datax=T)
qqline(ford.s[,2],datax=T)
grid()
if( save_plots ){ dev.off() }

# Part c:
#
shapiro.test(ford.s[,2])

# Part d (following the code in the Rlab):
#
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter4/chap_4_ex_1_t_dof_plots.eps", onefile=FALSE, horizontal=FALSE) }
logR = ford.s[,2]
#logR = logR[ -which.min(logR) ] # Remove the black Monday return
n = length(logR)
q.grid = (1:n)/(n+1)
df = c(1,4,6,10,20,30)
par(mfrow=c(3,2))
for(j in 1:6){
  qqplot(logR, qt(q.grid,df=df[j]), main=paste("F, df=", df[j]))
  abline(lm( qt(c(0.25,0.75),df=df[j]) ~ quantile(logR,c(0.25,0.75)) ))
}
if( save_plots ){ dev.off() }

# Part e:
#
F_inv_q = median( logR )
d = density( logR )
a = approx(d$x, d$y, F_inv_q, method = "linear")
q = 0.5

sqrt( ( q*(1-q) ) / ( length(logR) * a$y^2 ) ) 

sd(logR)/sqrt(length(logR))


# E 2:
# 
ford.s = read.csv('../../BookCode/Data/RecentFord.csv', header=TRUE)
logR = diff(log(ford.s[, 7]))

summary(logR)
sd(logR)

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter4/chap_4_ex_2_normal_qqplot.eps", onefile=FALSE, horizontal=FALSE) }
qqnorm(logR, datax=T)
qqline(logR, datax=T)
grid()
if( save_plots ){ dev.off() }

# Part c:
#
shapiro.test(logR)

# Part d (following the code in the Rlab):
#
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter4/chap_4_ex_2_t_dof_plots.eps", onefile=FALSE, horizontal=FALSE) }
#logR = logR[ -which.min(logR) ] # Remove the black Monday return
n = length(logR)
q.grid = (1:n)/(n+1)
df = c(1,4,6,10,20,30)
par(mfrow=c(3,2))
for(j in 1:6){
  qqplot(logR, qt(q.grid,df=df[j]), main=paste("F, df=", df[j]))
  abline(lm( qt(c(0.25,0.75),df=df[j]) ~ quantile(logR,c(0.25,0.75)) ))
}
if( save_plots ){ dev.off() }

# Part e:
#
F_inv_q = median( logR )
n = length(logR)
d = density( logR )
a = approx(d$x, d$y, F_inv_q, method = "linear")
q = 0.5

sqrt( ( q*(1-q) ) / ( n * a$y^2 ) ) 

sd(logR)/sqrt(n)


# E 3 EPage 77
# 
if( "Ecdat" %in% rownames(installed.packages()) == FALSE ){
    install.packages("Ecdat")
}
library(Ecdat)
diff_dy = diff(Garch$dy)

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter4/chap_4_ex_2_density_plot.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(1,1))
x = seq(min(diff_dy),max(diff_dy),length.out=100)
plot(density(diff_dy), lwd=2)
lines(x,dnorm(x,mean=mean(diff_dy),sd=sd(diff_dy)),lty=5,lwd=2)
lines(x,dnorm(x,mean=median(diff_dy),sd=mad(diff_dy)),lty=3,lwd=4)
legend( "topright", c("KDE","normal(mean,sd)","normal(median,mad)"),lwd=c(2,2,4),lty=c(1,5,3))
if( save_plots ){ dev.off() }


# E 4 EPage 77
# 
mk_p_plots = function(y){
  n = length(y)
  q.grid = (1:n)/(n+1)
  p.grid = c(0.25,0.1,0.05,0.025,0.01,0.0025)
  par(mfrow=c(3,2))
  for(j in 1:6){
    p_value = p.grid[j]
    qqplot(y, qnorm(q.grid), main=sprintf("p=%.4f", p_value))
    abline(lm( qnorm(c(p_value,1-p_value)) ~ quantile(y,c(p_value,1-p_value))))
  }
}

diffbp = diff(Garch$bp)
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter4/chap_4_ex_4_qq_diffbp.eps", onefile=FALSE, horizontal=FALSE) }
mk_p_plots(diffbp)
if( save_plots ){ dev.off() }

# Do the same thing as above but with n simulated N(0,1) random variables 
# 
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter4/chap_4_ex_4_qq_sim_normal.eps", onefile=FALSE, horizontal=FALSE) }
mk_p_plots(rnorm(length(diffbp)))
if( save_plots ){ dev.off() }

# Do the same but for diff(log(bp)):
#
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter4/chap_4_ex_4_qq_difflbp.eps", onefile=FALSE, horizontal=FALSE) }
difflbp = diff(log(Garch$bp))
mk_p_plots(difflbp)
if( save_plots ){ dev.off() }


# E 6:
#
print( c(qnorm(0.025), qnorm(0.975)) )
print( qnorm(0.975, mean=1, sd=sqrt(2)) )