Practical Portfolio Optimization
K V Fernando
NAG Ltd
Wilkinson House
Jordan Hill
Oxford OX2 8DR
United Kingdom

email: vince@nag.co.uk

i
Abstract
NAG Libraries have many powerful and reliable optimizers which
can be used to solve large portfolio optimization and selection problems
in the ﬁnancial industry. These versatile routines are also suitable foracademic research and teaching.

Key words
Markowitz, mean-variance analysis, optimal portfolios, minimum
variance portfolio, portfolio selection, portfolio allocation, portfolio
diversiﬁcation, portfolio optimization, eﬃcient frontier,mean-variance frontier, MV eﬃciency

iiContents
1 Introduction

12 NAG Routines for Optimization
2.1 A Selection of Library Routines . . . . . . . . . . . . . . . . .
2.2 Quadratic Programming with Linear Constraints . . . . . . .
2.3 Nonlinear Programming . . . . . . . . . . . . . . . . . . . . .
2.4 Routines for Sparse Matrix Problems . . . . . . . . . . . . . .
2.5 Forward and Reverse Communication . . . . . . . . . . . . .
2.6 Hardware, Operating Systems and Environments . . . . . . .

2
2
3
3
3
4
4

3 Interfaces to Routines
3.1 Portfolio Weights . . . . . . . . . . . . . . . . . . . . . . . . .
3.2 Primary Data . . . . . . . . . . . . . . . . . . . . . . . . . . .
3.3 General Linear Constraints . . . . . . . . . . . . . . . . . . .
3.4 Nonlinear Constraints . . . . . . . . . . . . . . . . . . . . . .
3.5 Cold and Warm Starts . . . . . . . . . . . . . . . . . . . . . .

4
4
4
5
6
6

4 The Optimization Problems

6

5 Processing of Raw Data
5.1 The Covariance Matrix in Factored Form . . . . . . . . . . .
5.2 Determination of the Singular Values of the Cholesky Factors
5.3 If the Covariance Matrix Already Exists . . . . . . . . . . . .
5.4 Eigenvalues of the Covariance Matrix . . . . . . . . . . . . . .
5.5 Missing Values . . . . . . . . . . . . . . . . . . . . . . . . . .

7
7
9
9
10
10

6 Numerical Examples: Selection of Equities

10

7 Numerical Example: Asset Allocation

15

8 Transactions Costs

15

9 An Example Program

17

10 Acknowledgements

20

1

1

Introduction

The selection of assets or equities is not just a problem of ﬁnding attractive
investments. Designing the correct portfolio of assets cannot be done byhuman intuition alone and requires modern, powerful and reliable mathematical programs called optimizers. The Numerical Algorithms Group Ltd(NAG) is world renowned for its work on numerical algorithms, and NAG
routines for optimization are being used extensively in industry, commerce
and academia. Many leading ﬁnancial companies and institutions employ
NAG optimizers to select, diversify and rebalance their portfolios. They arealso used by business and management schools for teaching and research.
Any investor would like to have the highest return possible from an investment. However, this has to be counterbalanced by the amount of risk
the investor is able or desires to take. The expected return and the risk measured by the variance (or the standard deviation, which is the square-root ofthe variance) are the two main characteristics of a portfolio. Unfortunately,
equities with high returns usually correlate with high risk.
The behaviour of a portfolio can be quite diﬀerent from the behaviour ofindividual components of the portfolio. The risk of a properly constructed
portfolio from equities in leading markets could be half the sum of the risks ofindividual assets in the portfolio. This is due to complex correlation patterns
between individual assets or equities. A good optimizer can exploit thecorrelations, the expected returns, the risk (variance) and user constraints
to obtain an optimized portfolio. NAG optimization routines can deliver
optimized and diversiﬁed portfolios to match investor expectations.
The mathematical problem of portfolio optimization was initiated byProfessor Harry Markowitz in the ﬁfties and he was rewarded with a Nobel
Prize in Economics in 1990 which he shared with Professors William Sharpeand Merton Miller [8]. NAG optimizers can handle the classical Markowitz
optimization problems [7], [9], [10] and many modern day extensions [4],[11], [13], [14], [15]. NAG also provides a consultancy service to the ﬁnancial sector to solve mathematical, numerical, programming and visualization
problems associated with portfolio optimization.
Portfolio optimization is often called mean-variance (MV) optimization.
The term mean refers to the mean or the expected return of the investment
and the variance is the measure of the risk associated with the portfolio. The
mathematical problem can be be formulated in many ways but the principal
problems can be summarized as follows.
1. Minimize risk for a speciﬁed expected return
2. Maximize the expected return for a speciﬁed risk
3. Minimize the risk and maximize the expected return using a speciﬁedrisk aversion factor

