Welcome back. If you read my last two posts, you may recall that I showed how to use QuantLib to do linear optimization in the context of a bond portfolio construction problem and mean-variance portfolio analysis by way of the concept of the Efficient Frontier, respectively. In this installment of my ‘Introducing QuantLib’ series, I will combine elements of each of these two posts in order to demonstrate how QuantLib can be used to generate the Efficient Frontier in the presence of constraints on a portfolio’s composition. Specifically, short sales will not be allowed. Finding the Efficient Frontier with no short sales will require the repeated solution of a non-linear portfolio optimization problem in which every asset in the portfolio must have a non-negative weighting.

As in previous posts, I will first model the problem in a spreadsheet and then replicate the solution in C++ using classes from the QuantLib library. The essentials of the model are borrowed from Chapter 11 of Simon Benninga’s excellent book *Financial Modeling, 2nd Edition. *Specifically*, *the weightings on each of four assets, AAPL, IBM, ORCL, GOOG, must be determined such that Sharpe Ratio of the portfolio is maximized subject to the following constraints:

- The portfolio’s asset weightings must sum to 1
- Every asset weighting in the portfolio must be greater than or equal to zero

A screen shot of the full spreadsheet is below, which can be downloaded from my Box account at https://www.box.com/s/7w4tvfqueqqu676vhz5l

As stated in my earlier post on linear optimization, all optimization problems feature a function that must be either maximized or minimized. This function is called the objective function. In this portfolio optimization problem, *Theta* is the value of the objective function that is to maximized. It is defined as the ratio of excess return to portfolio standard deviation, where excess return is the portfolio return over and above the risk-free rate, defined as the constant, c, in the spreadsheet. Theta, is then essentially, the same thing as the portfolio’s Sharpe Ratio for a given value of c. More intuitively, Theta is a measure of reward for a given amount of risk. In maximizing Theta, we seek to allocate our capital amongst the four stocks in the portfolio such that we accrue the greatest reward for the least amount of risk. This goal is consistent with the fact that most investors are risk-averse.

On a side note, you may have noticed that this spreadsheet was created using OpenOffice, not LibreOffice, which I employed in previous posts. I had to switch to OpenOffice in order to install the ‘Solver for Nonlinear Programming‘ OpenOffice extension, which is currently not supported by LibreOffice on my laptop’s operating system, Ubuntu 12.04. Without the non-linear solver extension, LibreOffice detects that the problem is non-linear and fails to find a solution:

Now let’s see how to set up and solve this problem in QuantLib. The C++ code is below:

```
#include <iostream>
#include <cstdlib>
#define BOOST_AUTO_TEST_MAIN
#include <boost/test/unit_test.hpp>
#include <boost/detail/lightweight_test.hpp>
#include <ql/quantlib.hpp>
#include <boost/format.hpp>
#include <functional>
#include <numeric>
#include <fstream>
namespace {
using namespace QuantLib;
//aggregates all constraint expressions into a single constraint
class PortfolioAllocationConstraints : public Constraint {
public:
PortfolioAllocationConstraints(const std::vector& expressions) : Constraint(boost::shared_ptr(new
PortfolioAllocationConstraints::Impl(expressions))) {}
private:
// constraint implementation
class Impl : public Constraint::Impl {
public:
Impl(const std::vector& expressions) :
expressions_(expressions) {}
bool test(const Array& x) const {
for (auto iter = expressions_.begin();
iter < expressions_.end(); ++iter) {
if (!(*iter)(x)) {
return false;
}
}
//will only get here if all constraints satisfied
return true;
}
private:
const std::vector expressions_;
};
};
class ThetaCostFunction: public CostFunction {
public:
ThetaCostFunction(const Matrix& covarianceMatrix,
const Matrix& returnMatrix) : covarianceMatrix_(covarianceMatrix),
returnMatrix_(returnMatrix) {
}
Real value(const Array& proportions) const {
QL_REQUIRE(proportions.size()==3, "Four assets in portfolio!");
Array allProportions(4);
allProportions[0] = proportions[0];
allProportions[1] = proportions[1];
allProportions[2] = proportions[2];
allProportions[3] = 1 - (proportions[0] + proportions[1] + proportions[2]);
return -1 * ((portfolioMean(allProportions) - c_)/portfolioStdDeviation(allProportions));
}
Disposable values(const Array& proportions) const {
QL_REQUIRE(proportions.size()==3, "Four assets in portfolio!");
Array values(1);
values[0] = value(proportions);
return values;
}
void setC(Real c) { c_ = c; }
Real getC() const { return c_;}
Real portfolioMean(const Array& proportions) const {
Real portfolioMean = (proportions * returnMatrix_)[0];
//std::cout << boost::format("Portfolio mean: %.4f") % portfolioMean << std::endl;
return portfolioMean;
}
Real portfolioStdDeviation(const Array& proportions) const {
Matrix matrixProportions(4,1);
for (size_t row = 0; row < 4; ++row) {
matrixProportions[row][0] = proportions[row];
}
const Matrix& portfolioVarianceMatrix = transpose(matrixProportions) * covarianceMatrix_ * matrixProportions;
Real portfolioVariance = portfolioVarianceMatrix[0][0];
Real stdDeviation = std::sqrt(portfolioVariance);
//std::cout << boost::format("Portfolio standard deviation: %.4f") % stdDeviation << std::endl;
return stdDeviation;
}
private:
const Matrix& covarianceMatrix_;
const Matrix& returnMatrix_;
Real c_;
};
BOOST_AUTO_TEST_CASE(testNoShortSales) {
Matrix covarianceMatrix(4,4);
//row 1
covarianceMatrix[0][0] = .1; //AAPL-AAPL
covarianceMatrix[0][1] = .03; //AAPL-IBM
covarianceMatrix[0][2] = -.08; //AAPL-ORCL
covarianceMatrix[0][3] = .05; //AAPL-GOOG
//row 2
covarianceMatrix[1][0] = .03; //IBM-AAPL
covarianceMatrix[1][1] = .20; //IBM-IBM
covarianceMatrix[1][2] = .02; //IBM-ORCL
covarianceMatrix[1][3] = .03; //IBM-GOOG
//row 3
covarianceMatrix[2][0] = -.08; //ORCL-AAPL
covarianceMatrix[2][1] = .02; //ORCL-IBM
covarianceMatrix[2][2] = .30; //ORCL-ORCL
covarianceMatrix[2][3] = .2; //ORCL-GOOG
//row 4
covarianceMatrix[3][0] = .05; //GOOG-AAPL
covarianceMatrix[3][1] = .03; //GOOG-IBM
covarianceMatrix[3][2] = .2; //GOOG-ORCL
covarianceMatrix[3][3] = .9; //GOOG-GOOG
std::cout << "Covariance matrix of returns: " << std::endl;
std::cout << covarianceMatrix << std::endl;
//portfolio return vector
Matrix portfolioReturnVector(4,1);
portfolioReturnVector[0][0] = .08; //AAPL
portfolioReturnVector[1][0] = .09; //IBM
portfolioReturnVector[2][0] = .10; //ORCL
portfolioReturnVector[3][0] = .11; //GOOG
std::cout << "Portfolio return vector" << std::endl;
std::cout << portfolioReturnVector << std::endl;
//constraints
std::vector<std::function > noShortSalesConstraints(2);
//constraints implemented as C++ 11 lambda expressions
noShortSalesConstraints[0] = [] (const Array& x) {Real x4 = 1.0 - ( x[0] + x[1] + x[2]); return (x[0] >= 0.0 && x[1] >= 0.0 && x[2] >= 0.0 && x4 >= 0.0);};
noShortSalesConstraints[1] = [] (const Array& x) { Real x4 = 1.0 - ( x[0] + x[1] + x[2]); return 1.0 - (x[0] + x[1] + x[2] + x4) < 1e-9;};
//instantiate constraints
PortfolioAllocationConstraints noShortSalesPortfolioConstraints(noShortSalesConstraints);
Size maxIterations = 100000;
Size minStatIterations = 100;
Real rootEpsilon = 1e-9;
Real functionEpsilon = 1e-9;
Real gradientNormEpsilon = 1e-9;
EndCriteria endCriteria (maxIterations, minStatIterations,
rootEpsilon, functionEpsilon, gradientNormEpsilon);
std::map<Rate, std::pair > mapOfStdDeviationToMeanNoShortSales;
Rate startingC = -.035;
Real increment = .005;
for (int i = 0; i < 40; ++i) {
Rate c = startingC + (i * increment);
ThetaCostFunction thetaCostFunction(covarianceMatrix, portfolioReturnVector);
thetaCostFunction.setC(c);
Problem efficientFrontierNoShortSalesProblem(thetaCostFunction, noShortSalesPortfolioConstraints, Array(3, .2500));
Simplex solver(.01);
EndCriteria::Type noShortSalesSolution = solver.minimize(efficientFrontierNoShortSalesProblem, endCriteria);
std::cout << boost::format("Solution type: %s") % noShortSalesSolution << std::endl;
const Array& results = efficientFrontierNoShortSalesProblem.currentValue();
Array proportions(4);
proportions[0] = results[0];
proportions[1] = results[1];
proportions[2] = results[2];
proportions[3] = 1.0 - (results[0] + results[1] + results[2]);
std::cout << boost::format("Constant (c): %.4f") % thetaCostFunction.getC() << std::endl;
std::cout << boost::format("AAPL weighting: %.4f") % proportions[0] << std::endl;
std::cout << boost::format("IBM weighting: %.4f") % proportions[1] << std::endl;
std::cout << boost::format("ORCL weighting: %.4f") % proportions[2] << std::endl;
std::cout << boost::format("GOOG weighting: %.4f") % proportions[3] << std::endl;
std::cout << boost::format("Theta: %.4f") % (-1 * efficientFrontierNoShortSalesProblem.functionValue()) << std::endl;
Real portfolioMean = thetaCostFunction.portfolioMean(proportions);
std::cout << boost::format("Portfolio mean: %.4f") % portfolioMean << std::endl;
Volatility portfolioStdDeviation = thetaCostFunction.portfolioStdDeviation(proportions);
std::cout << boost::format("Portfolio standard deviation: %.4f") % portfolioStdDeviation << std::endl;
mapOfStdDeviationToMeanNoShortSales[c] = std::make_pair(portfolioStdDeviation, portfolioMean);
std::cout << "------------------------------------" << std::endl;
}
//write efficient frontier with no short sales to file
std::ofstream noShortSalesFile;
noShortSalesFile.open("/tmp/noshortsales.dat", std::ios::out);
for (std::map<Rate, std::pair >::const_iterator i=mapOfStdDeviationToMeanNoShortSales.begin(); i != mapOfStdDeviationToMeanNoShortSales.end(); ++i) {
noShortSalesFile <first % i->second.first % i->second.second << std::endl;
}
noShortSalesFile.close();
//plot with gnuplot using commands below Run 'gnuplot' then type in:
/*
set terminal png
set output "/tmp/noshortsales.png"
set key top left
set key box
set xlabel "Volatility"
set ylabel "Expected Return"
plot '/tmp/noshortsales.dat' using 2:3 w linespoints title "No Short Sales"
*/
}
}
```

