##########################################################
# Finding Tangency Portfolio that maximizes Sharpe Ratio #
##########################################################
library(tseries)
miu1 = 5;s1 = 4; miu2=8; s2 = 5; miu3=7; s3 = 4.5; rho12 = -0.5; rho13 = -0.3; rho23 = 0.4
s12=rho12*s1*s2; s13=rho13*s1*s3; s23=rho23*s2*s3
miuv = c(miu1, miu2, miu3)
sigma = matrix(c(s1^2,s12,s13,s12,s2^2, s23, s13,s23,s3^2),nrow=3,ncol=3,byrow=T)

# analytical solution when risk-free return is 3
rf = 3
w = solve(sigma)%*%(miuv-rf)/sum(solve(sigma)%*%(miuv-rf))
w
cat("Expected Return of Tangency Portfolo is ", t(w)%*%miuv, "\n")
cat("Standard deviation of Tangency Portfolo return is ", sqrt(t(w)%*%sigma%*%w), "\n")

# grid search
targetmiu = seq(0,10,0.01)
sigp = rep(0, length(targetmiu))
for (i in 1:length(targetmiu)) {
op = portfolio.optim(t(miuv), pm=targetmiu[i],covmat=sigma, shorts=T)
sigp[i] = sqrt(t(op$pw)%*%sigma%*%op$pw)
}
sharpe.r = (targetmiu-rf)/sigp
cat("Expected Return of Tangency Portfolo is ", targetmiu[which.max(sharpe.r)], "\n")
op = portfolio.optim(t(miuv), pm=targetmiu[which.max(sharpe.r)],covmat=sigma, shorts=T)
op$pw

# drawing tangency portfolio 
a = rf
b = (t(w)%*%miuv-rf)/sqrt(t(w)%*%sigma%*%w)
plot(sigp, targetmiu, xlim=c(0,max(sigp)+1), main="Tangency Portfolio",type="l", lwd=1)
abline(a,b)