## R script for "Weight compared to risk" post 2011/04/18

# location of list of stock symbols
sp500.symbol.url <- "http://blog.quanttrader.org/wp-content/uploads/sp500.csv"

# get function 'pp.TTR.multsymbol' to read in data (and other functions to ignore)
source('http://www.portfolioprobe.com/R/pprobe_functions01.R')

# get character vector of stock symbols
sp500.symbols <- scan(url(sp500.symbol.url), what="")

# read data in
require(TTR)
sp500.close <- pp.TTR.multsymbol(sp500.symbols, 20070101, 20110415)

# delete stocks with missing data
# a real analysis would want to be more careful to include all data
sp500.closeok <- sp500.close[, colSums(is.na(sp500.close)) == 0]

# create returns from closing prices
sp500.ret <- diff(log(sp500.closeok))

# find the row of the last day of the desired quarter
match("2008-09-30", as.character(index(sp500.ret)))

# compute the variancce matrix
require(BurStFin) # available via burns-stat.com
sp500.var08Q3 <- var.shrink.eqcor(sp500.ret[seq(to=440, length=250),])

# get vector of prices at the end of the quarter
sp500.price08Q3 <- drop(as.matrix(sp500.closeok[440,]))

# generate random portfolios
require(PortfolioProbe)
rp.08Q3.w05 <- random.portfolio(1000, prices=sp500.price08Q3, gross=1e6, long.only=TRUE, max.weight=.05,   port.size=c(40,50))
rp.08Q3.rf05 <- random.portfolio(1000, prices=sp500.price08Q3, gross=1e6, long.only=TRUE, risk.fraction=.05,   port.size=c(40,50), variance=sp500.var08Q3)

# get weight vectors of random portfolios
rp.08Q3.w05.weights <- randport.eval(rp.08Q3.w05, FUN=function(x) valuation(x)$weight)
rp.08Q3.rf05.weights <- randport.eval(rp.08Q3.rf05, FUN=function(x) valuation(x)$weight)

# get risk fractions of random portfolios
rp.08Q3.rf05.rfrac <- randport.eval(rp.08Q3.rf05, keep="risk.fraction")
rp.08Q3.w05.rfrac <- randport.eval(rp.08Q3.w05, additional.args=list(risk.fraction=.05, variance=sp500.var08Q3),   keep="risk.fraction")

# get risk fractions into a more convenient form
rp.08Q3.rf05.rfracred <- unlist(rp.08Q3.rf05.rfrac)
rp.08Q3.rf05.rfracred <- rp.08Q3.rf05.rfracred[rp.08Q3.rf05.rfracred != 0]

rp.08Q3.w05.rfracred <- unlist(rp.08Q3.w05.rfrac)
rp.08Q3.w05.rfracred <- rp.08Q3.w05.rfracred[rp.08Q3.w05.rfracred != 0]

# sort of reproduce Figure 1
w05.1.w <-  rp.08Q3.w05.weights[[1]]
w05.1.rf <-  rp.08Q3.w05.rfrac[[1]][names(w05.1.w),]
w05.1.rf <-  rp.08Q3.w05.rfrac[[1]][names(w05.1.w)]
w05.1.rf <-  rp.08Q3.w05.rfrac[[1]][[1]][names(w05.1.w),]
plot(w05.1.w, w05.1.rf)

# Figure 2
plot(density(rp.08Q3.w05.rfracred * 100, from=0), xlim=c(0,12), lwd=3, 
   col='blue', ylab="", yaxt="n", xlab="Risk fraction (percent)", main="")
mtext("Density", side=2, line=1)
lines(density(rp.08Q3.rf05.rfracred * 100, from=0, to=5), lwd=3, col='gold')

# Figure 3
plot(density(unlist(lapply(rp.08Q3.w05.rfrac, function(x) max(x[[1]])))  * 100), lwd=3, col='blue', 
  ylab="", yaxt="n", xlab="Maximum risk fraction (percent)", main="")
mtext("Density", side=2, line=1)

# analysis for 2011Q1

match("2011-03-31", as.character(index(sp500.closeok)))
sp500.price11Q1 <- sp500.closeok[1070,]
sp500.var11Q1 <- var.shrink.eqcor(sp500.ret[seq(to=1070, length=250),])
sp500.var11Q1 <- var.shrink.eqcor(sp500.ret[seq(to=1070, length=250),])
sp500.price11Q1 <- drop(as.matrix(sp500.closeok[1070,]))
rp.11Q1.w05 <- random.portfolio(1000, prices=sp500.price11Q1, gross=1e6, long.only=TRUE, max.weight=.05,   port.size=c(40,50))
rp.11Q1.rf05 <- random.portfolio(1000, prices=sp500.price11Q1, gross=1e6, long.only=TRUE, risk.fraction=.05,   port.size=c(40,50), variance=sp500.var11Q1)
rp.11Q1.rf05.weights <- randport.eval(rp.11Q1.rf05, FUN=function(x) valuation(x)$weight)
rp.11Q1.w05.weights <- randport.eval(rp.11Q1.w05, FUN=function(x) valuation(x)$weight)
rp.08Q3.rf05.rfrac <- randport.eval(rp.08Q3.rf05, keep="risk.fraction")
rp.08Q3.rf05.rfrac <- randport.eval(rp.08Q3.rf05, keep="risk.fraction")
rp.11Q1.rf05.rfrac <- randport.eval(rp.11Q1.rf05, keep="risk.fraction")
rp.11Q1.w05.rfrac <- randport.eval(rp.11Q1.w05, additional.args=list(risk.fraction=.05,   variance=sp500.var11Q1),keep='risk.fraction')
rp.11Q1.w05.rfracred <- unlist(rp.11Q1.w05.rfrac)
rp.11Q1.w05.rfracred <- rp.11Q1.w05.rfracred[rp.11Q1.w05.rfracred != 0]
rp.11Q1.rf05.rfracred <- unlist(rp.11Q1.rf05.rfrac)
rp.11Q1.rf05.rfracred <- rp.11Q1.rf05.rfracred[rp.11Q1.rf05.rfracred != 0]