When run this code solves for the optimal weightings for all values of c between -.035 and .16. For example, the ouput for c = .05 is:

————————————

Solution type: StationaryPoint

Constant (c): 0.0500

AAPL weighting: 0.5659

IBM weighting: 0.1047

ORCL weighting: 0.3294

GOOG weighting: 0.0000

Theta: 0.1839

Portfolio mean: 0.0876

Portfolio standard deviation: 0.2046

————————————

This solution is in very close agreement with the solution arrived at by the OpenOffice spreadsheet with the non-linear solver extension. So, it looks like the QuantLib code is doing the right thing. Cool!

Towards the end of code listing, I provide the gnuplot commands that render the chart for the Efficient Frontier with no short sales, pictured below:

Lastly, I want to mention that I originally planned to use a QuantLib OptimizationMethod class other than Simplex, such as ConjugateGradient, because I mistakenly believed that QuantLib’s Simplex class could not handle non-linear optimization problems. However, when I looked at the comments for the Simplex class in the Simplex.hpp header, I learned that it is based on the multidimensional Nelder Mead Simplex method, which does support non-linear optimization in several variables.

That concludes this installment of my ‘Introducing QuantLib’ series. I hope you enjoyed it. In my next post, I will be steering away from portfolio theory and optimization, which have been the focus of my last three posts, and transitioning to a discussion of option pricing. Check back soon and, as always, I encourage you to submit comments and/or questions. In the meantime, have fun with QuantLib!

Reblogged this on Pink Iguana.

Mike,

Interesting post and thank you for posting.

I was experimenting with the Quant Lib optimizer and trying to do something similar. As an exercise I attempted to use XLW to create an XLL and provide several portfolio weights functions (e.g. GetWeightsERC(), GetWeightsMV(), GetWeightsMSR()).