2
4. Minimize the risk regardless of the expected return
5. Maximize the expected return regardless of the risk
6. Minimize the expected return regardless of the risk
The above problems could have linear or nonlinear constraints, equality andinequality constraints.
The ﬁrst three problems are essentially mathematically equivalent andthe solutions are called mean-variance (MV) eﬃcient. The eﬃcient points
in the Return-Risk graph are called the Eﬃcient Frontier.
The fourth problem gives minimum variance solutions which are for cautious investors. It is also used for comparison and benchmarking of otherportfolios. The ﬁfth problem gives the upper bound of the expected return
which can be attained; this is also useful for comparisons. The last problem
indicates a worst case scenario.
When market conditions or the investors risk preferences change, it isadvisable to rebalance the portfolio. Any of the above problems can besolved relative to an existing portfolio or a benchmark.
The transactions costs associated with purchasing a new portfolio or
rebalancing a portfolio could represent a signiﬁcant amount to the investor.
NAG optimization routines can handle transactions costs and they could
signiﬁcantly aﬀect the composition of the portfolio.

2

NAG Routines for Optimization

2.1

A Selection of Library Routines

NAG Libraries devote a complete Chapter to function optimization. Table 1
summarizes the routines available in NAG Libraries.
Routine
name
e04nc
e04nf
e04nk
e04uc
e04uf
e04ug

Constraints
linearlinearlinearnonlinear
nonlinear
nonlinear

Objectivefunction
quadratic
quadratic
quadratic
nonlinear
nonlinear
nonlinear

Dense orsparsedense
dense
sparsedense
dense
sparseTable 1: Optimization routines

3

2.2

Quadratic Programming with Linear Constraints

The e04nf routine is mainly for Linear Programming (LP) and Quadratic
Programming (QP), and the following objective functions can be minimized.
LP1: ct x. This option can be used for ﬁnding the maximum or the minimum
expected return regardless of the variance. The vector c is set tothe expected return vector a. The weight vector of the portfolio isrepresented by x.
QP1: 12 xt Hx. This option can be used to minimize the variance; the matrix
H is set to the covariance matrix V .
QP2: 12 xt Hx + ct x. This option can be used to minimize the variance fora speciﬁed risk aversion factor; the matrix H is set to the covariancematrix V and c is set to −λat x where λ is the risk aversion parameterand a is the expected return vector. Sometimes, the objective function
is deﬁned as 12 (xt V x − µat x); in that case, the new risk aversion factorµ is equal to 2λ.
QP3: 12 xt H t Hx. This is similar to QP1 except that the variance is supplied
via the Cholesky factor H of the covariance matrix V . See Section 5.1
for more details.
QP4: 12 xt H t Hx + ct x. This is similar to QP2 except that the covariance is
supplied via the Cholesky factor H of the covariance matrix V .
The routine e04nc can also handle least squares minimization or regression
problems. Many portfolio optimization problems can be formulated as leastsquares or regression problems. Hence the same routine e04nc can be used
for solving general portfolio optimization as well as portfolio problems which
are posed as least squares or regression problems.

2.3

Nonlinear Programming

The more general routines e04uc and e04uf can accommodate nonlinear
objective functions with nonlinear constraints.

2.4

Routines for Sparse Matrix Problems

