# GLPK/Portfolio Optimization

This example is taken from investment theory, a branch of financial economics. The aim is to achieve the best tradeoff between return and an investor's tolerance for risk when allocating capital among a set of potential investments with known historical performance.

## Contents

## Summary[edit]

This tutorial example shows how to use GLPK/GMPL to optimize investment portfolios where investment risk is measured with the Mean Absolute Deviation (MAD) criterion. In addition to portfolio optimization, the example highlights techniques that may be useful in other GMPL applications, including

- Cholesky decomposition of positive definite matrices.
- Generating samples from a multivariate Normal distribution for the simulation of historical returns.
- Seeding GLPK's random number generator from GMPL.
*Note*: a new`randseed`option is currently under evaluation.

Additional examples by the same author are available at Introduction to Operations Research.

## Portfolio Optimization for Mean Absolute Deviation (MAD)[edit]

Portfolio optimization refers to the process of allocating capital among a set of financial assets to achieve a desired tradeoff between risk and return. The classical Markowitz approach to this problem measures risk using the expected variance of the portfolio return. This criterion yields a quadratic program for the relative weights of assets in the optimal portfolio.

In 1991, Konno and Yamazaki ^{[1]} proposed a linear programming model for portfolio optimization whereby risk is measured by the mean absolute deviation (MAD) from the expected return. Using MAD as the risk metric produces portfolios with several desirable properties not shared by the Markowitz portfolio, including second order stochastic dominance.

As originally formulated by Konno and Yamazaki, one starts with a history of returns for every asset in a set **Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "/mathoid/local/v1/":): S**
of assets. The return at time is determined by the change in price of the asset, **Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "/mathoid/local/v1/":): {\displaystyle R_i(t_n)= {(P_i(t_n)-P_i(t_{n-1}))}/{P_i(t_{n-1})}}**
. For each asset, the expected return is estimated by

The investor specifies a minimum required return . The portfolio optimization problem is to determine the fraction of the total investment allocated to each asset, **Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "/mathoid/local/v1/":): {\displaystyle w_i}**
, that minimizes the mean absolution deviation from the mean

subject to the required return and a fixed total investment:

The value of the minimum required return, **Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "/mathoid/local/v1/":): {\displaystyle R_p}**
expresses the investor's risk tolerance. A smaller value for increases the feasible solution space, resulting in portfolios with lower values of the MAD risk metric. Increasing results in portfolios with higher risk. The relationship between risk and return is a fundamental principle of investment theory.

This formulation doesn't place individual bounds on the weights **Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "/mathoid/local/v1/":): {\displaystyle w_i}**
. In particular, corresponds to short selling of the associated asset. Constraints can be added if the investor imposes limits on short-selling, on the maximum amount that can be invested in a single asset, or other requirements for portfolio diversification. Depending on the set of available investment opportunities, additional constraints can lead to infeasible investment problems.

## Implementation Details[edit]

The purpose of this example is to demonstrate the modeling of portfolio optimization problems in GMPL. The main tasks are to recast the mean absolute deviation objective in linear form, and to simulate historical returns given the mean and covariances for a set of potential investments.

### Reformulation of the MAD Objective[edit]

The following formulation of the objective function is adapted from "Optimization Methods in Finance" by Gerald Curnuejols and Reha Tutuncu (2007) ^{[2]} which, in turn, follows Feinstein and Thapa (1993) ^{[3]}. The model is streamlined by introducing decision variables and , so that the objective becomes

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "/mathoid/local/v1/":): {\displaystyle \min_{w_i, y_n,z_n} \frac{1}{N}\sum_{n = 1}^{N} (y_n+z_n)}**

subject to constraints

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "/mathoid/local/v1/":): {\displaystyle \begin{array}{rcl} \sum_{i \in S} w_i\bar{R}_i & \geq & R_{p} \\ \\ \sum_{i \in S} w_i & = & 1 \\ \\ y_n-z_n & = & \sum_{i \in S} w_i(R_i(t_n)-\bar{R}_i) \quad \mbox{for}\quad n = 1,\ldots,N \end{array} }**

