Portfolio Optimization using Mean Absolute Deviation
Contents
Portfolio Optimization using Mean Absolute Deviation#
This IPython notebook demonstrates portfolio optimization using the Mean Absolute Deviation (MAD) criterion. A portion of these notes is adapted from GLPK Wikibook tutorial on the subject written by me.
Background#
Konno and Yamazaki (1990) proposed a linear programming model for portfolio optimization in which the risk measure is mean absolute deviation (MAD). This model computes a portfolio minimizing MAD subject to a lower bound on return.
In contrast to the classical Markowitz portfolio, the MAD criterion requires a data set consisting of returns on the investment assets. The data set may be an historical record or samples from a multivariate statistical model of portfolio returns. The MAD criterion produces portfolios with properties not shared by the Markowitz portfolio, including second degree stochastic dominance.
Below we demonstrate portfolio optimization with the MAD criterion where data is generated by sampling a multivariate normal distribution. Given mean return \(r\) and the Cholesky decomposition of the covariance matrix \(\Sigma\) (i.e., \(C\) such that \(CC^T = \Sigma\)), we compute \(r_t = r + C z_t\) where the elements of \(z_t\) are zero mean normal variates with unit variance.
The rest of the formulation is adapted from “Optimization Methods in Finance” by Gerald Curnuejols and Reha Tutuncu (2007) which, in turn, follows an implementation due to Fienstein and Thapa (1993).
A complete tutorial describing the implementation of this model is available on the GLPK wikibook.
Mean Absolute Deviation#
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 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 [[w:Stochastic_dominance|second order stochastic dominance]].
As originally formulated by Konno and Yamazaki, one starts with a history of returns \(R_i(t_n)\) for every asset in a set \(S\) of assets. The return at time \(t_n\) is determined by the change in price of the asset,
For each asset, the expected return is estimated by
The investor specifies a minimum required return \(R_{p}\). The portfolio optimization problem is to determine the fraction of the total investment allocated to each asset, \(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, \(R_p\) expresses the investor’s risk tolerance. A smaller value for \(R_p\) increases the feasible solution space, resulting in portfolios with lower values of the MAD risk metric. Increasing \(R_p\) 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 \(w_i\). In particular, \(w_i < 0\) 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.
Reformulation of the MAD Objective#
The following formulation of the objective function is adapted from “Optimization Methods in Finance” by Gerald Curnuejols and Reha Tutuncu (2007) which, in turn, follows Feinstein and Thapa (1993) .
The model is streamlined by introducing decision variables \(y_n \geq 0\) and \(z_n \geq 0\), \(n=1,\ldots,N\) so that the objective becomes
subject to constraints
As discussed by Feinstein and Thapa, this version reduces the problem to \(N+2\) constraints in \( 2N+\mbox{card}(S)\) decision variables, noting the card function.
Seeding the GLPK Pseudo-Random Number Generator#
Unfortunately, MathProg 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#
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 \(\bar{R}\) and covariance matrix \(\Sigma\) estimated by
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 \(\Sigma\), \(\Sigma=CC^T\) and \(C\) is a lower triangular matrix. The following MathProg 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 \(\Sigma\) 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 \(C\) has been computed, a vector of simulated returns \(R(t_n)\) is given by \(R(t_n) = \bar{R} + C Z(t_n)\) where the elements of \(Z(t_n)\) 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);
MathProg Model#
%%script glpsol -m /dev/stdin
# Example: PortfolioMAD.mod Portfolio Optimization using Mean Absolute Deviation
/* Stock Data */
set S; # Set of stocks
param r{S}; # Means of projected returns
param cov{S,S}; # Covariance of projected returns
param r_portfolio
default (1/card(S))*sum{i in S} r[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(cov[i,i]-(sum {k in S : k < i} (c[i,k]*c[i,k])))
else
(cov[i,j]-sum{k in S : k < j} c[i,k]*c[j,k])/c[j,j];
/* Because there is no way to seed the PRNG, a workaround */
param utc := prod {1..2} (gmtime()-1000000000);
param seed := utc - 100000*floor(utc/100000);
check sum{1..seed} Uniform01() > 0;
/* Normal random variates */
param N default 5000;
set T := 1..N;
param zn{j in S, t in T} := Normal(0,1);
param rt{i in S, t in T} := r[i] + sum {j in S : j <= i} c[i,j]*zn[j,t];
/* MAD Optimization */
var w{S} >= 0; # 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]*r[s] >= r_portfolio;
s.t. C2: sum {s in S} w[s] = 1;
s.t. C3 {t in T}: (y[t] - z[t]) = sum{s in S} (rt[s,t]-r[s])*w[s];
solve;
/* Report */
/* Input Data */
printf "Stock Data\n\n";
printf " Return Variance\n";
printf {i in S} "%5s %7.4f %7.4f\n", i, r[i], cov[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 ", cov[i,j];
printf "\n";
}
/* MAD Optimal Portfolio */
printf "\nMinimum Absolute Deviation (MAD) Portfolio\n\n";
printf " Return = %7.4f\n",r_portfolio;
printf " Variance = %7.4f\n\n", sum {i in S, j in S} w[i]*w[j]*cov[i,j];
printf " Weight\n";
printf {s in S} "%5s %7.4f\n", s, w[s];
printf "\n";
table tab0 {s in S} OUT "JSON" "Optimal Portfolio" "PieChart":
s, w[s]~PortfolioWeight;
table tab1 {s in S} OUT "JSON" "Asset Return versus Volatility" "ScatterChart":
sqrt(cov[s,s])~StDev, r[s]~Return;
table tab2 {s in S} OUT "JSON" "Portfolio Weights" "ColumnChart":
s~Stock, w[s]~PortfolioWeight;
table tab3 {t in T} OUT "JSON" "Simulated Portfolio Return" "LineChart":
t~month, (y[t] - z[t])~PortfolioReturn;
/* Simulated Return data in Matlab Format */
/*
printf "\nrt = [ ... \n";
for {t in T} {
printf {s in S} "%9.4f",rt[s,t];
printf "; ...\n";
}
printf "];\n\n";
*/
data;
/* Data for monthly returns on four selected stocks for a three
year period ending December 4, 2009 */
param N := 200;
param r_portfolio := 0.01;
param : S : r :=
AAPL 0.0308
GE -0.0120
GS 0.0027
XOM 0.0018 ;
param cov :
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;
GLPSOL: GLPK LP/MIP Solver, v4.52
Parameter(s) specified in the command line:
-m /dev/stdin
Reading model section from /dev/stdin...
Reading data section from /dev/stdin...
/dev/stdin:115: warning: final NL missing before end of file
115 lines were read
Checking (line 24)...
Generating MAD...
Generating C1...
Generating C2...
Generating C3...
Model has been successfully generated
GLPK Simplex Optimizer, v4.52
203 rows, 404 columns, 1608 non-zeros
Preprocessing...
140 rows, 142 columns, 698 non-zeros
Scaling...
A: min|aij| = 2.301e-04 max|aij| = 1.000e+00 ratio = 4.346e+03
GM: min|aij| = 9.139e-02 max|aij| = 1.094e+01 ratio = 1.197e+02
EQ: min|aij| = 8.625e-03 max|aij| = 1.000e+00 ratio = 1.159e+02
Constructing initial basis...
Size of triangular part is 140
0: obj = 1.502038445e-02 infeas = 7.709e+00 (0)
* 116: obj = 6.534903273e-02 infeas = 0.000e+00 (0)
* 157: obj = 4.760242137e-02 infeas = 0.000e+00 (0)
OPTIMAL LP SOLUTION FOUND
Time used: 0.0 secs
Memory used: 1.0 Mb (1076387 bytes)
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.0033
Weight
AAPL 0.2828
GE 0.0000
GS 0.0000
XOM 0.7172
Writing tab0...
Invalid table driver `JSON'
/dev/stdin:71: error on opening table tab0
Model postsolving error