#
# Written by:
# -- 
# John L. Weatherwax                2009-04-21
# 
# email: wax@alum.mit.edu
# 
# Please send comments and especially bug reports to the
# above email address.
# 
#-----

save_plots = F

set.seed(0)

source('util_fns.R')

# Test the function "efficient_frontier" using the example on EPage 297:
#
library(Ecdat)
data(CRSPday)
R = 100*CRSPday[,4:6]
mean_vect = apply(R, 2, mean)
cov_mat = cov(R)
sd_vect = sqrt(diag(cov_mat))

# The target portfolio returns:
muP = seq(0.05,0.14, length.out=300)
mu_free = 1.3/253 

result = efficient_frontier(R, muP, mu_free)

plot( result$sdP, result$muP )


# P 1: EPage 306
#
dat = read.csv("../../BookCode/Data/Stock_FX_Bond.csv", header=T)
prices = cbind(dat$GM_AC, dat$F_AC, dat$CAT_AC, dat$UTX_AC, dat$MRK_AC, dat$IBM_AC)
stock_names = c("GM", "F", "CAT", "UTX", "MRK", "IBM")
n = dim(prices)[1]
returns = 100 * ( prices[2:n,]/prices[1:(n-1),] - 1 ) # turns returns into percentages 
#pairs(returns)
mean_vect = apply(returns, 2, mean)
cov_mat = cov(returns)
sd_vect = sqrt(diag(cov_mat))

mu_free = 3.0/365

# The target portfolio return:
muP = seq(min(mean_vect), max(mean_vect), length.out=500)

result = efficient_frontier(returns, muP, mu_free=mu_free, w_lower_limit=-0.1, w_upper_limit=0.5)
sdP = result$sdP

# Plot the efficient frontier:
#
if( save_plots ){ postscript("../../WriteUp/Graphics/Chapter11/chap_11_rlab_prob_1_ef.eps", onefile=FALSE, horizontal=FALSE) }
plot( sdP, muP, type="l", lty=3, xlab="sigma portfolio", ylab="mu porfolio", xlim=c(0,1.1*max(sd_vect)), ylim=c(0,1.1*max(mean_vect)) )
points(0, mu_free, cex=4, pch="*")
print("optimal porfolio weights")
ind = result$max_sharpe_ind
print( result$weights[ind,] )
lines( c(0,2), mu_free + c(0,2)*( muP[ind] - mu_free )/sdP[ind], lwd=4, lty=2 )
points(sdP[ind], muP[ind], cex=4, pch="*") # show tangency porfolio
ind = result$min_variance_ind 
points(sdP[ind], muP[ind], cex=2, pch="+") # show min variance porfolio 
inds = muP > muP[result$min_variance_ind]
lines(sdP[inds], muP[inds], type="l", lwd=2, xlim=c(0, 2.5)) # plot the efficient frontier
# Plot some individual equities:
for( ii in 1:dim(returns)[2] ){
  text(sd_vect[ii], mean_vect[ii], stock_names[ii], cex=1.5)
}
grid()
if( save_plots ){ dev.off() }


# P 2: EPage 306
#
ind = result$max_sharpe_ind

E_R_P = 0.0007 * 100
E_R_T = result$muP[ind] # the return of the tangency portfolio
c( E_R_P, E_R_T, mu_free ) 

omega = ( E_R_P - mu_free ) / ( E_R_T - mu_free )
print(sprintf("omega= %10.6f", omega))

print("Total porfolio weights")
ind = result$max_sharpe_ind
print( omega*result$weights[ind,] )

S = 100000
S * omega * result$weights[ind,]