#
# 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-2: EPage 74
#
data(EuStockMarkets)
mode(EuStockMarkets)
class(EuStockMarkets)

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter4/chap_4_prob_1_eu_stocks.eps", onefile=FALSE, horizontal=FALSE) }
plot(EuStockMarkets) 
if( save_plots ){ dev.off() }

logR = diff(log(EuStockMarkets))

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter4/chap_4_prob_1_eu_stocks_diff.eps", onefile=FALSE, horizontal=FALSE) }
plot(logR)
if( save_plots ){ dev.off() }

# P 3: EPage 74
#
plot(as.data.frame(logR))

index.names = dimnames(logR)[[2]]

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter4/chap_4_prob_1_eu_stocks_qq_plots.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(2,2))
for(i in 1:4){
  qqnorm(logR[,i],datax=T,main=index.names[i])
  qqline(logR[,i],datax=T)
  print(shapiro.test(logR[,i]))
}
if( save_plots ){ dev.off() }

# P 5: EPage 74
#
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter4/chap_4_prob_5_DAX_qq_plots.eps", onefile=FALSE, horizontal=FALSE) }
n = dim(logR)[1]
q.grid = (1:n)/(n+1)
df = c(1,4,6,10,20,30)
i=1 # Just study the DAX index 
par(mfrow=c(3,2))
for(j in 1:6){
  qqplot(logR[,i], qt(q.grid,df=df[j]), main=paste(index.names[i], ", df=", df[j]))
  abline(lm( qt(c(0.25,0.75),df=df[j]) ~ quantile(logR[,i],c(0.25,0.75)) ))
}
if( save_plots ){ dev.off() }

# P 6: EPage 76
#
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter4/chap_4_prob_6_DAX_return_de.eps", onefile=FALSE, horizontal=FALSE) }
library("fGarch")
x = seq(-0.1,+0.1,by=0.001)
par(mfrow=c(1,1))
#plot( density(logR[,1]), lwd=2, ylim=c(0,60) )
plot( density(logR[,1]), lwd=2, ylim=c(0,20), xlim=c(-0.05,-0.0) )
lines(x,dstd(x,mean=median(logR[,1]),sd=mad(logR[,1]),nu=5),lty=5,lwd=2)
lines(x,dnorm(x,mean=median(logR[,1]),sd=mad(logR[,1])),lty=3,lwd=4)
legend( "topleft", c("KDE","t: df=5","normal"),lwd=c(2,2,4),lty=c(1,5,3))
if( save_plots ){ dev.off() }