##############################################################
### Section 4.2.3 Evaluating the methods on simulated data
##############################################################



###
# Functions
###

# Function to compute empirical portfolio weights for max-exp, min-var, and trade off problems

getempMVtotpos<-function(npos,nsample,meanvec,covmat,sigma2lim,mulim,c,R0)
{
	# npos = number of repeated samples
	# nsample = sample size of each sample
	# sigma2lim = variance constraint
	# mulim = expected return constraint
	# c = trade off parameter
	# R0 = return on risk free asset
	 
dimen<-dim(covmat)[1]
outvals<-matrix(0,4*dimen,npos)
R0vec<-matrix(R0,dimen,1)
tmpmean<-R0vec

for(i in (1:npos))
{
	tmpvals<-rmvn(nsample,meanvec,covmat) # observed sample 
	tmpmean<-matrix(rowMeans(tmpvals),dimen,1) # estimated mean
	tmpcovmat<-cov(t(tmpvals)) # estimated covariance

	# Compute portfolio weights for max-exp?
	loopvals<-solve(tmpcovmat)%*%(tmpmean-R0vec) # numerator
	tmpdenom<-sqrt(t(tmpmean-R0vec)%*%solve(tmpcovmat)%*%(tmpmean-R0vec)) # denominator
	outvals[1:2,i]<-loopvals*as.double(sqrt(sigma2lim)/tmpdenom) # max-exp portfolio weights

	# Compute portfolio weights for min-var
	tmpdenom<-t(tmpmean-R0vec)%*%solve(tmpcovmat)%*%(tmpmean-R0vec) #denominator 
	outvals[3:4,i]<-loopvals*as.double((mulim-R0)/tmpdenom)
	if(t(outvals[3:4,i])%*%tmpmean+(1-sum(outvals[3:4,i]))<1)
	{
		outvals[3,i]<-0
		outvals[4,i]<-0
	}
	# Compute portfolio weights for trade off 
	outvals[5:6,i]<-(1/as.double(c))*loopvals
	# Compute estimated expected excess return
	outvals[7:8,i]<-(tmpmean-R0vec)
	}
	outvals
}

# Function to compute the expected value and standard deviation of a portfolio

getempMVexpsds<-function(posvec,meanvec,covmat,R0)
{
outvals<-matrix(0,2,length(posvec[1,]))
for(i in (1:length(posvec[1,])))
{
outvals[1,i]<-sqrt(t(posvec[,i])%*%covmat%*%posvec[,i])
outvals[2,i]<-t(meanvec)%*%posvec[,i]+(1-sum(posvec[,i]))*R0
}	
outvals
}

# Generate samples from a multivariate normal distribution
rmvn<-function(n,meanvec,covmat)
{
library(MASS)
dimen<-dim(covmat)[1]
t(mvrnorm(n,meanvec,covmat))
}

###
# Set parameters
###

meanvec1 <- c(1.025,1.075)
# meanvec1 <- c(1.1,1.2)
sigma1 <- 0.3
sigma2 <- 0.5
sigmamat <- rbind(c(sigma1^2,0.5*sigma1*sigma2),c(0.5*sigma1*sigma2,sigma2^2))

meanvec1
sigmamat

R0 <- 1  # Return on risk free asset
sigma0 <- 0.3 # Variance constraint in max exp problem
ones <-matrix(1,2,1)
mulim <- R0+sigma0*sqrt(t(meanvec1-R0*ones)%*%solve(sigmamat)%*%(meanvec1-R0*ones)) # Mean constraint in min var problem
c <- sqrt(t(meanvec1-R0*ones)%*%solve(sigmamat)%*%(meanvec1-R0*ones))/sigma0 # Trade off parameter

###
# Theoretical portfolio weights
###

# Compute portfolio weights for max-exp
	R0vec <- R0*ones
	loopvals<-solve(sigmamat)%*%(meanvec1-R0vec) # numerator
	denom.max.exp <-sqrt(t(meanvec1-R0vec)%*%solve(sigmamat)%*%(meanvec1-R0vec)) # denominator
	max.exp.weights <-loopvals*as.double(sigma0/denom.max.exp) 
	max.exp.weights
	
	# Compute portfolio weights for min-var
	denom<-t(meanvec-R0vec)%*%solve(sigmamat)%*%(meanvec1-R0vec) 
	min.var.weights<-loopvals*as.double((mulim-R0)/denom)
	min.var.weights
	# Compute portfolio weights for trade off 
	trade.off.weights <-(1/as.double(c))*loopvals
	trade.off.weights


###
# Perform simulation
###

posveclong<-getempMVtotpos(3000,200,meanvec1,sigmamat,sigma0^2,mulim,c,R0)

# Plot portfolio weights for max-exp problem
plot(posveclong[1,],posveclong[2,],xlab="",ylab="",xlim=c(min(posveclong[1,]),max(posveclong[1,])),ylim=c(min(posveclong[2,]),max(posveclong[2,])))

# Theoretical weights with red dot
points(max.exp.weights[1],max.exp.weights[2],col="red",pch=19)

# Plot portfolio weights for min-var problem
plot(posveclong[3,],posveclong[4,],xlab="",ylab="",xlim=c(min(posveclong[3,]),max(posveclong[3,])),ylim=c(min(posveclong[4,]),max(posveclong[4,])))
# Theoretical weights with red dot
points(max.exp.weights[1],max.exp.weights[2],col="red",pch=19)

# Plot portfolio weights for trade-off problem
plot(posveclong[5,],posveclong[6,],xlab="",ylab="",xlim=c(min(posveclong[5,]),max(posveclong[5,])),ylim=c(min(posveclong[6,]),max(posveclong[6,])))
# Theoretical weights with red dot
points(max.exp.weights[1],max.exp.weights[2],col="red",pch=19)

###
# Compute the true mean and standard deviation of the portfolios
###

sdexpveclong1<-getempMVexpsds(posveclong[1:2,],meanvec1,sigmamat,R0)
sdexpveclong2<-getempMVexpsds(posveclong[3:4,],meanvec1,sigmamat,R0)
sdexpveclong3<-getempMVexpsds(posveclong[5:6,],meanvec1,sigmamat,R0)
# Efficient frontier
true.mean.sigma<-getempMVexpsds(max.exp.weights,meanvec1,sigmamat,R0)

# Plot the mean and standard deviation for max-exp
plot(sdexpveclong1[1,],sdexpveclong1[2,],xlab="",ylab="",xlim=c(min(sdexpveclong1[1,]),max(sdexpveclong1[1,])),ylim=c(min(sdexpveclong1[2,]),max(sdexpveclong1[2,])))
# Plot the true mean and std
points(true.mean.sigma[1],true.mean.sigma[2],col="red",pch=19)

# For min var
plot(sdexpveclong2[1,],sdexpveclong2[2,],xlab="",ylab="",xlim=c(min(sdexpveclong2[1,]),max(sdexpveclong2[1,])),ylim=c(min(sdexpveclong2[2,]),max(sdexpveclong2[2,])))
# Plot the true mean and std
points(true.mean.sigma[1],true.mean.sigma[2],col="red",pch=19)


# For trade-off
plot(sdexpveclong3[1,],sdexpveclong3[2,],xlab="",ylab="",xlim=c(min(sdexpveclong3[1,]),max(sdexpveclong3[1,])),ylim=c(min(sdexpveclong3[2,]),max(sdexpveclong3[2,])))
# Plot the true mean and std
points(true.mean.sigma[1],true.mean.sigma[2],col="red",pch=19)