arXiv:cs/0107007v1 [] 3 Jul 2001

Chapter 1
THE RISK PROFILE PROBLEM FOR
STOCK PORTFOLIO OPTIMIZATION
Ming-Yang Kao∗
Department of Computer Science
Yale University
New Haven, CT 06520
USA
kao-ming-yang@cs.yale.edu

Andreas Nolte
Department of Computer Science
Yale University
New Haven, CT 06520
USA
nolte@cs.yale.edu

Stephen R. Tate†
Department of Computer Science
University of North Texas
Denton, TX 76203
USA.
srt@cs.unt.edu

Abstract

This work initiates research into the problem of determining an optimal
investment strategy for investors with different attitudes towards the
trade-offs of risk and profit. The probability distribution of the return
values of the stocks that are considered by the investor are assumed to

∗ Research

supported in part by NSF grants CCR-9531028 and CCR-9974871. Current address: Department of Computer Science, Northwestern University, Evanston, IL 60201, USA.
† Research supported in part by Texas Advanced Research Program Grant 1997-003594-019.

1

The Risk Profile Problem forStock Portfolio Optimization

1

be known, while the joint distribution is unknown. The problem is to
find the best investment strategy in order to minimize the probability
of losing a certain percentage of the invested capital based on different
attitudes of the investors towards future outcomes of the stock market.
For portfolios made up of two stocks, this work shows how to exactly and quickly solve the problem of finding an optimal portfolio for
aggressive or risk-averse investors, using an algorithm based on a fast
greedy solution to a maximum flow problem. However, an investor looking for an average-case guarantee (so is neither aggressive or risk-averse)
must deal with a more difficult problem. In particular, it is ♯P -complete
to compute the distribution function associated with the average-case
bound. On the positive side, approximate answers can be computed by
using random sampling techniques similar to those for high-dimensional
volume estimation. When k > 2 stocks are considered, it is proved that
a simple solution based on the same flow concepts as the 2-stock algorithm would imply that P = N P , so is highly unlikely. This work gives
approximation algorithms for this case as well as exact algorithms for
some important special cases.
Keywords: risk management, portfolio optimization, computational hardness, approximation algorithms, greedy strategies, network flows, volume estimation, random walks.

1.

Introduction

This work initiates the study of the risk profile problem for stock
portfolio optimization. The problem has several variants depending on
a given investor’s preference toward the trade-off between risk and return
[Sharpe et al., 1995].
In the problem, the investor has a capital, which is normalized to
one dollar. She considers k different stocks S1 , . . . , Sk and wishes to
invest P
some xi dollars in each stock Si for a certain period of time,
where ki=1 xi = 1 and xi ≥ 0 for all i. The vector ~x = hxi iki=1 =
hx1 , x2 , . . . , xk i is called a portfolio. Let Pk be the set of all portfolios
for k stocks. The return of ~x is the ratio, expressed as a percentage, of
the worth of this portfolio at the end of the investment period to the
initial investment of one dollar. The return of stock Sj is the ratio of its
price at the end of the investment period to its initial price, which is the
same as the return of the portfolio hxi iki=1 with xj = 1 and all the other
xi = 0.
In mathematical finance, stock prices are often assumed to follow
geometric Brownian motions or its variants (e.g., see [Duffie, 1996, Elliott and Kopp, 1999, Fouque et al., 2000, Hull, 2000, Karatzas, 1997,
Karatzas and Shreve, 1998, Musiela and Rutkowski, 1997]). To comple-

2
ment this conventional approach with computer science methodologies
[Cormen et al., 1990], we assume that stock prices can move arbitrarily.
Let µ be a positive real number. Let m1 and m2 be integers with m1 <
m2 , and let m = m2 − m1 + 1. Let ∆ = {ℓµ | ℓ = m1 , . . . , m2 }. Each
stock Si is associated with a discrete probability distribution Si over ∆,
where Si (β) is the probability that the stock’s return is β%. For the
sake of technical convenience, we allow m1 and m2 to be negative. The
probability distributions S1 , . . . , Sk are part of the input in our problem
and are obtainable, e.g., by observing historical market data. We assume
that non-zero values satisfy S1 (β) ≥ 1/nc for some constant c, and
when representation is important we assume that these values can be
represented as fixed-point numbers with O(log n) bits. The parameters
µ, m1 , and m2 control the precision and range of such observations.
For instance, for µ = 1, m1 = 0, and m2 = 200, the set of possible
returns are 0%, 1%, . . . , 200%. The joint distribution of the k probability
distributions Si is usually unavailable for a variety of practical reasons.
In particular, a joint distribution consists of nk entries and thus would
require observing an exponential number of data points in k.
The investor’s goal is to find a portfolio ~x, which is optimal according
to her risk preference in six basic cases as follows. For a risk-averse
investor, minimizing loss is more important than maximizing win, while
an aggressive investor has the opposite priority. Each of these two investor types can be further classified into three subtypes, namely, bestcase, worst-case, and average-case, referring to whether the probability
of loss or win is estimated in the best, worst, or average case over the
feasible joint distributions. More precisely, for each of these six types,
the investor first chooses a target return α and then looks for such a
portfolio ~x that optimizes one of the following six probabilities:
RAb (α, ~x) (respectively, RAw (α, ~x) or RAa (α, ~x)) is the smallest
(respectively, largest or average) probability that the return of ~x
is at most α% over all joint distributions for S1 , . . . , Sk .
AGb (α, ~x) (respectively, AGw (α, ~x) or AGa (α, ~x)) is the largest (respectively, smallest or average) probability that the return of ~x is
at least α% over all joint distributions for S1 , . . . , Sk .
If the investor is best-case (respectively, worst-case or average-case)
risk-averse, she would choose ~x to minimize RAb (α, ~x) (respectively,
RAw (α, ~x) or RAa (α, ~x)). In contrast, if the investor is best-case (respectively, worst-case or average-case) aggressive, she would choose ~x to
maximize AGb (α, ~x) (respectively, AGw (α, ~x) or AGa (α, ~x)).
While the risk profile problem originates from a very applied field,
the corresponding mathematical model has a substantial combinatorial

