sp500<- read.table("sp500 tickers 10-17.csv",sep="\t",header=T,quote="",colClasses="character",strip.white=T)
colnames(sp500) <- c("name","ticker")
#replace / by -
sp500$ticker <- sub('/', '-', sp500$ticker)

#first don't take care about survivor bias
tickerlist <- list()
for(i in 1:length(sp500$ticker)){
  ticker <- sp500[i,"ticker"]
  tickerlist[[ticker]] <- getYahooData(ticker,"20000101")
  cat("ticker",i,ticker,"\n")
}

#compute matrix of returns
ticker <- sp500[1,"ticker"]
returns <- diff(log(tickerlist[[ticker]]$Close))
for(i in 2:length(sp500[,"ticker"])){
  ticker <- sp500[i,"ticker"]
  returns <- merge(returns,diff(log(tickerlist[[ticker]]$Close)))
}
colnames(returns) <- sp500[,"ticker"]
returns <- returns[-1,]
universe <- which(!is.na(returns[1,]))
returns.subset <- returnds[,universe]

#constant weight rebalancing
w <- rep(1/dim(returns.subset)[2],dim(returns.subset)[2])
N <- 1000000
W <- N
for(day in 1:length(returns[,1])){
  r <- as.vector(exp(returns.subset[day,])-1)
  r[is.na(r)] <- 0 #NA means we don't have quotes that day, but next return
                   #will catch up
  dN <- N*crossprod(w,r)
  N <- N+dN
  W <- c(W,N)
}
#no rebalancing
N <- 1000000
W <- N
R <- rep(1,dim(returns.subset)[2])
for(day in 1:length(returns[,1])){
  r <- as.vector(exp(returns.subset[day,]))
  r[is.na(r)] <- 1
  R <- R*r
  W <- c(W,N*mean(R))
}

indexchange <- read.table("IndexChange_500.csv",sep="\t",header=T,quote="",
                          colClasses="character")
#cleanup tickers
indexchange$Ticker.add <- sub('/', '-',indexchange$Ticker.add)
indexchange$Ticker.add <- sub('.wi', '',indexchange$Ticker.add,fixed=T)
indexchange$Ticker.add <- sub('wi', '',indexchange$Ticker.add)
indexchange$Ticker.add <- sub('.', '-',indexchange$Ticker.add,fixed=T)
indexchange$Ticker.del <- sub('/', '-',indexchange$Ticker.del)
indexchange$Ticker.del <- sub('.wi', '',indexchange$Ticker.del,fixed=T)
indexchange$Ticker.del <- sub('wi', '',indexchange$Ticker.del)
indexchange$Ticker.del <- sub('.', '-',indexchange$Ticker.del,fixed=T)

#rollback sp500 changes
sp500.historic <- sp500
for( i in 1:100){
  #Ticker.add empty means index is <500 for a while
  if(indexchange$Ticker.add[i]!=""){
    if(!is.na(match(indexchange$Ticker.add[i],sp500.historic$ticker))){
      sp500.historic <-
        sp500.historic[-match(indexchange$Ticker.add[i],sp500.historic$ticker),]
    }else{
      cat(paste("ticker",indexchange$Ticker.add[i],"not found\n"))
      next()
    }
  }
   #add ticker that was removed
  if(indexchange$Ticker.del[i]!=""){
    sp500.historic <- rbind(sp500.historic,
                            c(indexchange$Company.Deletions[i],
                              indexchange$Ticker.del[i]))
  }
}

n.stocks <- ncol(P)
N <- nrow(P)
current.wealth <- 1
w <- rep(1/n.stocks,n.stocks)
dP <- apply(P,2,diff)
L <- rep(0,n.stocks)
V <- 0
W <- w
pnl <- 0
for(i in 1:(N-1)){
  r <- dP[i,]/as.numeric(P[i,])
  #r <- c(r,-r)
  pnl <- cbind(pnl,current.wealth*crossprod(w,r))
  current.wealth <- 1+sum(pnl)
  L <- rbind(L,-r)
  V <- V+crossprod(w,r^2)
  T <- 2/3*sqrt(log(N)/V)
  #w <- exp(-T*(apply(1+L,2,prod)-1))
  w <- exp(-T*colSums(L))
  w <- w/sum(w)
  W <- rbind(W,w)
  if(i%%10==0){
    cat(current.wealth,sep="\n")
  }
}

price2return <- function(x){
  diff(x)/x[-length(x)]
}
#number of shares implementations
tc <- 0.005+0.02 
days <- nrow(data.bus)
init.capital <- 1000000
tickers <- memb(sp500.tickers,add,as.Date(time(data.bus)[1]))
tickers.index <- which(sp500.tickers.extended%in%tickers)
n.stocks <- length(tickers)
w <- rep(0,length(sp500.tickers.extended))
w[tickers.index] <- rep(1/n.stocks,n.stocks)
N <- round((capital*w)/as.numeric(data.bus[1,tickers.index]))
dP <- apply(data.bus,2,diff)
L <- rep(0,n.stocks)
V <- 0
W <- w
pnl <- 0
tcvec <- sum(N)*tc
for(i in 1:days){
  tickers <- memb(sp500.tickers,add,as.Date(time(data.bus)[1]))
  tickers.index <- which(sp500.tickers.extended%in%tickers)
  r <- dP[i,]/as.numeric(data.bus[i,])
  pnl <- cbind(pnl,crossprod(N,dP[i,]))
  capital <- init.capital+sum(pnl)-sum(tcvec)
  L <- rbind(L,-r)
  V <- V+crossprod(w,r^2)
  T <- 2/3*sqrt(log(days)/V)
  #w <- exp(-T*(apply(1+L,2,prod)-1))
  w <- exp(-T*colSums(L))
  w <- w/sum(w)
  newN <- round((capital*w)/as.numeric(data.bus[(i+1),tickers.index]))
  tcvec <- c(tcvec,sum(abs(newN-N)*tc))
  N <- newN
  if(i%%10==0){
    cat(capital,sep="\n")
  }
}

fixed.rebal <- function(P,delta){
  init.capital <- 1
  dP <- apply(P,2,diff)
  capital <- init.capital
  pnl <- c()
  for(i in 1:(nrow(P)-1)){
    r <- dP[i,]/P[i,]
    pnl <- c(pnl,capital*crossprod(delta,r))
    capital <- init.capital+sum(pnl)
  }
  return( pnl )
}