Most of the optimization routines in NAG Libraries are for dense matrix
problems; that is, the constraint matrix is assumed to have mostly non-zeroelements. However, the routines e04nk and e04ug are designed to work withsparse constraint matrices. The routine e04nk is for linear and quadratic
programming and the routine e04ug is for nonlinear programming. The
routine e04nk does not explicitly require the covariance matrix V but theuser has to supply the matrix-vector product V x.

4

2.5

Forward and Reverse Communication

Most of the optimization routines are based on forward communications.
In such programs, the routine is called only once to obtain the results andthe user supplies all the necessary information to the NAG routine via asubroutine.
However, in some circumstances, it is necessary to do the optimization
step by step and call the user routine repeatedly to get fresh information.
The NAG routine e04uc is a forward communication routine and e04uf isthe reverse communication equivalent. This reverse communication routine
is particularly useful when it is called from another language (i.e., Microsoft
VBA) which does not fully support procedure arguments in a way that iscompatible with NAG routines.

2.6

Hardware, Operating Systems and Environments

NAG Libraries are available for all popular hardware systems and operating
systems. There are implementations for PCs, workstations, parallel computers such as SGI/Cray T3E and IBM SP3, networks of workstations, andpowerful SMP machines. Various operating systems from Windows (2000,
98, NT), DOS and Linux for PCs to UNIX and VMS for larger systems aresupported.
All of the optimization routines are available as DLLs.

3

Interfaces to Routines

3.1

Portfolio Weights

The weight or the proportion by value of the holding of the asset i is denotedby xi .

3.2

Primary Data

NAG optimization routines require the following details about the portfolio.
n number of assets in the portfolio
a expected returns of the assets in a vector of length ns volatilities (standard deviations) of the assets in a vector of length n
C correlation coeﬃcients in an n by n symmetric matrix
li lower limit of the weight of asset iui upper limit of the weight of asset i

5However, instead of the volatility vector s and the correlation matrix C, thefollowing may be available:
V covariance matrix in an n by n symmetric matrix
Usually, NAG optimization routines require the covariance matrix V , andif it is not available then it can be easily computed from the standard deviation vector s and the correlation matrix C. However, in Section 5.1,
we advocate that the covariance matrix V should not be formed explicitly
to avoid information loss; what is required is the Cholesky factor R of thecovariance matrix V and not the matrix V itself.

3.3

General Linear Constraints

The primary equality constraint is that all the weights should add up to a
constant, usually to unity. That isn

xi = 1i=1and the above is known as the budget constraint or the investment constraint. In matrix notation
et x = 1
where the column vector x is the weight vector which holds the proportionsof the assets
x = (x1 , x2 , . . . , xn )tand e is a column vector with all elements equal to unity.
e = (1, 1, . . . , 1)t .
The expected return vector is denoted bya = (a1 , a2 , . . . , an )t
where ai contains the the expected return of the asset or the equity i.
NAG optimizers can handle many equality constraints, and inequality
constraints and they are of the form
l≤

x

≤u

L ≤ Ax ≤ U
where

l

=

(l1 , . . . , ln )tu=(u1 , . . . , un )tL=(L1 , . . . , Lm )t

U=(U1 , . . . , Um )t6li and Li are lower limits, ui and Ui represent the upper limits, and mdenotes the number rows of the matrix A. For equality constraints, lower
and upper limits are set to equal values. It is also possible to set the upper
limits to +∞ and the lower limits to −∞These constraints could represent sector, industry, country, exposure(momentum, QRQ, P/B) and other user-deﬁned ﬁnancial constraints.

3.4

Nonlinear Constraints

NAG optimization routines, in particular routines e04uc, e04uf and e04ug,
can handle nonlinear constraints; they also can support nonlinear objectivefunctions.

3.5

Cold and Warm Starts

Cold starts refer to solutions of the problem from scratch. However, if theroutines are called repeatedly then approximate solutions are available fromprevious solutions. In that case, the initial conditions for the next iteration
may be supplied from the previous. Such son-called warm start facilities areavailable for many NAG optimization routines.

4

The Optimization Problems

