#
# 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_plot = F

set.seed(0)

# Ex. 1:
#
library(Ecdat)
data(CRSPday)
dimnames(CRSPday)[[2]]

r = CRSPday[,5]
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter5/ex_1_returns_plot.eps", onefile=FALSE, horizontal=FALSE) }
plot(r)
if( save_plots ){ dev.off() }

mode(r)
class(r)

r2 = as.numeric(r)
class(r2)
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter5/ex_1_numeric_returns_plot.eps", onefile=FALSE, horizontal=FALSE) }
plot(r2)
if( save_plots ){ dev.off() }

cov(CRSPday[,4:6])
cor(CRSPday[,4:6])
apply(CRSPday[,4:6],2,mean)


# Ex. 5:
#
s2 = 10
ps = seq( 0, 1, length.out=100 )
numer = 3*(1-ps)*s2^2 + 3*ps
denom = ((1-ps)*s2 + ps)^2 
kur = numer / denom
plot( ps, kur, type="l" ) 
ps_max = s2*(1-s2) / ( 1 - s2^2 )
abline( v=ps_max, col="red" )

# Lets find values of p and sigma so that Kur(X) >= 10000
s2 = 20000
p = s2*(1-s2)/(1-s2^2) 
kur = (3*(1-p)*s2^2 + 3*p)/(((1-p)*s2 + p)^2)

# Ex. 6:
#
library("fGarch")
dat = read.csv("../../BookCode/Data/GasFlowData.csv",header=T)
fit_sstd = sstdFit( dat$Flow1 )

if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter5/ex_6_density_plots.eps", onefile=FALSE, horizontal=FALSE) }
plot( density( dat$Flow1 ), type="l", xlab="x", ylab="density" )
x_values = seq( min(dat$Flow1), max(dat$Flow1), length.out=100 )
y_values = dsstd( x_values, mean=fit_sstd$estimate[1], sd=fit_sstd$estimate[2], nu=fit_sstd$estimate[3], xi=fit_sstd$estimate[4] )
lines(x_values,y_values,lwd=2,lty=5)
legend("topleft",c("KDE","sstd"),lwd=2,lty=c(1,5))
if( save_plots ){ dev.off() }


# Ex. 8:
#
28.0834 + 1.8098 * qnorm(1-0.05/2) * c(-1,+1)  
0.6884 + 0.1638 * qnorm(1-0.05/2) * c(-1,+1)


# Ex. 9:
#
library(evir)
library(fGarch)
data(bmw)
start_bmw = c( mean(bmw), sd(bmw), 4)
loglik_bmw = function(theta){
    -sum(log(dstd(bmw,mean=theta[1],sd=theta[2],nu=theta[3])))
}

mle_bmw = optim(start_bmw, loglik_bmw, hessian=T) # find the ML parameters 
FishInfo_bmw = solve( mle_bmw$hessian )
sqrt( diag( FishInfo_bmw ) )


# Ex. 10:
#
data(siemens)
n = length(siemens)
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter5/ex_10_qqplots.eps", onefile=FALSE, horizontal=FALSE) }
par(mfrow=c(3,2))
qqplot(siemens, qt(((1:n)-0.5)/n,2), ylab="t(2) quantiles", xlab="data quantiles")
qqplot(siemens, qt(((1:n)-0.5)/n,3), ylab="t(3) quantiles", xlab="data quantiles")
qqplot(siemens, qt(((1:n)-0.5)/n,4), ylab="t(4) quantiles", xlab="data quantiles")
qqplot(siemens, qt(((1:n)-0.5)/n,5), ylab="t(5) quantiles", xlab="data quantiles")
qqplot(siemens, qt(((1:n)-0.5)/n,8), ylab="t(8) quantiles", xlab="data quantiles")
qqplot(siemens, qt(((1:n)-0.5)/n,12), ylab="t(12) quantiles", xlab="data quantiles")
if( save_plots ){ dev.off() }

library("MASS")
fitdistr(siemens,"t") # fits a t-model.  This function is in the MASS package