##################################################################
### Example 9.6: Estimation and fit of an Elliptical model
##################################################################

data <- read.csv("Data/Dow_Jones.csv") # Read data into R
numberofdays <- 501 # Selecting the size of the sample
data <- data[1:(1+numberofdays),]
# Converting calendar dates into R format
dates <- rev(as.character(data[,1]))
dates.dj <- as.Date(dates, "%m/%d/%y")
# plotting the dow jones index
index.dj <- rev(data[,7]) # Using the adjusted Dow Jones index
logret.dj <- diff(log(index.dj))

data <- read.csv("Data/Nasdaq.csv", header = TRUE)
data <- data[3:(3+numberofdays),]
dates.nq <- as.Date(data[,1], "m/%d/%y")
index.nq <- rev(data[,7])
logret.nq <- diff(log(index.nq)) # Calculating log returns
logret <- matrix(c(logret.dj-mean(logret.dj),logret.nq-mean(logret.nq)),length(logret.dj),2) 


S <- cov(logret)
n <- numberofdays-1
u <- (1:n)/(n+1)

###
# Selecting the random portfolios
###

N <- 100
z <- rnorm(2*N)
dim(z) <- c(N,2)
w <- z/sqrt(rowSums(z^2))
plot(w[,1],w[,2]) #plot the portfolio weights


###
# QQplots for all portfolios
###

portfolios <- w%*%t(logret)
pquants <- matrix(NA,N,n)
for(i in 1:N){
weight <- 1/sqrt(t(w[i,])%*%S%*%w[i,])
  for(j in 1:n){
  pquants[i,j] <- weight*empquant(portfolios[i,],u[j])
}
}

plot(qnorm(u), pquants[1,], ylim = c(-6,6), xlab="",ylab="")
for(i in 1:N){
  points(qnorm(u),pquants[i,])
}



gq <- colMeans(pquants)
lines(qnorm(u),gq,col="red")
x <- qnorm(u)
x3 <- qnorm(u)^3

y <- gq
mm <- glm(gq~1+x+x3)
inter <- mm$coefficients[[1]]
c1 <- mm$coefficients[[2]]
c3 <- mm$coefficients[[3]]
lines(x,inter+c1*x+c3*x3,col="yellow")

lstlsfit.A<-function(datavec)
{
n<-length(datavec)
sdatavec<-sort(datavec,dec=FALSE)
locscatfun<-function(params){
	sum((sqrt((params-2)/params)*qt((1:n)/(n+1),params)-sdatavec)^2)}
optimize(locscatfun, c(2,100))
}



lstlsfit.B<-function(datavec,initguess)
{
n<-length(datavec)
sdatavec<-sort(datavec,dec=F)
locscatfun<-function(params){sum((params[1]*qt((1:n)/(n+1),params[2])-sdatavec)^2)}
optim(initguess,locscatfun)
}


###
# Method A
###

res.A <- lstlsfit.A(y)
res.B <- lstlsfit.B(y,c(0.01,3))

df <- res.A$minimum
xx <- sqrt((df-2)/df)*qt(u,df)
lines(qnorm(u),xx,col="purple")
df
yy.A <- mean(logret.nq)+sqrt(S[2,2])*(xx)

# Comparing the Student's t and the polynomial Normal models

plot(qnorm(u),yy.A, type="l", xlab="", ylab="")
y <- (1.46+12.3*x-0.54*x^2+1.98*x^3)/1000
lines(qnorm(u),y,lty="dashed")

plot(y,yy,type="l")

###
# Method B
###

res.B <- lstlsfit.B(y,c(0.01,3))

df.B <- res.B$par[[2]]
scale.B <- res.B$par[[1]]

xx.B <- scale.B*qt(u,df.B)
lines(qnorm(u),xx.B,col="blue")
df.B
yy.B <- mean(logret.nq)+xx.B

# Comparing the Student's t and the polynomial Normal models

plot(qnorm(u),yy.B, type="l", xlab="", ylab="")
y <- (1.46+12.3*x-0.54*x^2+1.98*x^3)/1000
lines(qnorm(u),y,lty="dashed")

plot(y,yy.B,type="l")