Here we give some of the problems which can be solved using NAG optimizers. However, this list is not comprehensive and many more optimization problems can be solved using NAG software. For the mathematicallyminded, these problems can be stated as follows.
Problem 1 (Minimize the Risk)
minimize the variance xt V x
with a speciﬁed expected return r = at x subjected to linear and (or nonlinear)constraints.
Problem 2 (Maximize the Expected Return)maximize the expected return r = at x
with a speciﬁed variance ν = xt V x. This is also equivalent to minimizationof −r.
Problem 3 (Maximize Expected Return with Risk Aversion)
maximize λat x − xt V x
subjected to linear and (or nonlinear) constraints where λ is the risk aversion
parameter. This is also equivalent tominimize xt V x − λat x.

7
Problem 4 (Minimize Risk)
minimize the variance xt V x
subjected to linear and (or nonlinear) constraints.
Problem 5 (Maximize the Expected Return)maximize the expected return r = at x
subjected to linear and (or nonlinear) constraints. This is also equivalent tominimize − at x.
Problem 6 (Minimize the Expected Return)minimize the expected return r = at x
subjected to linear and (or nonlinear) constraints.
The above six problems can be solved with respect to a benchmark or the
existing portfolio. The following is an example and many more variants are
possible.
Problem 7 (Minimize Risk with Respect to a Benchmark)minimize (x − xb )t V (x − xb )with a speciﬁed expected relative return r deﬁned byr = at (x − xb )subjected to linear and (or nonlinear) constraints where xb represents abenchmark portfolio.

5

Processing of Raw Data

5.1

The Covariance Matrix in Factored Form

Let the element yi,j of the matrix Y denotes the return of the asset j attime i. If there are m time periods and n assets then the matrix Y is an mby n matrix. The sample mean for asset j is given bym
1 
yi,j .
aj =m i=1

Once the sample mean is removed from the returns, we get the deviationfrom the sample mean. This deviation Ŷ can be obtained asŶ = Y − eat8where e is a column vector of length m with each and every element equal
to unity,e = (1, . . . , 1)t .
The sample covariance matrix V is usually computed as
V =

1 t
Ŷ Ŷ
m

but some statisticians might prefer the formulaV =1
Ŷ t Ŷ .
m−1

However, it is well known to numerical analysts that such computation ofcovariance matrices lead to loss of information; see [1], [5] or [6].
Here is a well known example where information is lost in the formationof the covariance. Let Ŷ be




1 1 1
 θ



Ŷ = 
θθwhich is full rank provided that θ is non-zero. The covariance matrix is thengiven by

1 + θ21
1
1

1
1
1 + θ2Ŷ t Ŷ = 
.
3
2
1
1
1+θ
In IEEE machines with double precision ﬂoating point arithmetic the ma√chine precision  is approximately 2−53 . If θ is smaller than , then in
ﬂoating point arithmetic, the computer will determine the diagonal elementsof Ŷ t Ŷ as
1 + θ 2 −→ 1
due to rounding errors. In that case, the covariance matrix is evaluated as




1 1 1
1
1

t
Ŷ Ŷ =  1 1 1  = eet
3
3
1 1 1
which is a rank-one matrix and hence singular. Thus the explicit formulationof the covariance matrix should be avoided. One of the ways to avoid suchinformation loss is to compute the QR factorization of Ŷ , which is deﬁned
by
QR = Ŷ
where the m by n matrix Q contains orthogonal columns (i.e., Qt Q = I) and
R is an n by n upper triangular matrix. There are routines in NAG Libraries

9(e.g., f08ae or f01qc) to compute the QR factorization. The covariancematrix is then given by
V =

1
1 t t
R Q QR = Rt R.
mm

However, the above product does not have to be computed explicitly in√
portfolio optimization, and the Cholesky factor R/ m can be given as an
input to optimization routines. In this case, the optimization option QP3
or QP4 (see Section 2.2) should be used. Note that it is not necessary tocompute the matrix Q explicitly.

5.2

Determination of the Singular Values of the Cholesky
Factors

To check whether the covariance matrix is positive deﬁnite, the singularvalues of the Cholesky factor R may be computed. The eigenvalues of thecovariance matrix H are the squares of the singular values of the Cholesky
factors. Thus the eigenvalues of the covariance matrix may be determined
without explicitly forming the covariance matrix.