The Risk Profile Problem forStock Portfolio Optimization

3

structure. In the cases where the investor is highly risk-averse or highly
aggressive, we can model the problem as a network flow problem. Quite
surprisingly, in the two-stock case, this flow problem is solvable by a simple greedy algorithm in O(m) time. In contrast, for the three-stock case,
the applicability of a greedy flow-based algorithm would imply P = N P .
If the number k of stocks is part of the input, we give an exact algorithm
based on linear programming which takes time polynomial in the number of entries of a corresponding contingency table but exponential in
the input size. To supplement this algorithm, we also give a polynomialtime approximation algorithm based on linear programming. We further
present an exact polynomial-time algorithm in the practical case where
the capital can only be broken up into a fixed number of units (e.g.,
cents).
It remains open whether this problem is N P -complete if the number
of stocks is part of the input. We strongly suspect that this is indeed
the case.
In the case of an average-case investor we show ♯P -hardness of the
problem of computing the distribution function over various probability
bounds, a natural first-step in solving the average-case investor problem.
This hardness result holds even in two dimensions, and we describe an
approximation algorithm for this case. This algorithm uses a random
walk approach to sample from the feasible joint distributions, and is
closely related to volume computation and sampling from log-concave
distributions.
Section 2 defines some notation. Section 3 discusses the case where
there are only two stocks under consideration. Section 4 discusses the
case of general k. Due to page limitations, all figures are placed in the
appendix (these figures are helpful in understanding the material, but
are not strictly necessary).

2.

Notation
Let ~δ ∈ ∆k denote a vector hδ1 , . . . , δk i, where δi ∈ ∆. Let
M = [M~δ ]~δ∈∆k

denote a k-dimensional matrix indexed by ∆k . Let Mk denote the set of
k-dimensional matrices for all possible joint distributions of S1 , . . . , Sk ;
i.e., Mk consists of all matrices
M = [M~δ ]~δ∈∆k ,
where (1) M~δ is the probability that the return of stock Si is δi % for
i = 1, . . . , k, and (2) thus for all ~δ ∈ ∆k , M~δ ≥ 0 and for all β ∈ ∆ and

4
j = 1, . . . , k,
X

Sj (β) =

M~δ .

~
δ∈∆k ;δj =β

For instance, Mk contains the matrix M defined by
M~δ =

k
Y

Si (δi ).

i=1

Also, in the two-stock case, each M ∈ M2 is just a two-dimensional
m × m matrix, where for all δ1 , δ2 ∈ ∆, the entries of M in column δ1
sum up to S1 (δ1 ) and those in row δ2 sum up to S2 (δ2 ).
Given a portfolio ~x ∈ Pk and a target return α, let
L(α, ~x) =

(

∗∗

(

L (α, ~x) =
U (α, ~x) =
U ∗∗ (α, ~x) =

~δ ∈ ∆ |
k

k
X

)

xi δi ≤ α ,

i=1

~δ ∈ ∆ |
k

(

~δ ∈ ∆k |

(

~δ ∈ ∆k |

k
X

i=1
k
X

i=1
k
X

)

xi δi < α ,
)

xi δi ≥ α ,
)

xi δi > α ,

i=1

which are the sets of the indices of all entries in the matrices in Mk such
that the return of ~x is at most, less than, at least, and more than α%,
respectively. We further define the following functions on M ∈ Mk :
X

Lα,~x (M ) =

M~δ ,

~
δ∈L(α,~
x)
∗∗
Lα,~x (M )

X

=

M~δ ,

~
δ∈L∗∗ (α,~
x)

X

U α,~x (M ) =

M~δ ,

~
δ∈U (α,~
x)
∗∗

U α,~x (M ) =

X

M~δ ,

~
δ∈U ∗∗ (α,~
x)

