gpdfitml<-function(datavec,initguess)
{
	n <- length(datavec)
	optfun <- function(params)
	{
		n*log(params[2])+(1/params[1]+1)*sum(log(1+params[1]/params[2]*datavec))
	}
	optim(initguess,optfun)
}

qpareto <- function(p,gamma,beta)
{
	(beta/gamma) * ( (1-p)^(-gamma)-1 )
}

# Program

# Read data
data <- read.csv("Data/Nasdaq.csv",header=TRUE)
#data <- read.csv("Data/Dow_Jones.csv")
vec <- data[,7]
price <- rev(vec[1:2100])
return <- -diff(log(price))
N <- length(return)
z <- (-10:10)

#Fit data to Normal distribution
library(MASS)
nparams<-fitdistr(return,densfun="normal")$estimate
u <- (1:N)/(N+1)
qqplot(nparams[1]+nparams[2]*qnorm(u), return, xlab = "Normal quantiles", ylab = "sample quantiles")

#Choose the threshold so that we only consider [n] most extreme datapoints
n <- 50
return_new <- sort(return,dec=T)[1:n] - sort(return,dec=T)[n+1]
nparams <- fitdistr(return_new,densfun="normal")$estimate
pparams <- gpdfitml(return_new,c(0.1,0.1))$par
u <- 1-n/N + (1:n)/(N+1)
qqplot(nparams[1]+nparams[2]*qnorm(u), return_new, xlab = "Normal quantiles", ylab = "sample quantiles")
qqplot(qpareto(u,pparams[1],pparams[2]), return_new, xlab = "Pareto quantiles", ylab = "sample quantiles")