5.3

If the Covariance Matrix Already Exists

We have already indicated in a previous section that the covariance matrix
should not be formed explicitly. However, if it has been already formed thenit can be given as an input to NAG optimization routines directly using theoptions QP1 and QP2 (see Section 2.2 for further details). However, if NAG
optimizations are to be called more than once then it is more eﬃcient to use
QP3 and QP4 options. In that case, the Cholesky factors of the covariancematrixV = Rt R
should be determined prior to calling the optimization routines. This canbe done using f07fd or f03ae if the covariance matrix is positive deﬁnite. If
the factorization does not exist then these routines will indicate the failure
via an appropriate error handling ﬂag.
Covariance matrices are usually positive deﬁnite or at least positive semideﬁnite. However, due to rounding errors covariance matrices might not be
positive deﬁnite, and in that case the the factorization routines might fail
to compute the factors. However, there is no guarantee that the covariancematrix is positive deﬁnite even if the Cholesky factors exist. Consider theexample



1 1
1 0
1 1
=
= Rt R.
H=
1 1
1 0
0 0
The covariance matrix H has a Cholesky factor R but it is not positivedeﬁnite.

10

5.4

Eigenvalues of the Covariance Matrix

We have already indicated that to avoid information loss the covariance matrix should not be formed explicitly. However, if it has been already formed
then the best way to ﬁnd whether the covariance matrix is positive deﬁnite
or not is to compute the eigenvalues of the covariance matrix. If the eigenvalues are positive then the matrix is positive deﬁnite. However, even if all
eigenvalues are positive, it is advisable to compute the condition number ofthe matrix (which is deﬁned as the largest eigenvalue divided by the smallest eigenvalue). If the condition number is high then the matrix is nearly
positive semi-deﬁnite and this could lead to results which are unreliable.
The eigenvalues may be computed using f02aa but there are also manyother routines that could be used.

5.5

Missing Values

In Section 5.1, a matrix Y was deﬁned, and it contains the data to computethe covariance matrix and the Cholesky factors of the covariance matrix.
In practice, this matrix may contain missing (unrecorded) data. If this is
the case then a suitable data interpolation routine has to be used to assign
some values to the missing entries. More sophisticated users might ﬁnd
the missing values using the Brownian bridge EM algorithm [12] which isan adaptation of the classical EM algorithm [3]) for Brownian processes inﬁnance.

6

Numerical Examples: Selection of Equities

For illustrative purposes, we give a few examples of the results which could
be expected from NAG optimizers. The data was taken from [2] where
portfolios of shares in the Hang Seng (Hong Kong, 31 shares), DAX 100(Germany, 85 shares), FTSE 100 (UK, 89 shares), S & P 100 (USA, 98
shares) and the Nikkei (Japan, 225 shares) were studied using NAG e04
routines. Our eﬃcient frontiers are based on the routine e04nf.

The data refers to weekly returns and weekly variance. Figures 1 to 5
give the eﬃcient frontiers for this data. The dominant (red) curve denotes
the frontier if the weights are constrained between zero and unity (no shortsales are allowed). When tighter constraints are forced on the weights (0.0 ≤xi ≤ 0.1), a diﬀerent frontier appears. In general, with more constraints,
the expected returns become smaller for a speciﬁed risk value. Individual
points denote individual shares; the risk reduction due to optimization canbe judged by comparing these shares with the eﬃcient frontier in theseﬁgures.

11

0.012

0.01

Expected Return

0.008

0.006

0.004

0.002

0
0

0.001

0.002
0.003
Risk (variance)

0.004

0.005

Figure 1: The Eﬃcient Frontiers for the Hang Seng example

0.01

0.008

0.006

Expected Return

0.004

0.002

0

-0.002-0.004-0.006
0

0.001

0.002
0.003
Risk (variance)

0.004

0.005

Figure 2: The Eﬃcient Frontiers for the DAX 100 example

12

0.009

0.008

0.007

0.006

Expected Return

0.005

0.004

0.003

0.002

0.001

0

-0.001-0.002
0

0.001

0.002

0.003

Risk (variance)