I wanted to replicate the Equal Risk Contribution logic implemented in R on http://systematicinvestor.wordpress.com/ and discussed in http://www.thierry-roncalli.com/download/erc.pdf and then implement the analytical solution discussed in papers.ssrn.com/sol3/papers.cfm?abstract_id=2117303. Other implementations also discussed in http://www.wilmott.com/messageview.cfm?catid=34&threadid=38497.

Initially I attempted to create a simple function calculating the unconstrained MV portfolio both analytically and using the Simplex optimiser. In some cases my results differ from the analytical solution and native Excel solver solution and unsure why. Any suggestions or further references would be much appreciated.

Have attached relevant parts of code below.

Many thanks,

Tim

// mean variance objective function

class MeanVarianceFunction: public QuantLib::CostFunction{

private:

QuantLib::Matrix covariance_;

QuantLib::Size n;

public:

MeanVarianceFunction(const QuantLib::Matrix& covariance)

:covariance_(covariance)

{

n = covariance.rows();

}

QuantLib::Real value(const QuantLib::Array& x) const{

QL_REQUIRE(x.size()==n-1, “n – 1 weights required”);

QuantLib::Matrix weights(n,1);

QuantLib::Real sumWeights(0.);

for(QuantLib::Size i = 0; i < x.size(); ++i){

sumWeights += x[i];

weights[i][0] = x[i];

}

weights[n-1][0] = 1 – sumWeights;

QuantLib::Matrix var = transpose(weights) * covariance_ * weights;

return var[0][0];

}

QuantLib::Disposable values(const QuantLib::Array& x) const{

QuantLib::Array var(1,value(x));

return var;

}

};

class EqualRiskContributionFunction: public QuantLib::CostFunction{

private:

QuantLib::Matrix covariance_;

QuantLib::Size n;

public:

EqualRiskContributionFunction(const QuantLib::Matrix& covariance)

:covariance_(covariance)

{

n = covariance.rows();

}

QuantLib::Real value(const QuantLib::Array& x) const{

QuantLib::Array dsx = values(x);

QuantLib::Real res(0.);

for(QuantLib::Size i = 0; i < n; ++i){

for(QuantLib::Size j = 0; j < i; ++j){

res += (dsx[i] – dsx[j]) * (dsx[i] – dsx[j]);

}

}

return res;

}

QuantLib::Disposable values(const QuantLib::Array& x) const{

QuantLib::Array dsx(n,0.);

for(QuantLib::Size i = 0; i < n; ++i){

QuantLib::Array row(covariance_.row_begin(i), covariance_.row_end(i));

dsx[i] = QuantLib::DotProduct(row, x) * x[i];

double test = dsx[i];

}

return dsx;

}

};

// long only constraint

class LongOnlyConstraint : public QuantLib::Constraint {

private:

class Impl : public QuantLib::Constraint::Impl {

public:

Impl(){}

bool test(const QuantLib::Array& weights) const {

QuantLib::Real sumWeights(0.);

for (QuantLib::Size i=0; i<weights.size(); i++) {

if(weights[i] < 0 – GECKO::EPSILON)

return false;

sumWeights += weights[i];

}

// n'th asset weight equals 1 – sum(asset 1 …. assent n -1)

if(1-sumWeights < 0 – GECKO::EPSILON)

return false;

return true;

}

private:

};

public:

LongOnlyConstraint()

: QuantLib::Constraint(boost::shared_ptr(new LongOnlyConstraint::Impl)) {}

};

// fully invested constraint

class FullyInvestedConstraint : public QuantLib::Constraint {

private:

class Impl : public QuantLib::Constraint::Impl {

public:

Impl() {}

bool test(const QuantLib::Array& weights) const {

QuantLib::Real sumWeights(0.);

for (QuantLib::Size i=0; i<weights.size(); i++) {

sumWeights += weights[i];

}

QuantLib::Real weightN = 1. – sumWeights;

// check minimum leverage

if(1. – (sumWeights + weightN) < -GECKO::EPSILON)

return false;

return true;

}

private:

};

public:

FullyInvestedConstraint()

: QuantLib::Constraint(boost::shared_ptr(new FullyInvestedConstraint::Impl)) {}

};

Matrix Portfolio::WeightsMV(bool isConstrained)

