Annals of Operations Research 99, 189–206, 2000 2001 Kluwer Academic Publishers. Manufactured in The Netherlands.
The AURORA Financial Management System: Model and Parallel Implementation Design ∗ ´ ETANOWSKI G.CH. PFLUG and A. SWI ¸
[email protected],
[email protected] Department of Statistics and Decision Support Systems, Vienna University, Universitätsstr. 5, A-1010 Wien, Austria E. DOCKNER and H. MORITSCH
[email protected],
[email protected] Department of Management, Vienna University, Bruennerstr. 72, A-1210 Wien, Austria
Abstract. The AURORA financial management system under development at the University of Vienna is a modular decision support tool for portfolio and asset–liability management. It is based on a multivariate Markovian birth-and-death factor model for the economic environment, a pricing model for the financial instruments and an objective function which is flexible enough to express risk aversion. The core of the system is a large scale linear or convex program, which due to its size and structure is well suited for parallel optimization methods. As the system is still at an early stage of development, the results are preliminary in nature. Only a few types of financial instruments are handled and just two types of objectives are considered. The parallel optimization modules are still in the development phase. Keywords: financial management, asset–liability management, stochastic optimization, parallel algorithms AMS subject classification: primary 90C15, secondary 90B50
1.
Introduction
The classical model in financial management systems has been developed by Markowitz [19] more than forty years ago. He introduced a quadratic programming model for strategic asset allocation decisions of investors who plan their portfolio composition over asingle period. In this framework an investor chooses a portfolio of different assets, or different asset classes, so that the variance (the overall risk) of the portfolio is minimized subject to the constraints that the expected portfolio return should exceed a prespecified level. Since this quadratic programming problem is parameterized in the expected return it can be solved to yield the efficient frontier of optimal risk and return specifications for the assets or the asset classes. While the Markowitz model has become an industry standard, it suffers at least from two important shortcomings both of which relate to the lack of a dynamic structure. The asset allocation decision that is derived through the model is essentially a single ∗ This research is conducted within sub-project 6 of the Special Research Program SFB F011 AURORA
supported by the Austrian Science Fund.
190
PFLUG ET AL.
period myopic strategy. It therefore assumes that the investor has a one period horizon with a mean–variance objective function. As Brennan et al. [7] point out this is inconsistent with the data used in most practical applications of this model. The expected returns that are an input to the Markowitz model are typically not one period (short run) returns but rather long run expected returns or internal rates of returns of bonds with a maturity that goes way beyond the next period. If such long run estimates are used, the resulting efficient portfolio must be biased in one way or the other and does not yield the results aimed for by the decision maker. A second drawback of the traditional Markowitz framework pointed out by Brennan et al. relates to the choice of the objective function. With a one period variance criterion investors choose their asset allocation decisions in such a way that they minimize current risk and ignore hedging possibilities against shifts in future risks or future investment opportunities. The above argumentation suggests that multiperiod dynamic stochastic models are more suitable to realistically reflect the financial planning problem associated with tactical asset decisions or longer term risk management. The literature on dynamic financial planning problems can be classified into two different approaches. One in which the dynamic model is formulated as a continuous time stochastic optimal control model (Brennan et al.) and one in which a discrete time formulation is used that results in a multi-period stochastic optimization problem in the spirit of [8,20–22,26]. See also the books [39] and [25]. The difference between the two approaches rests on the decision structure and on the modeling of the stochastic process which is the input to the financial planning problem. The optimal control approach assumes that the returns are generated by different risk factors that are represented by stochastic differential equations. The (time continuous) decisions are modeled as functions which depend only on the actual state of the system expressed by the differential equations. Multistage dynamic planning models are based on the assumption that the risk factors are generated by a discrete stochastic process, and the decision variables may depend on the whole history of the process as well as on exogeneous parameters, like in- or outflow of capital. Brennan et al. apply the optimal control approach and demonstrate how such a model can be solved using stochastic dynamic programming techniques. In particular, they derive a set of necessary and sufficient optimality conditions that result in a system of partial differential equations in the state variables (the risk factors of the model), which they solve by applying numerical techniques. While they claim that this approach has several advantages over a multiperiod stochastic programming model, we argue in the opposite direction. Although a stochastic optimal control approach allows for a flexible specification of the state variables of the financial model, the partial differential equations need to be discretized as well as solved numerically which can be a complicated task and the solution procedures can be sensitive to the choice of the nonlinear objective function. Secondly, if the model is calibrated to actual data, a continuous time model is not consistent with the way observations are available to estimate the parameters of the modeled stochastic process. On the contrary, a multistage stochastic optimization
THE AURORA FINANCIAL MANAGEMENT SYSTEM
191
framework allows for a modeling of the economic environment as (for instance) a multifactor Markov model that can easily be estimated using actual market data and for which scenarios can be generated that characterize the different possible evolutions of the financial system. Moreover the dynamic stochastic programming model allows to include transaction costs, in- and outflow of capital and the dependence of the decision on the whole past. Multistage stochastic optimization is being increasingly used within the fields of portfolio planning, asset & liability management, risk management and asset allocation problems. Models of this type have been proposed by [8,12,18,24,39]. In this paper we present the AURORA financial management system, that is currently being developed at the University of Vienna. Its main features are: 1. A new technique to model and estimate the financial and economic environment within the framework of a multi-factor Markov model of the birth-and-death type. We go much beyond the classical binary lattice model, such as the Black et al. [6] model. Much emphasis is put on the statistical estimation: we apply principal components analysis to extract common factors that drive the dynamics of our financial variables such as asset returns and interest rates. We allow different specifications of factors which increases the flexibility of the system. 2. The specification of the objective function. Here we deviate from the traditional approach of expected terminal wealth maximization and introduce a parameterized risk functional. This approach integrates several advantages. Firstly, the parametrization provides us with the possibility of evaluating the sensitivity of the optimal solution with respect to the choice of the parameters (as it does in the case of, e.g., parametric utility functions). Secondly, changing the objective function itself reveals possible sensitivity to the form of the functional. This is important since it allows us to shed some light on the role of an investor’s objective when designing an optimal portfolio or risk management strategy. Clearly, all risk adjusted objective formulations allow us to drop the often unrealistic assumption of the decision maker’s risk neutrality. 3. A decomposed and parallel implementation. This part of the project is still in the design phase. We plan to exploit the special structure of the linear or convex optimization problem resulting from our model. An outline of one possible design of a parallel implementation is presented. Our paper is organized as follows. In the next section we present the structure of our financial model and point out its flexibility with respect to the specification of the financial contracts, the objective function, as well as the calibration of the data to the economic environment. In section 3 we present a short description of the features and estimation procedures for our multi-factor scenario model. Since our scenario approach rests on a discrete time Markov model we deal in section 4 with the approximation of a general diffusion process by a birth-and-death Markov chain. In section 5 we introduce our data parallel solution method’s design principles. Section 6 concludes the paper.
192
2.
PFLUG ET AL.
Modeling the system structure
In this section we present the building blocks of our financial planning model and how they are integrated to make up a complete system. We distinguish four different components of our model: 1. The financial instruments that are available to the investors. The set of financial instruments determine the feasible set of choices for an investor. In particular they determine in which way an investor can set up his portfolio or design a a dynamic hedging or an asset and liability strategy, and from which set of asset classes he can choose when designing asset allocation decisions. 2. The economic and financial environment in which the application (i.e., the asset allocation problem, the risk management or the portfolio optimization model) is suited. Specifying the economic environment means that we identify those variables that drive the dynamics of the prices and returns of the financial instruments. These variables are assumed to be exogeneous and in most of the applications they are the sources of risk that an investor has to deal with. Moreover, it is necessary to specify how the risk factors are mapped onto the prices and returns of the financial instruments. Hence, the financial planning model must include a valuation procedure for each instrument. 3. The objective function of the investor that reflects his understanding of risk as well as his risk preferences (see [31] for a glimpse at a case study done within this project). A careful modeling of the objective function is very important in financial optimization problems since many portfolio strategies are very sensitive with respect to the choice of the objective function. 4. The model equations that identify the constraints an investor faces as well as the balance equations for holding, buying or selling assets. The dynamic structure of the portfolio determines his current and future level of wealth. Figure 1 gives a graphical representation of our modeling system. In particular it identifies the four building blocks and how they comprise the financial planning model, which at the and is formulated as an large scale linear or nonlinear programming problem. In what follows we will discuss the four building blocks in detail and we give an analytical representation of them. 2.1. Modeling the financial instruments The AURORA financial planning problem allows modeling of a variety of different financial instruments such as bonds, stocks, forward contracts, futures, options, credits, loans and other contracts which may be bought and sold at some of the decision moments in a discrete time environment.
THE AURORA FINANCIAL MANAGEMENT SYSTEM
193
Figure 1. The structure of the AURORA financial management system: defining an optimization problem.
Each contract has a contract number j ∈ J = {1, . . . , J } and is characterized by four functions fb , fs , fv , fc , which in the most general case depend on the history of the economic environment X (and hence on time) and the contract number j . The four functions represent: • The buy-price fb (X , j ) of one unit of contract j .
194
PFLUG ET AL.
• The sell-price fs (X , j ) of one unit of contract j . • The value of the contract fv (X , j ). This value function is used to evaluate a contract at periods of time when it is not traded, e.g., the market value of a loan during the maturity of the contract. • The cash-flow fc (X , j ) generated by one unit of contract j , like dividends, coupon payments, interest payments etc. To determine the prices as well as the values of the contracts we apply valuation tools that determine the functions f (X , j ). Any contract that may be adequately described by these four functions and the time of its availability for trading, may in principle be handled within our model. This includes, but is not limited to, all the contract categories listed above. It should be stressed, however, that the current model does not include fixed transaction costs. In an environment in which such costs are an important part of the overall transaction costs and when they cannot be adequately modeled by proportional costs (i.e., the bid–ask spread), our model cannot be applied. At the moment we adopt a simple rule that determines the sell and the buy price of a contract. We use security pricing theory to determine the equilibrium price of an asset which is identical to its value at a given point in time. We assume that this equilibrium price is the average between the bid and the ask prices. While the bid–ask spread may arguably need modeling as a separate (correlated) stochastic process, in our current application study on asset alocation for pension funds (see [31]) it is treated as a constant. The modeling of bid–ask spread may be included in the model whenever a need should arise. 2.2. The economic environment The economic environment is characterized by a set of exogeneous variables that drive the dynamics of all prices and returns of the financial instruments. Examples of such variables include • spot interest rates, • exchange rates, • stock market indices. These financial variables are used as factors in multi-factor pricing models to determine the buy and sell prices of contracts, as well as the contracts’ value and generated cash flows. Clearly, the random variables listed above are known to be correlated. Having identified a minimal set of such variables needed to price all selected contracts, we proceed with a principal component analysis. The identified components are then used as the underlying “hidden” risk factors. It is more than likely that the principal components will not have any clear economic interpretation. After extracting their statistical properties from historic data, we are able to generate a scenario tree of their development (see sections 3 and 4).
THE AURORA FINANCIAL MANAGEMENT SYSTEM
195
The dynamics of risk factors is modeled by means of a discrete-time, discrete-state stochastic vector process X(t) = (X1 (t), . . . , Xd (t)); t = 1, 2, . . . , T . With the process X(t) we associate the history process Xt = (X(1), X(2), . . . , X(t)). Since the process X has only finitely many states and the starting state X(1) is deterministic, one may arrange all states of the history process in a finite tree, called the scenario tree. Its root is X(1). The nodes of the tree may be numbered 1 n N such that each node number n corresponds in a one-to-one way to a history path of the process X(·). To each node n the probability p(n) of reaching this node is associated. We consider the history as a function of the node number of the scenario tree, writing X (n) instead of Xt . 2.3. Modeling the objective function We assume that the objective of the investors is related to the level of his terminal wealth. Let W be the terminal wealth. This is a random variable with distribution depending on the decision variables. Then we specify the objective of the investor as max E(W ) − ρE W − E(W ) , where ρ measures the degree of risk aversion. The mean absolute deviation is a good measure of the risk associated with the decision and the risk aversion factor allows to adapt the objective to the specific views of the decision maker (see [1,17,28,30]. Moreover, it preserves linearity: if the terminal wealth W is linear in the decision variables x, then the objective is a piecewise linear function in x. As a variant, we also consider the lower-semivariance as a risk criterion, i.e., − 2 , max E(W ) − ρ E W − E(W ) where [·]− denotes the negative part. This objective is convex in the decisions x, if W is linear in x, as it was shown by Pflug in [30]. 2.4. Modeling the constraints The structural equations of our model consist of the accumulation equation for cash holdings, the evolution of the investor’s wealth and some bookkeeping identities that specify the number of contracts held in the current portfolio. The decision variables are the number of contracts bought xb (n, j ) or sold xs (n, j ) of contract type j at node n of the scenario tree. Let y(n, j ) the amount of contract j held right after the decisions are made at node n. The bookkeeping equations are ∀(n ∈ N , j ∈ J )
y(n, j ) = y(n− , j ) + xb (n, j ) − xs (n, j ),
where n− denotes the node number of the predecessor of n.
196
PFLUG ET AL.
Let c(n) denote the cash at node n. Cash holdings accumulate according to ∀(n ∈ N )
−
c(n) = c(n ) +
J
xs (n, j )fs X (n), j − xb (n, j )fb X (n), j .
j =1
Here no interest on cash is assumed. In cases where the cash (or part of it) is actually placed in a bank account, interest may be added to the above equation. Alternatively, assets held in banks may be treated as yet another category of contracts. Terminal wealth can then be calculated as follows. Let M be the set of terminal nodes in the scenario tree. Terminal wealth consists of cash held at the terminal node as well as the number of contracts of security j times the market value of one unit: ∀(n ∈ M)
W (n) = c(n) +
J
y(n, j ) fv X (n), j .
j =1
The decision variables are nonnegative and additional constraints may be added. The objective function (see section 2.3) is max p(m)W (m) W (n)p(n) − ρ W (n) − n∈M
m∈M
which can be written equivalently in the form of a linear objective max
p(n) W (n) − ρ V + (n) + V − (n)
n∈M
∀(n ∈ M)
(1)
W (n) − m∈M p(m)W (m) = V + (n) − V − (n), V + (n), V − (n) 0.
A convex nonlinear variant is presented in the conclusions (see (20)). 3.
Modeling and estimation of the scenarios
The estimation of the scenario tree is done in several steps: • the relevant economic time series (interest rates, exchange rates, stock indices) are identified and historical data of these time series are collected, • a factor model is extracted from historical data, • the factors are discretized and a multivariate birth-and-death type Markov chain is fitted to the data, • a reconstruction formula maps the states of the Markov chain to the values of factors, • the simulated factors are expanded back into the meaningful economic time series; in this step the expert opinion of the modeler may be built in.
THE AURORA FINANCIAL MANAGEMENT SYSTEM
197
Let Y(τ ) = (Y1 (τ ), . . . , Ym (τ )) be a multidimensional historical time series of length T . The modeler may or may not remove level and trend from this series, i.e., to consider Y(τ ) = Y(τ ) − a0 − b0 τ , in place of Y(τ ). This decision is dependent on the nature of the observed time series as well as on the decision maker’s opinion. The discussion of all possible cases is not possible within the scope of this article. Let C=
T T Y(τ ). Y(τ ) τ =1
Since the matrix C is symmetric and nonnegative definite, we may represent it as C = U DU T, where U = (u1 , . . . , um ) is the matrix of eigenvectors ui of C and D = diag(λ1 , . . . , λm ) is the diagonal matrix of (nonnegative) eigenvalues. We assume that λ1 λ2 · · · λm . Let d m and B be the [m × d] submatrix of U , which consists of the first d columns of U B = (u1 , . . . , ud ). The reduced times series is Z(τ ) = (Z1 (τ ), . . . , Zd (τ ))T = B T · Y(τ ). It is d-dimensional and the best approximation (in quadratic sense) of Y(τ ) by a d-dimensional series is
Y(τ ) = B · Z(τ ). The process Z(τ ) is discretized and transformed into a coarser time scale. The historic observations Y(τ ) are typically on a daily basis, whereas the time scale of the stochastic optimization model is months or years. Let the scale factor be κ. By discretizing the range of each component of Z(τ ) into classes and estimating the transition behavior of Z(τ ) to Z(τ +κ) we get a model X(t) = (X1 (t), . . . , Xd (t)). This discretization is described in the next section. Each component Xi takes only a finite number of values and is assumed to be independent of the other components. We assume that Xi is a birth-and-death Markov chain because these processes – as is shown in the next section – are flexible enough to allow approximations to general diffusion processes. The starting values X(1) are set in such a way that they are consistent with with the actual observation. Finally, the modeler may again have a guess for the trend for the future time series. Naturally, since X are mathematical artefacts, a meaningful intervention by the decision maker may only be made at the point at which the values of Y are recovered from X. 4.
Approximating a diffusion process by a birth-and-death Markov chain
There is a considerable gap between the modern theory of financial processes and the usual simple lattice structures proposed in the stochastic optimization literature. In this
198
PFLUG ET AL.
section, we show how to close this gap by defining birth-and-death (B&D) Markov chains, i.e., chains on the integer grid, which allow transitions only between neighbors. As we will see, B&D chains – as simple as they are – are universal in the sense, that they approximate the class of all diffusion processes which posses a stationary distribution. Conditions for the existence of a stationary distribution can be found in [16]. Suppose that a financial process is described by the stochastic differential equation dX(t) = f X(t) dt + σ X(t) dW (t) (2) such that X(t) is a recurrent Markov process. Examples include the Vasicek model [37] dX(t) = a µ − X(t) dt + σ dW (t),
(3)
the Cox–Ingersoll–Ross model [9] dX(t) = a µ − X(t) dt + σ X(t) dW (t)
(4)
and the exponential Ornstein–Uhlenbeck process 1 2 dX(t) = −aX(t) log X(t) + X(t)σ dt + σ X(t) dW (t). 2
(5)
The mean reverting exponential Brownian motion dX(t) = a µ − X(t) dt + σ X(t) dW (t), which is used extensively in theory (e.g., for the famous Black–Scholes formula), does not belong to this class, since it describes a transient process possessing no stationary distribution, which seems unnatural for a large class of financial processes. We will approximate (2) by a birth-and-death process. Recall that a birth-and-death process allows only three movements from every state: an upward jump, a downward jump and no jump (see figure 2). Let H be a strictly monotone, three times differentiable real function and define 1 h (H −1 (u))σ 2 (u) f (u) − h(H −1 (u)) 2 h3 (H −1 (u)) σ (u) . τ (u) = h(H −1(u))
g(u) =
and
(6) (7)
The function H must be chosen in such a way that |τ (x)| 1 for all x. Inverting the equations (6) and (7) we get 1 f (x) = h H −1 (x) g(x) + h H −1 (x) τ 2 (x) 2 σ (x) = h H −1 (x) τ (x).
and
(8) (9)
THE AURORA FINANCIAL MANAGEMENT SYSTEM
199
Figure 2. A birth-and-death Markov chain.
Define for a large integer N and any integer i the following probabilities: the probability of an upward jump 1 1 i 1 2 i (N) τ H + g H , (10) pi = 2 N N N 0 the probability of a downward jump 1 1 i 1 2 i τ H − g H qi(N) = 2 N N N 0
(11)
and the probability of no jump ri(N) = 1 − pi(N) − qi(N) . Here [a]10 = min(max(a, 0), 1). The approximating Markov chain Xn(N) is now defined on the states H (i/N) and has the transition probabilities pi(N) , qi(N) , ri(N) as just described. We have for large enough N i (N) − Xn(N) | Xn(N) = H E Xn+1 N i i i−1 i+1 (N) −H pi + H −H qi(N) = H N N N N i 1 i 1 1 i i 2 g + τ H +o = 2h h N N N 2N 2 N N N2 1 i 1 +o = 2 f H N N N2 and
i − | =H E N 2 2 i i i−1 i+1 (N) −H −H pi + H qi(N) = H N N N N 1 1 2 i i τ2 H +o = 2h N N N N2 1 i 1 +o . = 2 σ2 H N N N2
(N) Xn+1
2 Xn(N)
Xn(N)
200
PFLUG ET AL.
We know that setting X (N) = X[t(N) N 2] we have X (N)
converges weakly to X(t),
given by (2) (see [16], p. 169) . Birth-and-death processes are particularly simple, since the transition matrix is tridiagonal. As a consequence, the stationary distribution πi(N) (if existing) may be written in a closed form. We have πi(N) = c
i pj(N) −1 i=i0
qj(N)
(12)
,
where c is a norming constant (see [14]). If we set Y (t) = H −1 (X(t)), i.e., X(t) = H (Y (t)), we have by Ito’s formula dY (t) = g Y (t) dt + τ Y (t) dW (t). (13) We check here that the limiting stationary distribution of Yn(N) = H −1 (Xn(N) ) is the stationary distribution of (13). Inserting (10) and (11) into (12), we get πi(N)
i τ2
j −1
=c N 2 j −1 τ i=i0 N =c
i
1−
τ2
+ −
1 g j N−1 N 1 g j N−1 N
j N
i=i0
=c
i
1−
i=i0
= exp
i i=i0
− σ2
j −1
N τ 2 j N−1
1 g j N−1 − N1 g Nj N − N1 g j N−1
−
1 [τ 2 ] Nj − N2 g Nj N τ 2 j N−1 − N1 g j N−1
log 1 −
1 +o N
1 [τ 2 ] Nj − N2 g Nj N τ 2 j N−1 − N1 g j N−1
1 . +o N
If N tends to infinity, such that i/N tends to x, the expression (14) tends to x x c g(u) g(u) 2 du − log τ (x) = 2 exp 2 du c exp 2 2 2 τ (x) 0 τ (u) 0 τ (u) which is the stationary distribution of the diffusion process (13).
(14)
THE AURORA FINANCIAL MANAGEMENT SYSTEM
201
Examples The Vasicek model (3). Here we may choose H (u) = bu. Then h H −1 (x) = 0 h H −1 (x) = b, and therefore f (x) = bg(x),
σ (x) = bτ (x).
In case (3) it is sufficient to set b = σ . We get g(x) = σa (µ − x) and τ (x) = 1. The Cox–Ingersoll–Ross model (4). Here we choose H (u) = bu2 . Then √ h H −1 (x) = 2 bx,
h H −1 (x) = 2b
and we get g(x) =
σ2 a(µ − x) − √ √ 2 bx 8 bx
and
σ τ (x) = √ 2 b
and b must be chosen as σ 2 /4. The exponential Ornstein–Uhlenbeck model (5). Here we choose H (u) = v u . Then h H −1 (x) = x(log v)2 h H −1 (x) = x log v, and get g(x) = (−a log x)/log v. In order to have τ (x) = 1, one must set v = eσ . 5.
Parallel solution methods
Quite some effort has been made in order to develop parallel methods for large scale sparse optimization problems. For some of the more recent references see, e.g., [4,11, 24,27,32,35,36], until recently little had been known about data parallel approaches to those problems. An important subclass of algorithms exhibits the following common characteristics: • they are based on the use of augmented Lagrangians (see [2]) for the primal and/or dual problem, • they are iterative (with a varying number of loop levels), • their converge linearly with the speed of convergence increasing as the strength of the linking constraints decreases, • they only use elementwise dense vector operations, dense vector reductions and sparse matrix–dense vector products during the computations. In the following we shall present perhaps the simplest algorithm from this family.
202
PFLUG ET AL.
5.1. A simple algorithm example The method outlined below, due to [36], solves a decomposable problem with added linear equality linking constraints. It consists of two levels: the multiplier algorithm [2] at the top and the Diagonal Quadratic Approximation (DQA) [23] in the inner loop. The decomposable optimization problem with linking linear equality constraints is stated as: L fi (xi ) (15) min f (x) = i=1
s.t.
L
Ai xi = b,
i=1
xi ∈ Xi ,
i = 1, . . . , L,
where both data and variables are real, Ai ∈ ni ×m and other dimensions conform. With x T = [x1T . . . xLT ], A = [A1 . . . AL ] and a constant ρ > 0, the augmented Lagrangian is 1 (16) 1(x, y) = f (x) + y, b − Ax + ρb − Ax2 . 2 The use of an augmented Lagrangian is motivated by the ability to solve (15) by means of the multiplier algorithm. The simplicity of the multiplier algorithm makes it an attractive candidate for computer implementations. The algorithm below is the basic version of the multiplier method (as used in [36]): Algorithm 1. The basic multiplier algorithm Step 0. Choose the initial value of the multipliers y, required accuracy ε > 0 and fixed penalty ρ > 0; Step 1. For fixed y solve x := argmin x∈X 1(x, y); Step 2. If Ax − b < ε then stop (optimal solution found); Step 3. Update the multiplier vector: y := y + ρ(b − Ax); Step 4. Go to step 1. The minimization in step 1 of algorithm 1 is carried out approximately by means of a parallel method DQA, described below. Steps 2 and 3 use a combination of matrix– vector and vector–vector operations as well as one reduction (vector norm). The decomposable structure of the optimization problem (15) is lost due to the introduction of the quadratic penalty term in the augmented Lagrangian (16). The DQA method replaces the non-decomposable minimization of (16) by L independent (possibly parallel) minimizations of functions of the form 1 x , y) = fi (xi ) − y, Ai xi + ρ b − Ai xi − A x + Ai xi 2 , 1i (xi , 2
(17)
THE AURORA FINANCIAL MANAGEMENT SYSTEM
203
where x is an additional parameter, which is initially set equal to x and then iteratively updated. Algorithm 2. The parallel DQA method Step 0. Step 1. Step 2. Step 3. Step 4.
x := x; For i = 1, . . . , L choose step length factors τi > 0. Set x , y); For i = 1, . . . , L solve xi := argmin xi ∈Xi 1i (xi , If, for all i, Ai (xi − xi ) < ε, then stop (optimal solution found); xi + τi xi ; For i = 1, . . . , L set xi := (1 − τi ) Go to step 1.
Let Ni denote the number of neighbors of block i defined as the largest number of blocks appearing in the same constraint as block i. To eliminate trivially parallel subproblems, it is assumed that ∀ (i) Ni 1. According to [36] τi < 1/Ni is required for convergence, while setting τi = 1/(2Ni ) will maximize the convergence rate. Formulation (15) suggests that ∀(i) Ni = Li , and consequently the upper bound on τi only depends on L. For the method to converge reasonably well the linking constraints must be significantly sparse. Let us now proceed with an outline of the version specialized for linear programs. A linear problem (LP) in standard form may be written as min cT x s.t. Ax = b, x 0.
(18)
Let us denote ith column of matrix A with A∗i and j th row with Aj ∗ . With A∗i treated x and 5x = x − x, as a separate block, z = c − AT y, r = b − A 1 x , y) = zi 5xi + zi xi + ρ r − A∗i 5xi 2 . 1i (5xi , 2 Its closed form unconstrained solution is 5xi =
r T A∗i − zi /ρ . A∗i 2
(19)
xi + τi [ xi + 5xi ]+ , where [·]+ denotes a Step 3 of algorithm 2 becomes xi := (1 − τi ) projection on the nonnegative orthant. The existence of this closed form solution allows a data parallel implementation. It is easy to see that whole calculation inside the loop of algorithm 2 can be performed in parallel for all i = 1, . . . , L by a sequence of simple operations on matrix A and the relevant vectors. 5.2. Stochastic linear problem characteristics and data parallel implementation Assume that each element xi , i = 1, . . . , n one virtual processor is assigned. Let us denote this set of processors with P. The corresponding values xi , 5xi , ci , zi , Ni,i , τi
204
PFLUG ET AL.
and matrix column A∗i would reside on the same processor. Another set D of m virtual processors would hold the sets of values of yj , bj , rj and matrix row Aj ∗ . The duplicate storage of matrix A is likely to be the most efficient way of balancing computation with communication between processors. It is common knowledge that even sequential implementations of the simplex method for linear programming can benefit substantially from such double storage, even though the only efficiency gains come from simpler and less wasteful data access patterns (see, e.g., [5]). Interprocessor communication will result from the interactions between the primal and dual information stored separately. This interplay takes place when a linear operator A (or its adjoint AT ) transfers information from P to D. The sparse matrix–dense vector products AT y and AT r require communication, in which the whole vectors y and r, respectively, would be gathered from processor set D and broadcast to all processors in P, where only the local part of result vectors would be computed. Similarly, calculations of products A x , A( x − x) and Ax will force gathering of the appropriate vectors x, x −x and x from P and their broadcasting to D. Clearly, the sparse matrix–dense vector products will be the most important operations in an efficient parallel implementation. During the calculation of those products the communication is necessitated only by the presence of the nondiagonal blocks in the constraint matrix. In case of our stochastic dynamic problem this means that • only small fragments of dense vectors actually need to be sent over the communication links, • each processor needs at most one contiguous block of data to be sent from a small number of other processors. From the above we conclude that the stochastic dynamic linear programs arising from AURORA financial models will likely be good candidates for efficient data parallel solution. 6.
Conclusions
We have presented how we establish, estimate and optimize the AURORA financial management model. In its present form, it is a large scale linear program. Our objective allows to include a risk aversion factor into the objective. A variant of (1) would be to consider the lower semivariance, i.e., 2 p(n)W (n) − ρ p(n) V − (n) (20) max n∈M
∀(n ∈ M)
n∈M
W (n) − m∈M p(m)W (m) = V + (n) − V − (n), V + (n), V − (n) 0.
This leads to a large scale nonlinear convex program, which can be solved with exactly the same decomposition techniques as the linear one.
THE AURORA FINANCIAL MANAGEMENT SYSTEM
205
Our model may include proportional transaction costs, but not fixed transaction costs. All assets are assume to be divisible. Introducing fixed transaction costs and indivisibles would change the character of the problem into an mixed integer dynamic stochastic problem. Solution techniques for this class of problem are very little developed, the stochastic branch and bound technique (see [27]) seems to be a promising new technique for such problems. References [1] Ph. Artzner, F. Delbaen, J.-M. Eber and D. Heath, Coherent measures of risk, Mathematical Finance 9 (1999) 203–228. [2] D.P. Bertsekas, Constrained Optimization and Lagrange Multiplier Methods (Academic Press, 1982). [3] J.R. Birge, Decomposition and partitioning methods for multistage stochastic linear programs, Operations Research 33 (1985) 989–1007. [4] J.R. Birge and L. Qi, Computing block-angular Karmarkar projections with applications to stochastic programming, Management Science 34(12) (1990). [5] R.E. Bixby, Progress in linear programming, ORSA J. on Computing 6(1) (1994) 15–22. [6] F. Black, E. Derman and W. Toy, A one-factor model of interest rates and its application to tresury bond options, Financial Analysts Journal 3 (January/February 1990) 33–39. [7] M.J. Brennan, E.S. Schwartz and R. Lagnado, Strategic asset allocation, Journal of Economic Dynamics and Control 21 (1997) 1377–1403. [8] D.R. Cariño, T. Kent, D.H. Myers, C. Stacy, M. Sylvanus, A.L. Turner, K. Watanabe and W.T. Ziemba, The Russell–Yasuda Kasai model: an asset/liability model for a Japanese insurance company using multistage stochastic programming, Interfaces 24 (1994) 29–49. [9] J.C. Cox, J.E. Ingersoll and S.A. Ross, A theory of terms of interest rates, Econometrica 53 (1985) 385–407. [10] G.B. Dantzig and P. Wolfe, Decomposition principle for linear programs, Operations Research 8 (1960) 101–111. [11] M.A.H. Dempster and R.T. Thompson, Parallelization and aggregation of nested Benders decomposition, Working paper WP 01/95, The Judge Institute of Management Studies, Cambridge University (1995). [12] J. Dupacova, Stochastic programming models in banking, Working paper, International Institute for Applied Systems Analysis (1991). [13] A. Gupta, G. Karypis and V. Kumar, Highly scalable parallel algorithms for sparse matrix factorization, IEEE Transactions on Parallel and Distributed Systems 8(5) (1997). [14] P. Hoel, S. Port and Ch. Stone, Introduction to Stochastic Processes (Houghton-Mifflin, Boston, 1972). [15] P. Kall and S.W. Wallace, Stochastic Programming (Wiley, Chichester, 1994). [16] S. Karlin and H. Taylor, A Second Course in Stochastic Processes (Academic Press, Boston, 1993). [17] H. Konno and H. Yamazaki, Mean absolute deviation portfolio optimization model and its applications to Tokyo stock market, Management Science 37 (1991) 519–531. [18] M.I. Kusy and W.T. Ziemba, A bank asset and liability management model, Operations Research 34 (1986) 356–376. [19] H. Markowitz, Portfolio selection, Journal of Finance 7 (1952) 77–91. [20] J.M. Mulvey, Nonlinear network models in finance, in: Advances in Mathematical Programming and Financial Planning (JAI Press, 1987). [21] J.M. Mulvey, An asset–liability investment system, Interfaces 24 (1994) 22–33. [22] J.M. Mulvey, Multi-stage financial planning systems, in: Operations Research Models in Quantitative Finance, eds. R.L. D’Ecclesia and S.A. Zenios (Physica-Verlag, 1994) pp. 18–35.
206
PFLUG ET AL.
[23] J.M. Mulvey and A. Ruszczy´nski, A diagonal quadratic approximation method for large-scale linear programs, Operations Research Letters 12 (1992) 205–215. [24] J.M. Mulvey and A. Ruszczy´nski, A new scenario decomposition method for large-scale stochastic optimization, Operations Research 43 (1995) 477–490. [25] J.M. Mulvey and W.T. Ziemba, Worldwide Asset and Liability Modeling (Cambridge University Press, 1998). [26] J.M. Mulvey and H. Vladimirou, Stochastic network programming for financial planning problems, Management Science 38 (1992) 1642–1664. [27] V. Norkin, G.Ch. Pflug and A. Ruszczy´nski, A branch and bound method for stochastic global optimization, Mathematical Programming 83 (1998) 425–450. [28] W. Ogryczak and A. Ruszczy´nski, From stochastic dominance to mean-risk models: Semideviations as risk measures, European Journal of Operations Research 116 (1999) 33–50. [29] G.Ch. Pflug, Risk-reshaping contracts and stochastic optimization, Central European Journal of Operations Research 5(3–4) (1998) 205–230. [30] G.Ch. Pflug, How to Measure Risk? Festschrift to F. Ferschl (Physica-Verlag, 1999). ´ etanowski, Dynamic asset allocation under uncertainty for pension fund man[31] G.Ch. Pflug and A. Swi¸ agement, Control and Cybernetics 28(4) (1999). ´ etanowski, Accelerating the regularized decomposition method for two [32] A. Ruszczy´nski and A. Swi¸ stage stochastic linear problems, European Journal of Operations Research 101(2) (1997) 328–342. [33] A. Ruszczy´nski, An augmented Lagrangian decomposition method for block diagonal linear programming problems, Operations Research Letters 8 (1989) 287–294. [34] A. Ruszczy´nski, Interior point methods in stochastic programming, Technical Report WP–93–8, International Institute for Applied Systems Analysis, Laxenburg, Austria (1993). [35] A. Ruszczy´nski, Parallel decomposition of multistage stochastic programming problems, Mathematical Programming 58 (1993) 201–228. [36] A. Ruszczy´nski, On convergence of an augmented Lagrangian decomposition method for sparse convex optimization, Mathematics of Operations Research 20(3) (1995) 634–656. [37] O.A. Vasicek, An equilibrium characterization of the term structure, J. Financial Economics 5 (1977) 177–18. [38] H. Vladimirou and S.A. Zenios, Scalable parallel computations for large-scale stochastic programming, Annals of Operations Research 90 (1999) 87–129. [39] S.A. Zenios, Financial Optimization (Cambridge University Press, 1993). [40] S.A. Zenios and R.L. D’Ecclesia, eds., Operations Research Models in Quantitative Finance (Springer, 1994).