Figure 3: The Eﬃcient Frontiers for the FTSE 100 example

0.01

0.008

Expected Return

0.006

0.004

0.002

0

-0.002
0

0.001

0.002
0.003
Risk (variance)

0.004

Figure 4: The Eﬃcient Frontiers for the S & P 100 example

13

0.004

0.002

Expected Return0-0.002-0.004-0.006-0.008

-0.01
0

0.001

0.002

0.003
0.004
Risk (variance)

0.005

0.006

Figure 5: The Eﬃcient Frontiers for the Nikkei example

In practical problems, with more ﬁnancial constraints, the eﬃcient frontier reduces to a smaller curve or almost to a point.
For the Nikkei example, we have also computed the Ineﬃcient Frontier,which is obtained by using negative values for the risk aversion parameter
λ; see Problem 3 in Section 4. In Figure 6, the top half of the curve (in red)shows the Eﬃcient Frontier and the bottom part of the curve (in blue) gives
the Ineﬃcient Frontier. No short sales were allowed in this example.
Markowitz-type optimization is frequently used for asset allocation. The
asset manager could have equities, government and corporate bonds, international equities and bonds, real estate, venture capital and gold to selectthe portfolio of assets. Each asset could itself be a portfolio (with or without
optimization).
For our next example, we have used eight types of assets (US 3 monthtreasury bills, US government long bonds, SP 500, Wilshire 500, Nasdeq
composite corporate bond index, EAFE and Gold). The yearly returns arefrom 1973 to 1994. These could be found at
http://www.princeton.edu/~rvdb/ampl/nlmodels/markowitz/
Figure 7 shows two eﬃcient frontiers. The dominant (red) curve shows
the frontier with upper limits of assets set to unity and the lower limits tozero. When the lower limit is raised to 0.05 and the upper limit is lowered
to 0.5, we get the eﬃcient frontier shown with a dotted curve (in blue).

14

0.004

0.002

Expected Return0-0.002-0.004-0.006-0.008

-0.01
0

0.001

0.002

0.003
0.004
Risk (variance)

0.005

0.006

Figure 6: The Eﬃcient and Ineﬃcient Frontiers for the Nikkei example

1.15

1.14

1.13

Expected Return

1.12

1.11

1.1

1.09

1.08

1.07
0

0.02

0.04

0.06
0.08
Risk (variance)

0.1

0.12

0.14

Figure 7: The Eﬃcient Frontiers in asset allocation

15

7

Numerical Example: Asset Allocation

In this numerical example, we started with raw data (yearly returns) anddid not form the covariance matrix explicitly. Instead, as described in Section 5.1, only the Cholesky factor R of the covariance was computed. In
this case, the optimization option QP3 or QP4 (see Section 2.2) should be
used since only the Cholesky factor R and not the whole covariance matrix
V has to be given as an input. Figure 7 exhibits the eﬃcient frontier forthis problem.

8

Transactions Costs

In the classical work of Markowitz, transactions costs associated with buying
and selling of equities were not allowed. However, now the importance ofincorporating transactions costs in building new portfolios and also in rebalancing existing portfolios are well recognized. In general, transactions costs
are not trivial enough to be neglected, and the optimal portfolio depends
upon the total cost of transactions.
Let the transactions cost associated with buying the equity i be pi andthen it may be modelled as(xi − x̄i )gi
0

pi =

for
for

xi > x̄i
xi ≤ x̄iwhere xi is the new portfolio weight of the equity i, x̄i is the original weightof the equity i, and gi is a constant associated with buying the equity i.
Similarly, let qi be the cost of selling the equity i,(x̄i − xi )ii
0

qi =

for
for

xi < x̄i
xi ≥ x̄iwhere the constant hi is associated with selling. Note that both pi and qi
cannot be simultaneously non-zero since the buying cost is zero when there
is selling and vice versa.
Let φ(x) be the objective function for minimization without transaction
costs. The new objective function with transaction costs is then given by
φ(x) +n

pi + qi = φ(x) +

i=1

n


max{pi , qi }.

i=1