which are the probabilities in the joint distribution M that the return
of ~x is at most, less than, at least, and more than α%, respectively.
Formally, if uMk (M ) is a uniform density over Mk ,
RAb (α, ~x) =

min Lα,~x (M );

M ∈Mk

(1.1)

The Risk Profile Problem forStock Portfolio Optimization

RAw (α, ~x) =
RAa (α, ~x) =
AGb (α, ~x) =
AGw (α, ~x) =
AGa (α, ~x) =

max Lα,~x (M );

M ∈Mk

Z

Mk

Lα,~x (M )uMk (M )dM ;

5
(1.2)
(1.3)

max U α,~x (M );

(1.4)

min U α,~x (M );

(1.5)

M ∈Mk
M ∈Mk

Z

Mk

U α,~x (M )uMk (M )dM.

(1.6)

For example, in the two-stock case, L(α, hx1 , x2 i) is the set of all indices
in a two-dimensional table M in M2 on or below the line x1 δ1 +x2 δ2 = α,
and RAw (α, hx1 , x2 i) maximizes the sum of the entries in this region
under the condition that M has the given column and row sums of
S1 (m1 ), . . . , S1 (m2 ), S2 (m1 ), . . . , S2 (m2 ).
For technical convenience, we also define the following terms:
RA∗∗
x) =
b (α, ~
RA∗∗
x) =
w (α, ~
RA∗∗
x) =
a (α, ~
AGb∗∗ (α, ~x) =
∗∗
AGw
(α, ~x) =

AGa∗∗ (α, ~x) =

∗∗

(1.7)

∗∗

(1.8)

min Lα,~x (M );

M ∈Mk

max Lα,~x (M );

M ∈Mk

Z

∗∗

Mk

Lα,~x (M )dM ;
∗∗

(1.10)

∗∗

(1.11)

max U α,~x (M );

M ∈Mk

min U α,~x (M );

M ∈Mk

Z

Mk

(1.9)

∗∗

U α,~x (M )dM.

(1.12)

Lemma 1 The following statements hold.
min RAb (α, ~x) = 1 − max AGb∗∗ (α, ~x)

(1.13)

∗∗
(α, ~x)
min RAw (α, ~x) = 1 − max AGw

(1.14)

min RAa (α, ~x) = 1 − max AGa∗∗ (α, ~x)

(1.15)

x)
max AGb (α, ~x) = 1 − min RA∗∗
b (α, ~

(1.16)

x)
max AGw (α, ~x) = 1 − min RA∗∗
w (α, ~

(1.17)

x)
max AGa (α, ~x) = 1 − min RA∗∗
a (α, ~

(1.18)

~
x∈Pk
~
x∈Pk

x∈Pk
~

x∈Pk
~
~
x∈Pk

~
x∈Pk

~
x∈Pk

~
x∈Pk

~
x∈Pk

~
x∈Pk

~
x∈Pk

~
x∈Pk

6
Proof: Straightforward.
In light of Lemma 1, to solve the risk profile problem, it suffices to
show how to compute
min~x∈Pk RAb (α, ~x), min~x∈Pk RAw (α, ~x), min~x∈Pk RAa (α, ~x),
x).
x), min~x∈Pk RA∗∗
x), min~x∈Pk RA∗∗
min~x∈Pk RA∗∗
a (α, ~
w (α, ~
b (α, ~
The techniques for computing the latter three expressions are essentially the same as those for computing the former three. Furthermore,
the techniques for computing the first expression are almost identical
to those for computing the second. For these reasons, the remainder
of our discussion focuses on how to compute min~x∈Pk RAw (α, ~x) and
min~x∈Pk RAa (α, ~x).

3.

The Two-Stock Case

This section assumes that k = 2, i.e., there are only two stocks under
consideration. In the case of two stocks, we can visualize the problems
under consideration as in Figure 1.1. The discrete and finite set of
possible return pairs for the two stocks in the portfolio are shown as
the dots in this picture – each pair has a probability (from the joint
distribution) associated with it, with the given restrictions on column
and row sums. A given portfolio and target return α defines a half-space
on the set of return pairs, with the shaded area in Figure 1.1 giving
the area in which the total return is ≤ α. The problem of computing
RAw (α, ~x) then is the problem of determining which feasible assignment
of joint probabilities places the highest total probability in the shaded
region.

3.1.

A Worst-Case or Best-Case Investor

Given a target return α, this section focuses on how to compute an
optimal portfolio for a worst-case risk-averse investor. The cases of a
best-case risk-averse investor, a worst-case aggressive investor, and a
best-case aggressive investor can be solved similarly.
We first present a basic algorithm to compute RAw (α, ~x) by computing a worst-case joint distribution matrix M for S1 and S2 . For
convenience, we index the entries of M with {(i, j) | i, j = m1 , . . . , m2 },
where row i (respectively, column i) corresponds to return iµ of S1 (respectively, jµ of S2 ). We model the problem of computing M as a
network flow problem on the graph G defined below:
G has 2(m + 1) vertices, namely, a source s, a sink t, and
vm1 , . . . , vm2 , wm1 , . . . , wm2 , where vi (respectively, wi ) corresponds to return iµ of stock S1 (respectively, stock S2 ).

The Risk Profile Problem forStock Portfolio Optimization

Figure 1.1.

7

Visualization of two stock case

For all i, j = m1 , . . . , m2 , G has (1) edge (vi , wj ), which has capacity c(vi , wj ) = 1 if x1 iµ + x2 jµ ≤ α or 0 otherwise; (2) the edge
(s, vi ) with capacity c(s, vi ) = S1 (iµ); and (3) the edge (wj , t) with
capacity c(wj , t) = S2 (jµ).
Geometrically, we wish to push as much probability as possible into
the region of M defined by x1 i + x2 j ≤ αµ . In other words, the value of a
maximum s − t flow of G equals RAw (α, ~x). Thus, it is tempting to use
a maximum flow algorithm to solve this maximum flow problem. The
fastest known algorithm for this problem is due to Goldberg and Rao
2
[Goldberg and Rao, 1998] and runs in O ∗ (m2 3 ) time1 for our application
(note that m in this bound is as defined in this work, not as the number
of edges which is typical in general flow discussion). Instead of using this
algorithm, we exploit some structural properties of G to solve the flow
problem using a simple greedy algorithm in O(m) arithmetic operations.
Note that since G may have Ω(m2 ) edges with positive capacity, we
cannot afford to construct the whole G explicitly. The idea of our O(m)time algorithm can be described as follows.
Starting with vm2 , we try to push a flow of c(s, vm2 ) through
G. Assume c(vm2 , wm1 ) = 1 for simplicity. We consider the path
formed by edges (s, vm2 ), (vm2 , wm1 ), (wm1 , t) first. We can push flow
1 We use O ∗ (f (n)) for the “soft-O” notation, which ignores polylogarithmic factors. In bounds
for the approximation algorithms, this notation also ignores factors that depend only on the
approximation bound ǫ.

8
min(c(s, vm2 ), c(wm1 , t)) through this path, saturating either (s, vm2 )
or (wm1 , t). If we saturated (s, vm2 ) then we next consider the path
(s, vm2 −1 ), (vm2 −1 , wm1 ), (wm1 , t) for pushing additional flow; however, if we had saturated (wm1 , t) we will next consider the path
(s, vm2 ), (vm2 , wm1 +1 ), (wm1 +1 , t). We continue in this fashion until we
can push no more flow. The only complication is that if at some point we
are considering the path (s, vi ), (vi , wj ), (wj , t), and c(vi , wj ) = 0, then
obviously we can’t saturate either (s, vi ) or (vj , t), and we simply decrease i to next consider the path (s, vi−1 ), (vi−1 , wj ), (wj , t). The details
of this O(m) time algorithm are given in Figure 1.2.
procedure Greedy-Flow
F ←0
i ← m2
cv ← c(s, vi )
j ← m1
cw ← c(wj , t)
loop
if c(vi , wj ) = 1 and cw ≤ cv then
F ← F + cw
cv ← cv − cw
j ←j +1
if j > m2 then return F
cw ← c(wj , t)
else
if c(vi , wj ) = 1 then
F ← F + cv
cw ← cw − cv
end if
i←i−1
if i < m1 then return F
cv ← c(s, vi )
end if
end loop
Figure 1.2.

The procedure Greedy-Flow

Theorem 2 Given S1 , S2 , a valid portfolio vector ~x, and α as input,
Greedy-Flow computes the value of a maximum flow of G in O(m) arithmetic operations.

The Risk Profile Problem forStock Portfolio Optimization

9

Proof: As a first step we prove that the algorithm computes the
maximal flow. Let ℓ be the minimal index such that (wℓ , t) is not saturated after termination of the algorithm and k be the minimal index
such that c(vk , wℓ ) = 0. We define a partition V1 ∪ V2 of the nodes by
V1 = {s, vk , . . . , vm2 , wm1 , . . . , wℓ−1 },

V2 = V¯1 .

It is trivial from the definition of j that the edges e = (wi , t), i =
{m1 , . . . , ℓ − 1} are saturated.
Since x1 , x2 ≥ 0, and k is the minimal value such that c(vk , wℓ ) = 0,
we have c(vi , wℓ ) = 1 for i = m1 , . . . , k − 1. Since (wℓ , t) is not saturated,
all edges (s, vi ), i ∈ {m1 , . . . , k − 1} must be saturated.
From the definition of k and the non-negativity of the portfolio vector it is easy to see that edges e = (vi , wj ) for i ∈ {k, . . . , m2 },
j ∈ {ℓ, . . . , m2 } and positive capacity cannot exist. Thus, every edge
e = (x, y) with x ∈ V1 and y ∈ V2 is saturated. The Max-Flow-Min-Cut
Theorem then implies that the algorithm indeed computes a maximal
flow.
Observing the fact that in each loop iteration either index i is decremented or index j is incremented, and that there are only m different
values that either i or j can take on before the algorithm terminates,
there are at most 2m − 1 loop iterations, and the linear running time
bound follows.
To compute inf{RAw (α, ~x)| xi = 1} we have to compute RAw (α, ~x)
for all possible portfolios hx1 , x2 i. However, each feasible portfolio corresponds to a half-space (as in Figure 1.1) defined by a line that goes
through the point (α, α) (x1 α + x2 α = α, since x1 + x2 = 1), so we only
need to consider the O(m2 ) distinct subsets of return pairs that can be
defined by a line going through (α, α). We can identify each such portfolio with a different (non-positive) slope s1 , . . . , sm2 , which we assume
to be sorted in descending order. By using a suitable data structure it
is possible to compute the best portfolio much faster than the obvious
O(m3 ) algorithm that starts the greedy algorithm for each slope.
P

Theorem 3 Given S1 , S2 , and α, we can compute in O(m2 log m) arithmetic operations a portfolio hx1 , x2 i for a worst-case risk-averse investor
which minimizes equation (1.2).
Proof: Starting with the first slope s1 we build up a binary tree.
Each is labeled with a pair of two real entries (e1 , e2 ). The leaves of the
tree correspond to the rows and the columns in the following way.
Starting from column m2 we add leaves from left to right. We add
leaves with labels (0, S2 (m1 µ)), (0, S2 ((m1 + 1)µ)), . . ., (0, S2 (jm µ)), until we reach a row index jm such that x1 m2 µ + x2 (jm + 1)µ > α,

10
i.e., this index is the last under the crucial line. To be precise we
let jm = ⌊ α−xx21µm2 µ ⌋; note that it may be the case that jm < m1 ,
so this sequence of leaves may be empty. Then we add the leaf
(−S1 (m2 µ), 0). Next, we consider column m2 − 1 and add leaves
(0, S2 ((jm + 1)µ), . . . , (0, S2 ((jm−1 )µ)), until we reach an index jm−1 ,
such that x1 (m2 − 1)µ + x2 (jm−1 + 1)µ > α. Then we add the leaf
(−S1 ((m2 − 1)µ), 0) and proceed similarly with column m2 − 2. Note
that the order of adding leaves is crucial to this data structure and the
correctness of the algorithm is based on that. Starting from left to right
we group the leaves in pairs of 2 and build a parent node for each pair
according to the following rule
parent[(e1 , e2 ), (f1 , f2 )] = (e1 + min{e2 + f1 , 0}, max{e2 + f1 , 0} + f2 ).
We build O(log m) layers iteratively, until we reach a single root node
(r1 , r2 ). It is easy to see that this tree based algorithm imitates the
greedy algorithm described before and that 1 + r1 = 1 − r2 is exactly the
flow value. Building this tree structure takes constant time per tree node,
and since there are O(m) nodes we have a total time of O(m), which is
no better than the time bound of the greedy algorithm. The advantage
is that we can dynamically update this data structure efficiently.
We will first sort all of the m2 possible return pairs by their slope
with the point (α, α), so that as the slope determined by our portfolio
increases we can quickly (in constant time per pair) determine which
pairs are added and which are removed from our half-space of interest.
This takes O(m2 log m) time. To update our data structure for each
point insertion/removal, all that is required is swapping the position of
two neighboring leaves. With obvious techniques, the positions of these
two leaves can be found in O(1) time, and we can update the tree by
looking at the path from the two leaves to the root and update each
node on that path. Each update step requires O(1) operations and the
length of the path is bounded by O(log m). Since there are at most
m2 point additions and removals, each taking O(log m) time, it takes at
most O(m2 log m) time to consider all possible portfolios.

3.2.

The Average-Case Investor

For the average-case investor (RAa or AGa ), we are not interested
in the extremes of the joint distributions, but rather the distribution of
the feasible tables. In this section we consider Q = Lα,~x (M ) a random
variable where M is drawn from a uniform distribution over the feasible
tables Mk . The definition of RAa (α, ~x), from (1.3), is then E[Q]. We
will see that computing the distribution function of Q is a computation-

11

The Risk Profile Problem forStock Portfolio Optimization

ally difficult problem to solve exactly, but can be approximated within
a reasonable (polynomial) amount of time.
Theorem 4 Let γ ∈ [0, 1] be an n-bit rational. It is ♯P -hard to compute
the fraction of feasible tables M ∈ M2 with
Lα,~x (M ) =

X

Mδ ≤ γ

δ∈L(α,~
x)

(the integration of the corresponding indicator function, or the distribution function for Q).
Proof: Given positive integers a1 , . . . , an , b, it is shown in [Dyer and
Frieze, 1991] that computing the n-dimensional volume of the polyhedron P
n
X

aj y j ≤ b

0 ≤ yj ≤ 1 (j = 1, . . . , n)

Pn

and consider the polyhedron

j=1

is ♯P -hard. Let d =
n+1
X

j=1 aj

aj y j = d

0 ≤ yj ≤ 1 (j = 1, . . . , n + 1),

(1.19)

j=1

where an+1 = d. Note P
that for any valid assignment of values to
y1 , y2 , . . . , yn we have 0 ≤ nj=1 aj yj ≤ d, so there is a yn+1 ∈ [0, 1] that
will satisfy (1.19). Now let a′i = ai /(2d) and define a 2 × (n + 1) contingency table by t1j = a′j yj , t2j = a′j (1 − yj ), with row sums (1/2, 1/2)
and column sums (a′1 , . . . , a′n+1 ).
To completely define our stock problem, we must also give values for
µ, α, the portfolio ~x = hx1 , x2 i, and the threshold γ, which we do as
follows:
µ = 1,

x1 =

1
,
n+1

x2 =

n
,
n+1

α=

2n
,
n+1

γ=

b
.
2d

It is straightforward to verify from these values that the return pairs
in the critical region (the shaded region in Figure 1.1) are exactly the
entries t1j for j = 1, . . . , n. Therefore, the tables that satisfy our criteria,
that Lα,~x (M ) ≤ γ, are precisely those with
n
X

j=1

t1j ≤ γ

⇐⇒

n
X

j=1

a′j yj ≤ γ

⇐⇒

n
X

j=1

aj yj ≤ γ · 2d = b.

12
Therefore the feasible tables that meet our criteria are exactly those that
correspond to points in polyhedron P , and so the fraction of tables that
meet the criteria is exactly the volume of P .
Following the notation of Dyer, Kannan and Mount [Dyer et al., 1997],
who describe a sampling procedure for contingency tables with integer
entries and large row and column sums (≥ Ω(m3 )), we define
V (r, c) =




x ∈ Rm×m |



X

xij = ri for i = 1, . . . , m

j

and

X

xij = cj for j = 1, . . . , m

i

)

and
P (r, c) = V (r, c) ∩ {x|xij ≥ 0 for i = 1, . . . , m, j = 1, . . . , m}
as the contingency polytope. Thus, V (r, c) is the set of matrices with
row and column sums specified by r and c respectively. In our case
ri = S1 (iµ), ci = S2 (iµ) and P (r, c) is the set of joint distributions Mk .
Let U be the lattice
{x ∈ Zm×m |

X

xij = 0 for i = 1, . . . , m,

X

xij = 0 for j = 1, . . . , m}.

i

j

For 1 ≤ i ≤ m − 1 and 1 ≤ j ≤ m − 1, let b(ij) be the vector in Rm×m
given by b(ij)i,j = 1, b(ij)i+1,j = −1, b(ij)i,j+1 = −1, b(ij)i+1,j+1 = 1
and b(ij)k,ℓ = 0 for all other indices k, ℓ. Any vector x in V (0, 0) can be
expressed as linear combination of the b(ij)’s as follows
x=

m−1
X
X m−1
k=1 ℓ=1



k X
ℓ
X

xij  b(kℓ).
i=1 j=1

It is easy to see that the b(ij) are all linearly independent and the the
dimension of V (r, c) and P (r, c) for positive row and column sum vectors
r and c is (m − 1)2 [Dyer et al., 1997]. We will apply the sampling
algorithm pioneered by Dyer, Frieze and Kannan [Dyer et al., 1991] and
later refined in a sequence of papers (see [Kannan, 1994] for an overview)
to sample uniformly at random in P (r, c).
We sample in the space V (r, c). As mentioned in the introduction, we
know a starting point z0 in P (r, c) (multiplication of rows and column
sums). It is easy to see that a ball of radius b2 is inside P (r, c), if every
component of r and c is at least b. Since in our case r and c sum up

The Risk Profile Problem forStock Portfolio Optimization

13

to one, P (r, c) ⊂ B(0, 1). The following theorem is a corollary of the
analysis of the fastest sampling algorithm in convex bodies known so far
by Kannan, Lovász and Simonovits [Kannan et al., 1997].
Theorem 5 We can generate a point in P (r, s), which is almost uniform in the sense that its distribution is at most ǫ away from the uni6
form in total variation distance. The algorithm uses O ∗ ( mb4 ) membership
queries of P (r, s) (each requires O(m2 ) arithmetic operations).
procedure Estimate(x)
S←0
N = 100
ǫ2 δ
for ℓ = 1, . . . , N do
ζi ← result from sample procedure started at x
S ← S + Lα,~x (ζi )
end for
S ← S/N
return S
Figure 1.3.

The approximation algorithm

Theorem
 8 6 Procedure Estimate (in Figure 1.3) computes a number S
∗
in O b4mǫ2 δ arithmetic operations, which approximates RAa (α, ~x) (i.e.,
RAa (α, ~x) − ǫ ≤ S ≤ RAa (α, ~x) + ǫ) with probability 1 − δ.
Thus, E(Sk ) =
Proof:
Let Sk = k1 ki=1 Lα,~x (ζi ).
R
Lα,~x (M )w(M )dM , where w is the density produced by the random
walk. Since 0 ≤ Lα,~x (M ) ≤ 1 for all M ∈ M2 , it is easy to see that
σ 2 (S1 ) ≤ 1 and so σ 2 (Sk ) ≤ k1 . By Chebychev’s inequality,
P

P (|Sk − E(Sk )| ≥ ǫ/2) ≤

4
σ 2 (Sk )
≤ 2 .
2
(ǫ/2)
ǫ k

Since the samples are not entirely uniform, we must consider the error
introduced by the approximately uniform sampling distribution as well.
Let uMk (M ) denote a uniform density over the set Mk , and then approximating a uniform distribution within bound ǫ/4, Theorem 5 implies
|E(Sk ) − RAa (α, ~x)|
=

Z

Lα,~x (M )w(M )dM −

Z

Lα,~x (M )uMk (M )dM

14
≤

Z

w>uMk

+

Z

(w(M ) − uMk (M )) dM

w≤uMk

(uMk (M ) − w(M )) dM

≤ ǫ/2.
Setting k =

4.

4
ǫ2 δ

the theorem follows.

The k-Stock Case

In this chapter we consider the general case of more than two stocks.
Since the problem of estimating the probability distribution for the
average-case investor is already ♯-P complete in the two stock case, we
do not consider it any more and concentrate on a worst-case investor.
We start with a complexity result for three stocks, which implies that a
greedy or flow based portfolio is quite unlikely to exist.
Theorem 7 The existence of a greedy or flow based portfolio for the
problem with 3 or more stocks implies P = N P .
Proof:
We prove this result by reduction from NUMERICAL3-DIM-MATCHING. Consider an instance of NUMERICAL-3-DIMMATCHING, i.e., disjoint sets X1 , X2 , X3 , each containing m elements,
a size s(a) ∈ Z+ for each element a ∈ X1 ∪ X2 ∪ X3 and bound B ∈ Z.
We would like to know if X1 ∪ X2 ∪ X3 can be partitioned into m disjoint
sets such that each of these sets contains exactly one element from each
of X1 , X2 , and X3 , and the sum of the elements is exactly B (we can
change this requirement to ≤ B without difficulty). This problem is
NP-complete in the strong sense, so we restrict the sizes to be bounded
by a polynomial, s(a) ≤ nc for some constant c.
We construct an instance of the problem of computing
RAw (α, h1/3, 1/3, 1/3i) by making a contingency table in which
Sk (i) = ck,i /m, where ck,i is the number of items in set Xk with value i.
The existence of a greedy or flow based algorithm implies the existence
of a solution in which all entries in the solution table are multiples of
1/m, and such a solution exists with Lα,~x (M ) = 1 if and only if there
is a valid partition of X1 ∪ X2 ∪ X3 . If such a partition exists, we can
find it by simply taking all of the triples “selected” (with multiplicity
determined by the integer multiple of 1/m), and use elements from X1 ,
X2 , and X3 as determined by the three coordinates of each selected
point.
While this proof shows that it is unlikely that a fast and simple greedy
or flow-based algorithm exists, as it does for 2 stocks, we can indeed solve

The Risk Profile Problem forStock Portfolio Optimization

15

the problem for a fixed number of stocks in polynomial time using a more
time-consuming procedure based on linear programming. This is stated
in a general setting in the following theorem.
Theorem 8 If the number of stocks k is part of the input, the problem
of determining the best portfolio for a worst-case investor can be solved
in time polynomial in the number of entries of the contingency table (but
exponential in k).
Proof:
The problem can be modeled as linear program with a
number of variables, that corresponds to the number of entries of the
contingency table, and km inequalities.

4.1.

An Approximation Algorithm

In this section we describe an approximation algorithm, that solves
the problem of determining the worst case probability for a given portfolio within a given error ǫ ∈ R+ in polynomial time. Additionally, we
describe an important, non-trivial special case, where the problem can
be solved exactly in polynomial time.
Theorem 9 Suppose that a portfolio hxi iki=1 and a target return α are
given. The worst-case probability can be approximated (i.e., we compute
a value W with RAw (α, ~x)−ǫ ≤ W ≤ RAw (α, ~x)+ǫ) in time polynomial
in k and n. The number of steps is dominated by solving a linear program
in O(km2 /ǫ2 ) variables and O(km/ǫ) constraints.
Proof: We consider the first pair of stocks S1 and S2 as in the
1
and
two dimensional case and define a new portfolio as x˜1 = x1x+x
2
x2
1
x˜2 = x1 +x2 . We divide the two dimensional plane in ℓ = ǫ m log k
regions by ℓ parallel lines x˜1 x + x˜2 y = const of constant distance. Thus,
we divide the entries of the joint distribution matrix into ℓ different sets
(see Figure 1.4).
Each entry in the matrix corresponds to a variable and the variables
satisfy the row sum and column sum condition of the joint distribution.
Next, we sum up the entries in the ℓ different sets and assign the sums
to ℓ new variables. By combining these sum variables from two different pairs of stocks, we get a new table with new row and column sum
conditions, resulting again in ℓ new sum variables.
Repeating combinations in this manner, we stop after log k iterations
and the creation of O(km2 log k/ǫ2 ) variables and O(km log k/ǫ) constraints, leaving just one table with 2 border distributions (expressed as
variables). Assuming, that the variables of the border distributions cor-

16

Figure 1.4.

Striping idea used in worst-case approximation construction

respond to the distribution of the stocks S1 , . . . , Sk/2 and Sk/2+1 , . . . , Sk ,
we do the following.
x
x +···+x
P+···+xn for our
We define a portfolio x˜1 = 1 P x k/2 and x˜2 = k/2+1
x
i

i

last table and consider the line x˜1 x + x˜2 y = α, dividing our last table in
two sets. The variables below that line are summed up and we solve a
linear program by maximizing this sum subject to the constraints created
before. Since we reduced the number of entries in each table from Ω(m2 )
to only ℓ, that are considered in the next table, we lost some precision
during the combination. But, after the first pairing in the lowest level
of the binary tree, each sum variable represents a loss probability of the
combination of the two stocks within an error of logǫ k %. Furthermore,
it is easy to see that during the repeated combination of the stocks the
error accumulates linearly in each iteration. Thus, the theorem follows.

Theorem 10 Suppose that a portfolio hxi iki=1 and a target return probability p is given. Under the assumption, that the dollar, that has to be
invested, can only be broken into a fixed number c of equal units (cents),
the worst-case probability can be computed exactly in time polynomial in
k and m.
Proof: The proof is based on a similar construction as the approximation algorithm and is omitted for brevity.

REFERENCES

17

Acknowledgments
The authors wish to thank the anonymous referees for very helpful
comments.
A preliminary version of this work appeared in Proceedings of the 32nd
Annual ACM Symposium on Theory of Computing, pages 228–234, 2000.

References
[Cormen et al., 1990] Cormen, T. H., Leiserson, C. L., and Rivest, R. L.
(1990). Introduction to Algorithms. MIT Press, Cambridge, MA.
[Duffie, 1996] Duffie, D. (1996). Dynamic Asset Pricing Theory. Princeton University Press, Princeton, NJ, 2nd edition.
[Dyer and Frieze, 1991] Dyer, M. and Frieze, A. (1991). Computing the
volume of convex bodies: a case where randomness provably helps. In
Probabilistic combinatorics and its applications, pages 123–169. American Mathematical Society, Providence, RI.
[Dyer et al., 1991] Dyer, M., Frieze, A., and Kannan, R. (1991). A random polynomial-time algorithm for approximating the volume of convex bodies. Journal of the ACM, 38(1):1–17.
[Dyer et al., 1997] Dyer, M., Kannan, R., and Mount, J. (1997).
Sampling contingency tables. Random Structures & Algorithms,
10(4):487–506.
[Elliott and Kopp, 1999] Elliott, R. J. and Kopp, P. E. (1999). Mathematics of financial markets. Springer-Verlag, New York, NY.
[Fouque et al., 2000] Fouque, J. P., Papanicolaou, G., and Ronnie, S. K.
(2000). Derivatives in Financial Markets with Stochastic Volatility.
Cambridge University Press, London.
[Goldberg and Rao, 1998] Goldberg, A. V. and Rao, S. (1998). Beyond
the flow decomposition barrier. Journal of the ACM, 45(5):783–797.
[Hull, 2000] Hull, J. C. (2000). Options, Futures, and Other Derivatives.
Prentice Hall, Upper Saddle River, NJ, 4 edition.
[Kannan, 1994] Kannan, R. (1994). Markov chains and polynomial time
algorithms. In Proceedings of the 35th Annual IEEE Symposium on
Foundations of Computer Science, pages 656–671.
[Kannan et al., 1997] Kannan, R., Lovász, L., and Simonovits, M.
(1997). Random walks and an O ∗ (n5 ) volume algorithm for convex
bodies. Random Structures & Algorithms, 11(1):1–50.
[Karatzas, 1997] Karatzas, I. (1997). Lectures on the mathematics of
finance. American Mathematical Society, Providence, RI.

18
[Karatzas and Shreve, 1998] Karatzas, I. and Shreve, S. E. (1998).
Methods of Mathematical Finance, volume 39 of Applications of Mathematics. Springer-Verlag, New York, NY.
[Musiela and Rutkowski, 1997] Musiela, M. and Rutkowski, M. (1997).
Martingale methods in financial modelling. Springer-Verlag, Berlin.
[Sharpe et al., 1995] Sharpe, W. F., Alexander, G. J., and Bailey, J. V.
(1995). Investments. Prentice Hall, Upper Saddle River, NJ, 5th
edition.

