Introducing QuantLib: Portfolio Optimization

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 EditionSpecificallythe 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.

efnscalc

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:

notlinear

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:

noshortsales

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!

Advertisements

About Mick Hittesdorf

Financial Systems Architect, Analyst and Developer
This entry was posted in QuantLib and tagged , , , , , , , , , , , . Bookmark the permalink.

5 Responses to Introducing QuantLib: Portfolio Optimization

  1. E.L. Wisty says:

    Reblogged this on Pink Iguana.

  2. timson75 says:

    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

  3. Michel says:

    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

Leave a Reply

Fill in your details below or click an icon to log in:

WordPress.com Logo

You are commenting using your WordPress.com account. Log Out / Change )

Twitter picture

You are commenting using your Twitter account. Log Out / Change )

Facebook photo

You are commenting using your Facebook account. Log Out / Change )

Google+ photo

You are commenting using your Google+ account. Log Out / Change )

Connecting to %s