##################################################################
### Example 7.8: Thinning vs. historical simulation
##################################################################
#
# Version: 2013-11-07
#

data <- read.csv("Data/Dow_Jones.csv") # Read data into R
numberofdays <- 801 # Selecting the size of the sample
data <- data[1:numberofdays,]
# Converting calendar dates into R format
dates <- rev(as.character(data[,1]))
dates <- as.Date(dates, "%m/%d/%y")

# plotting the dow jones index
index <- rev(data[,7]) # Using the adjusted Dow Jones index
plot(dates, index, type="l", main="Dow-Jones Index")

# plotting log returns
logret <- diff(log(index))
plot(dates[2:length(index)], logret, type="h", xlab="dates", main="Dow-Jones logreturns")

# histogram of logreturns
h <- hist(logret,breaks = 50, probability = TRUE) # Constructing histogram
m <- mean(logret) # Estimating the mean
sigma <- sd(logret) # Estimating the standard deviations
lines(h$breaks,dnorm((h$breaks-m)/sigma)/sigma) # Plotting the fitted normal dist

n <- numberofdays-1
u <- (1:n)/(n+1)
qqplot(qnorm(u), logret, xlab = "Normal quantiles", ylab = "sample quantiles")

# Constructing 20 day log-returns by throwing away data.
newindex <- c()
for(k in 0:40)
{
  newindex <- c(newindex, index[20*k+1])
}

newlogret <- diff(log(newindex))

n <- length(newlogret)
networth <- 100*(exp(newlogret)-1)

snw <- sort(networth)
stair <- c(1:n)/n

plot(c(snw[1]-1,snw[1:n]),c(0,stair[1:n]), xlab = " ", ylab = " ", ylim = c(0,1), xlim = c(snw[1]-1,0), type = "s")  


xlim= c(snw[1]-1,0)

# Historical simulation
N <- 5000
sample <- c()
for(j in 1:N){
  tmpsample <- c()
  for(k in 1:20){
    tmpsample <- c(tmpsample, logret[floor(runif(1,1,length(logret)+1))])
  }
  sample <- c(sample, sum(tmpsample))
}

hlosses <- 100*(exp(sample)-1)
hsnw <- sort(hlosses)
hstair <- c(1:N)/N
lines(hsnw,hstair, type = "s")



# Plot of VaR

plot(c(0,stair),c(-snw[1]+1,-snw),  xlab = " ", ylab = " ", xlim = c(0,0.2), ylim= c(0,26), type = "s")
lines(hstair,-hsnw,  xlab = " ", ylab = " ", xlim = c(0,0.2), type = "s")