The above objective function is not smooth (non-diﬀerentiable). However,this can be transformed into a smooth problem by including a new variable
yi for each equity i. The new objective is then given by
φ(x) +

n

i=1

yi

16
0.012

0.01

Expected Return

0.008

0.006

0.004

0.002

0
0

0.001

0.002
0.003
Risk (variance)

0.004

0.005

Figure 8: The Eﬃcient Frontiers with and without Transactions Costs

subjected to the constraintspi ≤ y iandqi ≤ yifor i = 1, . . . , n, where yi is the transaction cost for equity i. The above two
inequalities can be written in the format−∞ ≤ gi xi − yi

≤ gi x̄ihi x̄i ≤ yi + hi xi ≤ +∞Unfortunately, this approach doubles the number of variables for optimization.
We have recomputed the eﬃcient frontier for the Hang Seng example byincluding transactions costs. We have assumed that pi (which determines
the buying costs) is equal to 0.0001 and qi (which determines the selling
costs) is equal to 0.0002 for each equity i. Furthermore, it was assumed thatthe initial portfolio had an equal weighting, xi = 1/38, for each equity i.
Figure 8 shows the eﬃcient frontiers with and without transactions costs.
For a given variance, the diﬀerence between the expected returns (for a given
variance) reﬂects the transactions costs.

17

9

An Example Program

It is fairly easy to use NAG routines to do practical portfolio optimization.
The following example routine shows how to compute the optimal portfolio
when the expected return is minimized (i.e., the problem 1 in Section 1).
This program may be run in conduction with data used in [2]). In addition,
the user has to specify the expected return of the portfolio. The program
outputs the optimal weights and the variance of the portfolio.
/* Example Program for Portfolio Optimization** Copyright 2000 Numerical Algorithms Group Ltd.
** This example program computes the optimal weights* of the equities for a specified expected return* of the portfolio */#include <Nag/nag.h>#include <stdio.h>#include <Nag/nag_stdlib.h>#include <Nag/nag_string.h>#include <Nag/nage04.h>
static void ex1(void);static void ex1();#define MAXN 50#define MAXLIN 2#define MAXBND MAXN+MAXLIN
main(){
Vprintf("Portfolio optimization using e04nfc.\n");ex1();
exit(EXIT_SUCCESS);}
static void ex1(){
double x[MAXN], cvec[MAXN], s[MAXN];
double a[MAXLIN][MAXN], h[MAXN][MAXN];
double bl[MAXBND], bu[MAXBND];
double objf,up,low,ret,tem,tex,hh;
Integer tda, tdh;

18
Integer i, j, n, nclin, nbnd, ii, jj;
Boolean print;
Nag_E04_Opt options;
static NagError fail, fail2;
Vscanf(" %*[^\n]"); /* Skip heading in data file */Vscanf(" %*[^\n]"); /* Skip heading in data file */fail.print = TRUE;
fail2.print = TRUE;/* Read the actual problem dimension.
* n = the number of equities. */Vscanf("%ld",&n);/* nclin = the number of general linear constraints .
nclin = 2;
nbnd = n + nclin;
tda = MAXN;
tdh = MAXN;
/* a
* bl
* bu
* x
* s

= the linear constraint matrix A.
= the lower bounds on x and A*x.
= the upper bounds on x and A*x.
= the initial estimate of the solution.
= the standard deviation of return

For simplicity we assume thatlower bound on x is low
upper bound for x is up
for all x
Read low and up */Vscanf(" %*[^\n]"); /* Skip heading in data file */Vscanf("%lf%lf",&low,&up);/*

Set the initial weights to 1/n
set bl and bu
set the second row of A to unity
set the linear constraints bl[n+1] and bu[n+1] */for (i = 0; i < n; ++i) {x[i] = 1.0/n;*/

19
bl[i] = low;bu[i] = up;a[1][i] = 1.0;}
bl[n+1] = 1.0;
bu[n+1] = 1.0;/* Read the the expected return and sd
The expected return is in first column of A */Vscanf(" %*[^\n]"); /* Skip heading in data file */for (i = 0; i < n; ++i) {Vscanf("%lf%lf",&a[0][i],&s[i]);}/* Read the correlations
and form the covariance matrix */Vscanf(" %*[^\n]"); /* Skip heading in data file */for (i = 0; i < n; ++i)for (j = i; j < n; ++j) {Vscanf("%ld%ld%lf",&ii,&jj,&hh);h[ii-1][jj-1] = hh*s[i]*s[j];h[jj-1][ii-1] = hh*s[i]*s[j];}/* Read the required expected return of the portfolioand set that in bl[n+1] and bu[n+1] */Vscanf(" %*[^\n]"); /* Skip heading in data file */Vscanf("%lf",&ret);
bl[n] = ret;bu[n] = ret;
e04xxc(&options); /* Initialise options structure *//* Set one option directly* Bounds >=inf_bound will be treated as plus infinity.
* Bounds <= -inf_bound will be treated as minus infinity.
options.inf_bound = 1.0e21;/* Obtain remaining settings */fail.print = TRUE;*/

20
print = TRUE;
e04xyc("e04nfc", "stdin", &options, print, "stdout", &fail);/* Set the problem typeand other optional settings*/

options.prob = Nag_QP1;/* Solve the problem from a cold start.

*/

if (fail.code == NE_NOERROR)
e04nfc(n, nclin, (double *)a, tda, bl, bu, cvec, (double *)h, tdh,
NULLFN, x, &objf, &options, NAGCOMM_NULL, &fail);
Vprintf("Expected Return= %lf \n",ret);
Vprintf("Expected Variance = %lf \n",2.0*objf);/* Free memory allocated by e04nfc to pointers in options */e04xzc(&options, "all", &fail2);
if (fail.code != NE_NOERROR || fail2.code != NE_NOERROR) exit(EXIT_FAILURE);} /* ex1 */.