As discussed by Feinstein and Thapa, this version reduces the problem to constraints in decision variables, noting the card function.

### Seeding the GLPK Pseudo-Random Number Generator[edit]

Unfortunately, GMPL does not provide a method to seed the pseudo-random number generator in GLPK. Instead, the following GMPL code fragment uses the function gmtime() to find the number of seconds since 1970. Dropping the leading digit avoids subsequent overflow errors, and the square returns a number with big changes every second. Extracting the lowest digits produces a number between 0 and 100,000 that determines how many times to sample GLPK's pseudo-random number generator prior to subsequent use. This hack comes with no assurances regarding its statistical properties.

/* A workaround for the lack of a way to seed the PRNG in GMPL */ param utc := prod {1..2} (gmtime()-1000000000); param seed := utc - 100000*floor(utc/100000); check sum{1..seed} Uniform01() > 0;

### Simulation of the Historical Returns[edit]

For this implementation, historical returns are simulated assuming knowledge of the means and covariances of asset returns. We begin with a vector of mean returns and covariance matrix estimated by

**Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "/mathoid/local/v1/":): {\displaystyle \Sigma_{ij} \approx \frac{1}{N-1}\sum_{n=1}^N(R_i(t_n)-\bar{R}_i)(R_j(t_n)-\bar{R}_j) }**

Simulation of historical returns requires generation of samples from a multivariable normal distribution. For this purpose we compute the Cholesky factorization where, for a positive semi-definite **Failed to parse (MathML with SVG or PNG fallback (recommended for modern browsers and accessibility tools): Invalid response ("Math extension cannot connect to Restbase.") from server "/mathoid/local/v1/":): {\displaystyle \Sigma}**
, and is a lower triangular matrix. The following GMPL code fragment implements the Cholesky-Crout algorithm.

/* Cholesky Lower Triangular Decomposition of the Covariance Matrix */ param C{i in S, j in S : i >= j} := if i = j then sqrt(Sigma[i,i]-(sum {k in S : k < i} (C[i,k]*C[i,k]))) else (Sigma[i,j]-sum{k in S : k < j} C[i,k]*C[j,k])/C[j,j];

Without error checking, this code fragment fails unless is positive definite. The covariance matrix is normally positive definite for real-world data, so this is generally not an issue. However, it does become an issue if one attempts to include a risk-free asset, such as a government bond, in the set of investment assets.

Once the Cholesky factor has been computed, a vector of simulated returns is given by where the elements of are independent samples from a normal distribution with zero mean and unit variance.

/* Simulated returns */ param N default 5000; set T := 1..N; param R{i in S, t in T} := Rbar[i] + sum {j in S : j <= i} C[i,j]*Normal(0,1);

## PortfolioMAD.mod[edit]

Save this model as `PortfolioMAD.mod`.