{

Matrix weights(n,1);

if(isConstrained){

MeanVarianceFunction mvFunc(Cov());

QuantLib::Size maxIterations=10000; //end search after 1000 iterations if no solution

QuantLib::Size minStatIterations=10; //don’t spend more than 10 iterations at a single point

Real rootEpsilon=1e-10; //end search if absolute difference of current and last root value is below epsilon

Real functionEpsilon=1e-10; //end search if absolute difference of current and last function value is below epsilon

Real gradientNormEpsilon=1e-6; //end search if absolute difference of norm of current and last gradient is below epsilon

EndCriteria myEndCrit(maxIterations, minStatIterations, rootEpsilon,

functionEpsilon, gradientNormEpsilon);

//constraints

LongOnlyConstraint longOnly;

FullyInvestedConstraint fullyInvested;

CompositeConstraint allConstraints(longOnly, fullyInvested);

Problem myProb(mvFunc, allConstraints, Array(n-1, 1./n));

Simplex solver(.05);

EndCriteria::Type solution=solver.minimize(myProb, myEndCrit);

switch (solution) {

case EndCriteria::None:

case EndCriteria::MaxIterations:

case EndCriteria::Unknown:

throw(“#err: optimization didn’t converge – no solution found”);

default:

;

}

Array x = myProb.currentValue();

double sumWeights(0);

for(QuantLib::Size i = 0; i < x.size(); ++i){

weights[i][0] = x[i];

sumWeights += weights[i][0];

}

weights[n-1][0] = 1- sumWeights;

}

else {

//closed form solution exists for unconstrained

Matrix covInv = inverse(Cov());

Matrix l(n,1,1);

weights = (covInv * l) / ( transpose(l) * covInv * l )[0][0];

}

return weights;

}

Matrix Portfolio::WeightsERC()

{

Matrix weights(n,1);

if (ValidateIdenticalCorrelation(correlation)){

weights = WeightsInvVol();

}

else

{

EqualRiskContributionFunction ercFunc(Cov());

QuantLib::Size maxIterations=10000; //end search after 1000 iterations if no solution

QuantLib::Size minStatIterations=100; //don't spend more than 10 iterations at a single point

Real rootEpsilon=1e-10; //end search if absolute difference of current and last root value is below epsilon

Real functionEpsilon=1e-10; //end search if absolute difference of current and last function value is below epsilon

Real gradientNormEpsilon=1e-6; //end search if absolute difference of norm of current and last gradient is below epsilon

EndCriteria myEndCrit(maxIterations, minStatIterations, rootEpsilon,

functionEpsilon, gradientNormEpsilon);

//constraints

LongOnlyConstraint longOnly;

FullyInvestedConstraint fullyInvested;

CompositeConstraint allConstraints(longOnly, fullyInvested);

//would be better to start with variance weighted array but easier to debug and less inputs

Problem myProb(ercFunc, allConstraints, Array(n, 2./n));

LevenbergMarquardt solver;

EndCriteria::Type solution=solver.minimize(myProb, myEndCrit);

switch (solution) {

case EndCriteria::None:

case EndCriteria::MaxIterations:

case EndCriteria::Unknown:

throw("#err: optimization didn't converge – no solution found");

default:

;

}

Array x = myProb.currentValue();

for(QuantLib::Size i = 0; i < x.size(); ++i)

weights[i][0] = x[i];

}

return weights;

}

Tim,

Thanks for posting your question and the example source code.

My experience with the QuantLib Simplex solver is that it is sensitive to both the initial guess, which is the third argument to the Problem class constructor (e.g. Array(n-1, 1/n)) and lambda, the Simplex constructor parameter. You might try varying these inputs to see if the solution becomes more stable if you haven’t tried that already.

When I get a chance, I’ll compile your example code and investigate a bit more. I’ll follow up with anything I find.

Thanks for reading my blog and keep the questions/comments coming! Mick

Hi Mick,

Great article, unfortunately i’m afraid that the code presents some formatting issues (vectors missing template type). Could you please upload share the original .cpp file or send it to me?

Many thanks

Michel

Sorry. I have to markup the template parameters or they get misinterpreted in the HTML. All the code for my posts is available on GitHub at https://github.com/mhittesdorf/allthingsfintech. You should find what you are looking for there.

Thanks for reading my blog. I hope you found it helpful. Mick