10

Acknowledgements

Discussions with Anne Trefethen, David Sayers, George Levy and Neil Swindells
are gratefully acknowledged.

References[1] A. Björk. Numerical Methods for Least Squares. SIAM, Philadelphia,
1996.
[2] T. J. Chang, N. Meade, J. E. Beasley, and Y. M. Sharaiha. Heuristics for cardinality constrained portfolio optimization. Computers andOperations Research, 2000, to appear.
[3] A. P. Dempster, N. M. Laired, and D. B. Rubin. Maximum likelihoodfrom incomplete data via the EM algorithm (with discussion). Journal of the Royal Statistical Society, Series B (Methodological), 39:1–38,
1977.
[4] E. J. Elton and M. J. Gruber. Modern Portfolio Theory and Investment
Analysis. Wiley, New York, 1995.

21[5] G. H. Golub and C. F. Van Loan. Matrix Computations. Johns Hopkins
University Press, Baltimore, 1996.
[6] C. L. Lawson and R. J. Hanson. Solving Least Squares Problem. SIAM,
Philadelphia, 1995.
[7] H. Markowitz. Portfolio selection. Journal of Finance, 7:77–91, 1952.
[8] H. Markowitz, W. F. Sharpe, and M. Miller. Founders of Modern Finance: Their Prize Winning Concepts and 1990 Nobel Lectures. AIMR,
Charlottesville VA, 1991.
[9] H. M. Markowitz. Mean-Variance Analysis in Portfolio Choice andCapital Markets. Blackwell, Oxford, 1987.
[10] H. M. Markowitz. Portfolio Selection: Eﬃcient Diversiﬁcation of Investments. Blackwell, Oxford, 1991.
[11] R. O. Michaud. Eﬃcient Asset Management. Harvard Business School
Press, Boston, 1998.
[12] W. Morokoﬀ. The Brownian bridge e-m algorithm for covaraince estimation with missing data. Journal of Computational Finance, 2:75–100,
1998.
[13] W. F. Sharpe. Portfolio Theory and Capital Markets. McGraw-Hill,
New York, 1999.
[14] W. F. Sharpe, G. J. Alexander, and J. V. Bailey. Investments. Prentice
Hall, Upper Saddle River, NJ, 1998.
[15] S. A. Zenios. Financial Optimization. Cambridge University Press,
Cambridge, 1993.