##############################
# Markowitz Portfolio Theory #
##############################
# Part 1
# Two Asset Portfolio 
p1 = c(100, 110, 108, 120, 117)
p1lag = c(NA, 100, 110, 108, 120)
r1 = (p1 - p1lag)/p1lag
mean(r1, na.rm=T)
sd(r1, na.rm=T)

mu1 = mean(r1, na.rm=T); sig1 = sd(r1, na.rm=T)
mu2 = 0.03; sig2 = 0.05
rho12 = 0
w1 = seq(0,1,0.01)
w2 = 1 - w1
mup = w1*mu1 + w2*mu2 
sigp = sqrt(w1^2*sig1^2+w2^2*sig2^2+2*w1*w2*rho12*sig1*sig2)
plot(sigp,mup,type="l",xlab="sigmap",ylab="mup",main="mu-sigma diagram")

f_sigp = function(rho) {
return(sqrt(w1^2*sig1^2+w2^2*sig2^2+2*w1*w2*rho*sig1*sig2))
}

plot(f_sigp(-1),mup,type="l",xlab="sigmap",ylab="mup",main="mu-sigma diagram",col="red")
lines(f_sigp(0),mup,col="blue")
lines(f_sigp(1),mup,col="black")

# Part 2
# Three assets---constrainted portfolio mean
# Rational investor should let targmiu be greater than mean of min-var portfolio
# Quadratic Programing
library(quadprog)
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)
targmiu = 7

Dmat = 2*sigma
dvec = c(0,0,0)
Amat = matrix(c(1,1,1,miu1,miu2,miu3),nrow=3,ncol=2)
bvec = c(1, targmiu)
ou = solve.QP(Dmat,dvec,Amat,bvec,meq=2)
cat("Optimal weight is ", ou$solution, "\n")
cat("Min-standard deviation is ", sqrt(ou$value), "\n")

library(tseries)
op = portfolio.optim(t(miuv), pm=7,covmat=sigma)
op$pw

* allow for shorting
portfolio.optim(t(miuv), pm=9,covmat=sigma,shorts=T)


# Part 3
# Set of all three-asset portfolios (three-asset bullet) 
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)

miup = NULL
sigp = NULL
wmatrix = NULL
w1 = 0
while (w1<=1) {
	w2 = 0
	while (w2<=1-w1) {
		w3 = 1 - w1 - w2
		w = c(w1,w2,w3)
		miup = c(miup,w%*%miuv)
		sigp = c(sigp,sqrt(t(w)%*%sigma%*%w))
                wmatrix = cbind(wmatrix,w)
		w2 = w2 +0.01
	}
	w1 = w1 + 0.01
}

plot(sigp,miup, main="Set of Portfolios with Three Assets")
text(4+0.15, 5-0.02,"Asset 1")
text(5-0.1, 8+0.02,"Asset 2")
text(4.5+0.2, 7+0.02,"Asset 3")

# GMVP obtained by grid search
min(sigp)
miup[which.min(sigp)]
wmatrix[,which.min(sigp)]



# Efficiency Frontier 
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)

targetmiu = seq(5,8,0.1)
sigp = rep(0, length(targetmiu))
for (i in 1:length(targetmiu)) {
op = portfolio.optim(t(miuv), pm=targetmiu[i],covmat=sigma)
sigp[i] = sqrt(t(op$pw)%*%sigma%*%op$pw)
}

# GMVP obtained by quadratic programming or portfolio.optim
targetmiu[which.min(sigp)]
min(sigp)
portfolio.optim(t(miuv), pm=6.29,covmat=sigma)

# Plot efficiency frontier
plot(sigp, targetmiu, main="Red Efficiency Frontier", col = ifelse(targetmiu>6.3, "red", "black"))


# Part 4
# complete picture 
miu1 = 0.05; miu2=0.10; miu3=0.15
miuv = c(miu1, miu2, miu3)
sigma = matrix(c(0.25,0.15,0.17,0.15,0.21, 0.09, 0.17,0.09,0.28),nrow=3,ncol=3,byrow=T)

combination.line = function(muv, Sig) {
w1 = seq(-2,2,0.02)	
w2 = 1-w1
mupf = w1*muv[1] + w2*muv[2] 
sigpf = sqrt(w1^2*Sig[1,1]+w2^2*Sig[2,2]+2*w1*w2*Sig[1,2])	
return(cbind(mupf,sigpf))	
}

cl12 = combination.line(miuv[-3],sigma[-3,-3])
cl13 = combination.line(miuv[-2],sigma[-2,-2])
cl23 = combination.line(miuv[-1],sigma[-1,-1])

sigv = sqrt(c(sigma[1,1],sigma[2,2],sigma[3,3]))

minimumvarianceset = function(muv, Sig, targetrange) {
targetmiu = seq(targetrange[1],targetrange[2],0.02)
sigp = rep(0, length(targetmiu))
for (i in 1:length(targetmiu)) {
op = portfolio.optim(t(muv), pm=targetmiu[i],covmat=Sig, shorts=TRUE)
sigp[i] = sqrt(t(op$pw)%*%Sig%*%op$pw)
}
return(cbind(targetmiu,sigp))
}

mvs = minimumvarianceset(miuv,sigma,c(min(cl13[,1]),max(cl13[,1])))
plot(mvs[,2],mvs[,1], xlab="sigmap",ylab="mup",type="l",lwd=2, main="Minimum Variance Set and Combination Line")
lines(cl12[,2],cl12[,1],col="red",lty=2,lwd=1.5)
lines(cl13[,2],cl13[,1], col="blue",lty=2,lwd=1.5)
lines(cl23[,2],cl23[,1], col="green",lty=2,lwd=1.5)
points(sigv, miuv,col="red",pch=4)
text(sigv[1]-0.04, miuv[1]-0.02,"Asset 1")
text(sigv[3]-0.03, miuv[3]+0.01,"Asset 3")


# straight combination line with risk free asset 
miuz = 2; miut = 6; sigz = 0; sigt = 5; sigzt = 0
w = seq(0,2,0.01)
miuzt = w*miut + (1-w)*miuz
sigzt = sqrt(w^2*sigt^2+(1-w)^2*sigz^2+2*w*(1-w)*sigzt)
plot(sigzt,miuzt,main="combination line with risk free asset")