#'@title Implementation of PSR Portfolio Optimization
#'@description
#'Maximizing for PSR leads to better diversified and more balanced hedge fund allocations compared to the concentrated
#'outcomes of Sharpe ratio maximization.We would like to find the vector of weights that maximize the expression
#'
#'\deqn{\hat{PSR}(SR^{*}) = Z\bigg[\frac{(\hat{SR}-SR^{*})\sqrt{n-1}}{\sqrt{1-\hat{\gamma_3}SR^{*} + \frac{\hat{\gamma_4}-1}{4}\hat{SR^2}}}\bigg]}
#'
#'where \eqn{\sigma = \sqrt{E[(r-\mu)^2]}} ,its standard deviation.\eqn{\gamma_3=\frac{E\biggl[(r-\mu)^3\biggr]}{\sigma^3}} its skewness,
#'\eqn{\gamma_4=\frac{E\bigg[(r-\mu)^4\bigg]}{\sigma^4}} its kurtosis and \eqn{SR = \frac{\mu}{\sigma}} its Sharpe Ratio.
#'Because \eqn{\hat{PSR}(SR^{*})=Z[\hat{Z^{*}}]} is a monotonic increasing function of
#'\eqn{\hat{Z^{*}}} ,it suffices to compute the vector that maximizes \eqn{\hat{Z^{*}}}
#'
#'This optimal vector is invariant of the value adopted by the parameter \eqn{SR^{*}}.
#'Gradient Ascent Logic is used to compute the weights using the Function PsrPortfolio
#'@aliases PsrPortfolio
#'
#'@useDynLib noniid.pm
#'@param R The return series
#'@param refSR The benchmark Sharpe Ratio
#'@param bounds The bounds for the weights
#'@param MaxIter The Maximum number of iterations
#'@param delta The value of delta Z
#'
#'@author Pulkit Mehrotra
#'@seealso \code{\link{ProbSharpeRatio}} \code{\link{table.PSR}} \code{\link{MinTrackRecord}}
#'@references Bailey, David H. and Lopez de Prado, Marcos, \emph{The Sharpe Ratio
#'Efficient Frontier} (July 1, 2012). Journal of Risk, Vol. 15, No. 2, Winter
#'2012/13
#'
#'@keywords ts multivariate distribution models
#'
#'@examples
#'
#'data(edhec)
#'PsrPortfolio(edhec)
#'@export
PsrPortfolio<-function(R,refSR=0,bounds=NULL,MaxIter = 1000,delta = 0.005){
# DESCRIPTION:
# This function returns the weight for which Probabilistic Sharpe Ratio
# is maximized.
#
# INPUT:
# The return series of the portfolio is taken as the input
# refSR , bounds of the weights for each series in the portfolio ,
# max iteration for the optimization , delta by which the z value will
# change is also taken as the input.
#
# OUTPUT:
# The weights is given as the output.
x = checkData(R)
columns = ncol(x)
n = nrow(x)
columnnames = colnames(x)
k<-rep(0,5)
if(is.null(bounds)){
message("Bounds not given assuming bounds to be (0,1) for each weight")
bounds = matrix(rep(c(0,1),columns),nrow = columns,byrow = TRUE)
}
z = 0
iter = 0
w = rep(1,columns)/columns
d1z = 0
#Optimization Function
optimize<-function(){
weights = w
mean = NULL
for(column in 1:columns){
mean = c(mean,get_Moments(x[,column],1))
}
while(TRUE){
if(iter == MaxIter) break
dZ = get_d1Zs(mean,weights)
if(dZ$z>z && checkBounds(weights)==TRUE){
z = dZ$z
d1z = dZ$d1Z
w = weights
}
iter = iter + 1
weights = stepSize(weights,dZ$d1Z)
if(is.null(weights)) break
}
return(w)
}
# To Check the bounds of the weights
checkBounds<-function(weights){
flag = TRUE
for(i in 1:columns){
if(weights[i] < bounds[i,1]) flag = FALSE
if(weights[i] > bounds[i,2]) flag = FALSE
}
return(flag)
}
#Calculate the step size to change the weights
stepSize<-function(weights,d1Z){
if(length(which(d1Z!=0)) == 0){
return(NULL)
}
index = which(abs(d1Z) ==max(abs(d1Z)))
weights[index] = weights[index]+delta/d1Z[index]
weights = weights/sum(weights)
return(weights)
}
#To get the first differentials
get_d1Zs<-function(mean,weights){
d1Z = numeric(columns)
m = numeric(4)
x_portfolio = x%*%weights
m[1] = get_Moments(x_portfolio,1)
for(i in 2:4){
m[i] = get_Moments(x_portfolio,i,m[1])
}
stats = get_Stats(m)
#mu = mean(x_portfolio)
#sd = StdDev(x_portfolio)
#sk = skewness(x_portfolio)
#kr = kurtosis(x_portfolio)
#stats = c(mu,sd,sk,kr)
#m = c(stats[1],stats[2]^2,stats[3]*(stats[2]^3),stats[4]*(stats[2]^4))
SR = get_SR(stats,n)
meanSR = SR$meanSR
sigmaSR = SR$sigmaSR
if(refSR>meanSR){
stop("The Reference Sharpe Ratio should be less than the Observred Sharpe Ratio")
}
for(i in 1:columns){
d1Z[i] = get_d1Z(stats,m,meanSR,sigmaSR,mean,weights,i)
}
dZ = list("d1Z"=d1Z,"z"=(meanSR-refSR)/sigmaSR)
return(dZ)
}
#To get the dZ/dw for each weight
get_d1Z<-function(stats,m,meanSR,sigmaSR,mean,weights,index){
d1Mu = get_d1Mu(mean,index)
d1Sigma = get_d1Sigma(stats[2],mean,weights,index)
d1Skew = get_d1Skew(d1Sigma,stats[2],mean,weights,index,m[3])
d1Kurt = get_d1Kurt(d1Sigma,stats[2],mean,weights,index,m[4])
d1meanSR = (d1Mu*stats[2]-d1Sigma*stats[1])/stats[2]^2
d1sigmaSR = (d1Kurt * meanSR^2+2*meanSR*d1meanSR*(stats[4]-1))/4
d1sigmaSR = d1sigmaSR-(d1Skew*meanSR+d1meanSR*stats[3])
d1sigmaSR = d1sigmaSR/(2*sigmaSR*(n-1))
d1Z = (d1meanSR*sigmaSR-d1sigmaSR*(meanSR-refSR))/sigmaSR^2
return(d1Z)
}
get_d1Mu<-function(mean,index){
return(mean[index])
}
get_d1Sigma<-function(sigma,mean,weights,index){
return(get_dnMoments(mean,weights,2,1,index)/(2*sigma))
}
get_d1Skew<-function(d1Sigma,sigma,mean,weights,index,m3){
d1Skew = get_dnMoments(mean,weights,3,1,index)*sigma^3
d1Skew = d1Skew - 3*(sigma^2)*d1Sigma*m3
d1Skew = d1Skew/sigma^6
return(d1Skew)
}
get_d1Kurt<-function(d1Sigma,sigma,mean,weights,index,m4){
d1Kurt = get_dnMoments(mean,weights,4,1,index)*sigma^4
d1Kurt = d1Kurt - 4*(sigma^3)*d1Sigma*m4
d1Kurt = d1Kurt/sigma^8
return(d1Kurt)
}
#To get the differential of the moments
get_dnMoments<-function(mean,weights,mOrder,dOrder,index){
sum = 0
x0 = 1
for(i in 0:(dOrder-1)){
x0 = x0*(mOrder-i)
}
x_mat = as.matrix(na.omit(x))
sum = 0
output = .Call("sums",mat = x_mat,index,mean,dOrder,weights,mOrder,sum)
#for(i in 1:n){
# x1 = 0
# x2 = (x_mat[i,index]-mean[index])^dOrder
#if(index == 1){
# print(x2)
#}
# for(j in 1:columns){
# x1 = x1 + weights[j]*(x_mat[i,j]-mean[j])
# }
# sum = sum + x2*x1^(mOrder-dOrder)
# }
return(x0*(output)/n)
}
# TO get meanSR and sigmaSR
get_SR<-function(stats,n){
meanSR = stats[1]/stats[2]
sigmaSR = ((1-meanSR*stats[3]+(meanSR^2)*(stats[4]-1)/4)/(n-1))^0.5
SR<-list("meanSR"=meanSR,"sigmaSR"=sigmaSR)
return(SR)
}
get_Stats<-function(m){
return(c(m[1],m[2]^0.5,m[3]/(m[2]^1.5),m[4]/(m[2]^2)))
}
get_Moments<-function(series,order,mean = 0){
sum = 0
mat = na.omit(as.matrix(series))
sum = .Call("sums_m",mat,mean,order)
# for(i in series){
# sum = sum + (i-mean)^order
# }
return(sum/n)
}
weights = optimize()
result = matrix(weights,nrow = columns)
rownames(result) = columnnames
colnames(result) = "weight"
return(result)
}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.