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


# 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]))
}
par(mfrow=c(1,1))
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)) ))
}
par(mfrow=c(1, 1))
if( save_plots ){ dev.off() }


# P 6: EPage 76
#
# Install some needed packages if they are not already installed:
# 
if( "fGarch" %in% rownames(installed.packages()) == FALSE ){
    install.packages("fGarch")
}
library("fGarch")

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter4/chap_4_prob_6_DAX_return_de.eps", onefile=FALSE, horizontal=FALSE) }

x = seq(-0.1,+0.1,by=0.001)
col_indx = 1 # study the DAX index returns
df = 4
mad_t = mad(logR[, col_indx], constant=sqrt( df / (df - 2) ) / qt(0.75, df) )

par(mfrow=c(1, 3))

plot( density(logR[,col_indx]), lwd=2, ylim=c(0,60), main='both tails' )
lines(x, dstd(x, mean=mean(logR[,col_indx]), sd=mad_t, nu=df),lty=5, lwd=2, col='red')
lines(x, dnorm(x, mean=mean(logR[,col_indx]), sd=sd(logR[,col_indx])), lty=3, lwd=4, col='blue')
legend( 'topleft', c('KDE', paste('t: df= ', df),'normal'), lwd=c(2,2,4), lty=c(1,5,3), col=c('black', 'red', 'blue'))

plot( density(logR[,col_indx]), lwd=2, ylim=c(0,20), xlim=c(-0.05,-0.0), xlab='', main='left tail' )
lines(x, dstd(x, mean=mean(logR[,col_indx]), sd=mad_t, nu=df),lty=5, lwd=2, col='red')
lines(x, dnorm(x, mean=mean(logR[,col_indx]), sd=sd(logR[,col_indx])), lty=3, lwd=4, col='blue')

plot( density(logR[,col_indx]), lwd=2, ylim=c(0,20), xlim=c(0.0, 0.05), xlab='', main='right tail' )
lines(x, dstd(x, mean=mean(logR[,col_indx]), sd=mad_t, nu=df),lty=5, lwd=2, col='red')
lines(x, dnorm(x, mean=mean(logR[,col_indx]), sd=sd(logR[,col_indx])), lty=3, lwd=4, col='blue')

par(mfrow=c(1, 1))

if( save_plots ){ dev.off() }


# P 8: Try to fit a t-distribution to the CAC index returns
#
# First use QQ plots to estimate a good value for the degrees of freedom:
#
n = dim(logR)[1]
q.grid = (1:n)/(n+1)
df = c(1,4,6,10,20,30)
col_index=3 # the CAC index returns
par(mfrow=c(3,2))
for(j in 1:6){
  qqplot(logR[,col_index], qt(q.grid,df=df[j]), main=paste(index.names[col_index], ", df=", df[j]))
  abline(lm( qt(c(0.25,0.75),df=df[j]) ~ quantile(logR[,col_index],c(0.25,0.75)) ))
}
par(mfrow=c(1,1))

# With that value compare the parametric fit with that dof to the true density:
#
x = seq(-0.1,+0.1,by=0.001)
col_indx = 3 # study the CAC index returns
df = 6
mad_t = mad(logR[, 1], constant=sqrt( df / (df - 2) ) / qt(0.75, df) )
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter4/chap_4_prob_8_CAC_return_de.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(1, 3))

plot( density(logR[,col_indx]), lwd=2, ylim=c(0,60), main='both tails' )
lines(x, dstd(x, mean=mean(logR[,col_indx]), sd=mad_t, nu=df),lty=5, lwd=2, col='red')
lines(x, dnorm(x, mean=mean(logR[,col_indx]), sd=sd(logR[,col_indx])), lty=3, lwd=4, col='blue')
legend( 'topleft', c('KDE', paste('t: df= ', df),'normal'), lwd=c(2,2,4), lty=c(1,5,3), col=c('black', 'red', 'blue'))

plot( density(logR[,col_indx]), lwd=2, ylim=c(0,20), xlim=c(-0.05,-0.0), xlab='', main='left tail' )
lines(x, dstd(x, mean=mean(logR[,col_indx]), sd=mad_t, nu=df),lty=5, lwd=2, col='red')
lines(x, dnorm(x, mean=mean(logR[,col_indx]), sd=sd(logR[,col_indx])), lty=3, lwd=4, col='blue')
plot( density(logR[,col_indx]), lwd=2, ylim=c(0,20), xlim=c(0.0, 0.05), xlab='', main='right tail' )
lines(x, dstd(x, mean=mean(logR[,col_indx]), sd=mad_t, nu=df),lty=5, lwd=2, col='red')
lines(x, dnorm(x, mean=mean(logR[,col_indx]), sd=sd(logR[,col_indx])), lty=3, lwd=4, col='blue')
par(mfrow=c(1, 1))
if( save_plots ){ dev.off() }


# P 9
#
data = read.csv('../../BookCode/Data/MCD_PriceDaily.csv', header=TRUE)
data$Date = as.Date(data$Date, '%m/%d/%Y')
aderice = data[, 7]
plot(data$Date, aderice, type='l', lwd=2)
grid()


# P 10
#
LogRet = diff(log(aderice))
plot(data$Date[-1], LogRet, type='l', lwd=2, xlab='Date', ylab='logreturn')
grid()


# P 11
#
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter4/chap_4_prob_11_plots.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(1, 2))
hist(LogRet, 80, freq=FALSE)
qqnorm(LogRet, datax=T)
qqline(LogRet, datax=T)
par(mfrow=c(1, 1))
if( save_plots ){ dev.off() }

print(shapiro.test(LogRet))