/* Portfolio Optimization using Mean Absolute Deviation Jeff Kantor December 4, 2009 Revised: December 6, 2009 to fix problem with random variate generation. Revised: December 7, 2009 to add a 'seeding' of the PRNG Revised: July 8, 2010 reformatted for GLPK Wikibook */ /* Stock Data */ set S; # Set of stocks param Rbar{S}; # Means of projected returns param Sigma{S,S}; # Covariance of projected returns param Rp default (1/card(S))*sum{i in S} Rbar[i]; # Lower bound on portfolio return /* Generate sample data */ /* Cholesky Lower Triangular Decomposition of the Covariance Matrix */ param C{i in S, j in S : i >= j} := if i = j then sqrt(Sigma[i,i]-(sum {k in S : k < i} (C[i,k]*C[i,k]))) else (Sigma[i,j]-sum{k in S : k < j} C[i,k]*C[j,k])/C[j,j]; /* A workaround for the lack of a way to seed the PRNG in GMPL */ param utc := prod {1..2} (gmtime()-1000000000); param seed := utc - 100000*floor(utc/100000); check sum{1..seed} Uniform01() > 0; /* Simulated returns */ param N default 5000; set T := 1..N; param R{i in S, t in T} := Rbar[i] + sum {j in S : j <= i} C[i,j]*Normal(0,1); /* MAD Optimization */ var w{S}; # Portfolio Weights with Bounds var y{T} >= 0; # Positive deviations (non-negative) var z{T} >= 0; # Negative deviations (non-negative) minimize MAD: (1/card(T))*sum {t in T} (y[t] + z[t]); s.t. C1: sum {s in S} w[s]*Rbar[s] >= Rp; s.t. C2: sum {s in S} w[s] = 1; s.t. C3 {t in T}: (y[t] - z[t]) = sum{s in S} (R[s,t]-Rbar[s])*w[s]; solve; /* Report */ /* Input Data */ printf "\n\nStock Data\n\n"; printf " Return Variance\n"; printf {i in S} "%5s %7.4f %7.4f\n", i, Rbar[i], Sigma[i,i]; printf "\nCovariance Matrix\n\n"; printf " "; printf {j in S} " %7s ", j; printf "\n"; for {i in S} { printf "%5s " ,i; printf {j in S} " %7.4f ", Sigma[i,j]; printf "\n"; } /* MAD Optimal Portfolio */ printf "\n\nMinimum Absolute Deviation (MAD) Portfolio\n\n"; printf " Return = %7.4f\n",Rp; printf " Variance = %7.4f\n\n", sum {i in S, j in S} w[i]*w[j]*Sigma[i,j]; printf " Weight\n"; printf {s in S} "%5s %7.4f\n", s, w[s]; printf "\n"; data; /* Data for monthly returns on four selected stocks for a three year period ending December 4, 2009 */ # Simulation Horizon param N := 5000; # Minimum acceptable investment return param Rp := 0.01; # Historical returns on assets param : S : Rbar := AAPL 0.0308 GE -0.0120 GS 0.0027 XOM 0.0018 ; # Covariance on asset returns param Sigma : AAPL GE GS XOM := AAPL 0.0158 0.0062 0.0088 0.0022 GE 0.0062 0.0136 0.0064 0.0011 GS 0.0088 0.0064 0.0135 0.0008 XOM 0.0022 0.0011 0.0008 0.0022 ; end;

## Results[edit]

Typical output follows. The results show that a well designed portfolio can exhibit a substantially improved risk/return performance compared to the risk/return performance of individual assets.

Stock Data Return Variance AAPL 0.0308 0.0158 GE -0.0120 0.0136 GS 0.0027 0.0135 XOM 0.0018 0.0022 Covariance Matrix AAPL GE GS XOM AAPL 0.0158 0.0062 0.0088 0.0022 GE 0.0062 0.0136 0.0064 0.0011 GS 0.0088 0.0064 0.0135 0.0008 XOM 0.0022 0.0011 0.0008 0.0022 Minimum Absolute Deviation (MAD) Portfolio Return = 0.0100 Variance = 0.0036 Weight AAPL 0.2794 GE 0.0002 GS 0.1120 XOM 0.6084

## Possible Extensions[edit]

There are number of extensions to this example that would make it more useful in real world settings:

- Add error checking for the Cholesky factorization.
- Add upper and lower (for example, no short selling) bound constraints on weighting coefficients.
- Add the ability to specify risk-free assets for the portfolio.
- Add constraints to impose diversification requirements.
- Add parametric analysis for the risk/return tradeoff (will require external scripting)
- Develop a better means of seeding GLPK's random number generator from within GMPL.
- Add an ODBC database interface for empirical (rather than simulated) historical returns data.

## References[edit]

- ↑ Konno, Hiroshi; Yamazaki, Hiroaki (1991). "Mean-absolute deviation portfolio optimization model and its applications to Tokyo stock market".
*Management Science***37**: 519-531. - ↑ Curnuejols, Gerald; Tutuncu, Reha (2007).
*Optimization Methods in Finance*. Cambridge University Press. - ↑ Feinstein, Charles D.; Thapa, Mukund N. (1993). "A Reformulation of a Mean-Absolute Deviation Portfolio Optimization Model".
*Management Science***39**: 1552-1553.