\documentclass{llncs}
\usepackage{amsmath, amsfonts, amssymb, dsfont}
\usepackage{graphicx}

\title{Algorithmic aspects of scenario-based multi-stage decision process optimization}

\author{Ronald Hochreiter}
\institute{Department of Statistics and Decision Support Systems, University of Vienna \email{ronald.hochreiter@compmath.net}}

\begin{document}

\maketitle

\begin{abstract}
Multi-stage decision optimization under uncertainty depends on a careful numerical approximation of the underlying stochastic process, which describes the future uncertain values on which the decision will depend on. The quality of the scenario model severely affects the quality of the solution of the optimization model. Various approaches towards an optimal generation of discrete-state approximations (represented as scenario trees) have been suggested in the literature. Direct scenario tree sampling based on historical data or econometric models, as well as scenario path simulation and optimal tree approximation methods are discussed  from an algorithmic perspective. A multi-stage financial asset management decision optimization model is presented to outline strategies to analyze the impact of various algorithmic scenario generation methodologies.\\ \par
{\bf Key words:} Optimization under uncertainty, scenario generation, scenario optimization, financial decision theory, risk management

\end{abstract}

\section{Introduction}
\label{sec:introduction}

In this paper, we consider algorithmic issues of the computation of multi-stage scenario-based stochastic decision processes for decision optimization models under uncertainty. We consider that given a specific discrete-time stochastic process on the decision horizon $t = 1, \ldots, T$, a decision maker observes the realization of this random process $\xi_t$ at each decision stage $t$, and takes a decision $x_t$ based on all observed values up to $t$ ($\xi_1, \ldots, \xi_t$). Let there now be a sequence of decisions $x_1, \ldots, x_T$. At the terminal stage $T$ we observe a sequence of decisions $x = (x_1, \ldots, x_T)$ with realizations $\xi = (\xi_1, \ldots, \xi_T)$, which lead to cost $f(x, \xi)$ (or likewise profit). The stochastic optimization task is to find the sequence of decisions $x(\xi)$, which minimizes some probability functional (most commonly the expectation, a risk measure, or a combination of these two) of the respective cost function $f(x(\xi), \xi)$ - see e.g. \cite{EichhornRoemisch2005} for a classification of risk measures in this context. 

We consider the multi-stage case in such a way that there is at least one intermediary stage between root and terminal stage, i.e. $T>2$. See \cite{RuszczynskiShapiro2003} for a recent overview of the area of stochastic programming, and especially \cite{WallaceZiemba2005} for stochastic programming languages, environments, and applications. Unfortunately, the topic of scenario generation is not sufficiently treated in most text books on stochastic programming. Formally speaking, we consider multi-stage stochastic programming problems as defined in Eq.~\ref{equ:mssp}.
\begin{eqnarray}
\begin{array}{lll}
\mbox{minimize } x: & \mathds{F} \big ( f(x(\xi), \xi) \big ) & \\
\mbox{subject to} & (x(\xi), \xi) \in \mathcal{X} & \\
& x \in \mathcal{N} &
\end{array}
\label{equ:mssp}
\end{eqnarray}
The multi-variate, multi-stage stochastic process $\xi$ describes the future uncertainty, i.e. the subjective part of the stochastic program, and the constraint set $\mathcal{X}$ defines feasible combinations of $x$ and $\xi$. This constraint set is used to model the underlying real-world decision problem. Furthermore, a set of non-anticipativity constraints $\mathcal{N}$, consisting of functions $\xi \mapsto x$ which make sure that $x_t$ is only based on realizations up to stage $t$ $(\xi_1, \ldots, \xi_t)$, is necessary. To solve multi-stage programs with numerical optimization solvers, the underlying stochastic process has to be discretized into a scenario tree, and this scenario tree approximation will inherently fulfill these non-anticipativity constraints.

Two issues negatively affect the application of multi-stage stochastic programs for real-world decision problems: First, modeling the underlying decision problem is a non-trivial task. Multi-stage optimization models and stochastic scenario models require stable scenario tree handling procedures, which are considered to be too cumbersome to be applied to real-world applications, and the communication of tree-based models to non-experts is complicated. Secondly, modeling the underlying uncertainty is complex and messy. As discussed above, a good discrete-time, discrete-space scenario tree approximation of the underlying stochastic process has to be generated in order to numerically compute a solution, and the quality of the scenario model severely affects the quality of the solution (garbage in $\rightarrow$ garbage out).

In this paper, we focus on the second problem. While different tree generation methods have been proposed, there are hardly any comparisons of different methodologies, see e.g. \cite{DupacovaEtAl2000}, \cite{Pennanen2005}, or \cite{Kuhn2008} for an overview of various methodologies proposed so far. Numerical comparisons as shown by \cite{HilliPennanen2008} and \cite{KautWallace2007} are often focused on a small set of methods and on specific numerical questions. In this paper, a summary of various approaches towards an optimal generation of decision processes is shown to serve as an overview and an outline for various ways to extend and improve the currently applied approaches. 

This paper is organized as follows: In Section \ref{sec:model}, a straight-forward multi-stage financial asset management problem is presented, before any scenario generation methodologies are discussed in order to introduce the application area of stochastic programming to non-experts. This example will be used for numerical studies. Section \ref{ref:treegen} outlines important issues regarding scenario tree generation, while Section \ref{sec:methods} summarizes five different approaches towards tree generation, both direct tree creation techniques as well as  optimal tree approximations of pre-sampled scenario paths. Section \ref{sec:results} presents selected numerical results, while Section \ref{sec:conclusion} concludes the paper.

\section{Multi-stage financial asset management}
\label{sec:model}

To analyze the impact of different scenario generation methodologies we consider a multi-stage asset management decision optimization model. This model is an extension of the classical risk-return optimization approach which was introduced by Markowitz in the early 1950s \cite{Markowitz1952}. In \cite{Hochreiter2008} this approach has been generalized to single-stage stochastic portfolio optimization using a wide range of risk measures (e.g. all coherent risk measures as shown by \cite{ArtznerEtAl1999} among others) to overcome the problem of using the Variance as the only risk measure, and to use scenarios instead of using a pre-estimated correlation matrix, which can be problematic if the number of assets under consideration is large. Furthermore, a detailed multi-stage stochastic formulation of the classical Markowitz Mean-Variance optimization has been reported by \cite{Steinbach2001}. 

Combining the latter two approaches, our decision taker faces a discrete-time decision horizon $t = 1, \ldots, T$, and a set of investment assets $\mathcal{A}$ with uncertain future returns $V_a$, i.e. a stochastic process represented as a multi-variate, multi-stage scenario tree. There is some investment budget $B$ available at each stage up to $T-1$. This amount is deterministically determined in advance. The objective function consists of the aforementioned risk-return bi-criteria functional, whereby the aim is to maximize the expected terminal wealth and to maximize\footnote{Maximizing CVaR is equivalent to minimizing the anticipated financial risk of the stochastic future wealth.} the Conditional Value-at-Risk\footnote{CVaR is also called Expected Shortfall or Tail Value-at-Risk.} (CVaR at confidence level $\alpha$) of the wealth at the terminal stage $T$. Both criteria are weighted using a risk-aversion parameter $\kappa$, which can be adapted to the needs of the investor and to the current market situation. CVaR can be reformulated as a linear programming model conveniently as shown by \cite{RockafellarUryasev2000}, which simplifies the numerical solution procedure.
The main decision is concerned with the amount of budget $b_a$ to be invested into each asset $a$, as the portfolio is rebalanced at the stages $T = 2, \ldots, T-1$. There is no rebalancing at terminal stage $T$. One important constraint is that the amount of purchases $p$ in each stage cannot exceed the sum of the amount of sales $s$ plus the additional budget available at the respective stage.

Given the above problem specification, we may formulate our multi-stage stochastic programming model as shown in Eq. \ref{equ:msam}. The objective function as well as the constraints are shown in stage-based formulation, i.e. the numbers in square brackets represent the stage at which the respective constraint is active. \begin{eqnarray}
\begin{array}{llll}
\mbox{maximize} & \mathds{E} ( \sum_a b_a, T ) +  \kappa \text{CVaR}_\alpha ( \sum_a b_a, T ) & &  \\
\mbox{subject to} & \sum_a b_{a} = B & & [1] \\
&  b_{a} \leq V_a b_{a}^{(-1)} + p_{a} - s_{a} & \forall \mathcal{A} & [2, \ldots, T-1] \\
& \sum_a p_{a} \leq \sum_a s_{a} + B & & [2, \ldots, T-1] \\
&  b_{a} \leq V_a b_{a}^{(-1)}  & \forall \mathcal{A}  & [T] \\
& b_a, p_a, s_a \geq 0 & \forall \mathcal{A} & [1, \ldots, T] \\
\end{array}
\label{equ:msam}
\end{eqnarray}
The multi-stage recourse can be observed in the second and the fourth constraint, as (the scenario tree) $V_a$ represents the future asset return of asset $a$ in the respective stage and is multiplied by the invested budget $b_{a}^{(-1)}$ of the previous stage, which is notated using the $(-1)$ superscript. 

The parameters which have to be specified by the decision taker are the asset returns $V_a$, which are stochastic and need to be approximated, as well as the budget $B$, which is deterministic in this model. The stochastic decision variables, which will be calculated by a numerical optimization solver are the amounts of current investment budget $b_a$, purchases $p_a$, as well as sales $s_a$ of each asset $a$ out of the given investment universe $\mathcal{A}$ at each stage $t$.

\section{Multi-stage scenario tree generation}
\label{ref:treegen}

The quality of the scenario tree severely affects the quality of the solution of the multi-stage stochastic decision model, such that any approximation should be done in consideration of optimality criteria, i.e. before a stochastic optimization model is solved, a scenario optimization problem has to be solved independently of the optimization model. It should be noted, that there are also scenario generation approaches where the generation is not decoupled from the optimization procedure, see e.g. \cite{CaseySen2005}. A major drawback is that these methods are often limited to a certain restricted set of models, and cannot be generalized easily.

In the context of separated scenario optimization, optimality can be defined as the minimization of the distance between the original (continuous or highly discrete) stochastic process and the approximated scenario tree. Choosing an appropriate distance may be based on subjective taste, e.g. Moment Matching as proposed by \cite{HoylandWallace2001}, selected due to theoretical stability considerations (see \cite{RachevRomisch2002} and \cite{HeitschEtAl2006}), which leads to probability metric minimization problems as shown by \cite{Pflug2001} and \cite{DupacovaEtAl2003}, or it may be predetermined by chosen approximation method, e.g. by using different sampling schemes like QMC in \cite{PennanenKoivu2005} or RQMC in \cite{Koivu2005}, see also \cite{Pennanen2009}. It is important to remark that once the appropriate distance has been selected, an appropriate heuristic to approximate the chosen distance has to be applied, which affects the result significantly. 

Single-stage scenario generation, i.e. an optimal approximation of a multi-variate probability distribution without any tree structure can be done via various sampling as well as clustering techniques. The real algorithmic challenge of multi-stage scenario generation is caring about the tree structure while still minimizing the distance. Only in rare cases, this problem can be solved without heuristics. 

Different methods for two main strategies to design multi-stage scenario tree generation heuristics will be presented - direct scenario tree sampling, and optimal tree approximation of pre-simulated scenario paths.

\paragraph{Direct scenario tree sampling.} Based on historical data or some econometric model, a sampling procedure or a node-wise approximation method is applied to build the scenario tree iteratively from the root node to the terminal stage.

\paragraph{Scenario path simulation and optimal tree approximation.} It may be more convenient to use pre-sampled scenario paths and build an optimal approximation out of this data. The time-dependency of the underlying process is preserved, see Fig. \ref{fig:simpath}, and even more importantly any simulation model can be used, which is convenient for real-world applications - either plain econometric time-series models ((V)AR(I)MA, (G)ARCH, \ldots), specialized models (e.g. Wilkie model for actuarial use), or custom (company-specific) scenario generators may be used as input scenario path generator. The optimization task is to find the optimal tree node and arc links. 

\begin{figure}
\begin{center}
\begin{tabular}{cc}
\scalebox{0.3}{\includegraphics{sp1}} &
\scalebox{0.3}{\includegraphics{sp2}} \\
\end{tabular}
\end{center}
\caption{GARCH(1,1) estimation of financial stock returns, sampling $200$ paths.}
\label{fig:simpath}
\end{figure}

One issue, which is often disregarded due to the fact that most research papers on multi-stage scenario tree generation focus on one specific methodology is that different scenario tree methods are restricted to one of the following scenario tree structures, which leads to problems for numerical comparisons: Either the the scenario structure is represented nodes-by-node with a vector denoting the successors of each node in the respective stage, e.g. $N(5, 3, 2)$ results in a scenario tree with $30$ scenarios and $1 + 5 + 15 + 30 = 51$ nodes, or it is represented nodes-at-stage-wise with a vector containing the total number of nodes in each stage, e.g. $S(10, 20, 30)$ results in a scenario tree with $30$ scenarios and $1 + 10 + 20 + 30 = 61$ nodes. The first type is common for techniques, which build the tree from the root to the terminal stage node per node, while the second is inherently used for scenario generation methods, which try to find the optimal per-node branching structure in addition to approximating tree values. 

Another issue regarding multi-stage scenario tree generation is the absence of a common tree format, i.e. no common standard for representing discretized stochastic processes is available, mainly because of the lack of commercial interest. To accommodate all different tree structures, and the implicit formulation of non-anticipativity constraints, a node-based vector/matrix data format of scenario trees is proposed in this paper, which is shown in Table \ref{tab:scentreeformat}. Every scenario generation method presented below produces exactly this representation.
\begin{table}
\begin{center}
\begin{tabular}{ll}
$V(n, d)$ & $d$-dimensional value of node $n$ \\
$A(n)$ & ancestor node of node $n$ \\
$T(n)$ & stage of node $n$ \\
$P(n)$ & probability to reach node $n$ from its ancestor \\
\end{tabular}
\end{center}
\caption{Vector/matrix data format for multi-stage scenario trees}
\label{tab:scentreeformat}
\end{table}

\section{Multi-stage scenario tree generation methods}
\label{sec:methods}

\subsection{Direct scenario tree sampling}

\subsubsection{Iterative node-by-node data resampling} 

If we aim at building a scenario tree with a node-based tree structure $N(n_1, \ldots, n_T)$ given some (multi-variate) times-series, the most straightforward approach is to start at the root node, which represents the current (i.e. deterministic) value, and then to proceed as follows: for each stage $t = 1, \ldots, T$, and for each node in stage $t$ sample $n_t$ random values from the time series using time-series data as well as data from already sampled ancestor nodes.

\paragraph{Example.} If a time-series $v_{a, t}$ consists of three years of daily asset return data, and we want to sample a scenario tree reflecting possible developments for the next three years in three stages, we may sample $n_1$ nodes for the first stage. For each of this $n_1$ nodes, we sample $n_2$ values using a shorter time series, e.g. only the last two years and include each previously sampled value from the first stage into the initial data set for each second stage node. For the third stage we use only one year of data and both ancestor nodes from stage 1 and 2.

While the advantage of this approach is that it is simple, and historical dependence of even highly multi-variate structures can be kept cheaply, it bears the disadvantage that if small trees are needed, two runs might generate completely different trees, and it may be used for a huge amount of stages, if the time-series are not long enough.

\subsubsection{Iterative node-by-node data approximation}

To improve the quality of the solution of the previous method, random sampling of $n_t$ node values can be replaced with an optimal single-stage approximation, e.g. Principal Component Analysis (see e.g. \cite{TopaloglouEtAl2002}), K-Means clustering (see e.g. \cite{HochreiterPflug2007}), or moment matching. The result is affected by the chosen approximation method - both a distance and a heuristic needs to be specified.

\subsubsection{Distribution fitting and direct node sampling}

Sometimes the decision taker wants to use some specific stochastic process for the future development of uncertain values, and is aware of how to estimate the parameters of the respective process. Given a node-based tree structure $N(n_1, \ldots, n_T)$, the tree can be generated from the root node iteratively.

\paragraph{Example.} A univariate process $\xi$ might be normally distributed and follows an additive recursion, i.e. 
\begin{equation}
\xi_0 = \mu_0, \quad \xi_1 \sim N(\mu_0, \sigma_0^2), \quad \xi_t = b \xi_{t-1} + \epsilon_t, \quad \epsilon_t \sim N(\mu, \sigma^2),
\label{eq:process}
\end{equation}
where
$$\mu = \mu_0(1-b), \sigma = \sigma_0(1-b^2),$$
and $\xi_t$ and $\epsilon_t$ are independent. $\xi$ is a stationary Gaussian Markov process. This process is used for numerical results shown below.

In general, this approach is simple and extremely flexible, as any stochastic process which can be estimated and simulated can be used. One disadvantage is the loss of time any dependency, which leads to bad convergence behavior.

\subsection{Optimal tree approximations of scenario paths}

\subsubsection{Scenario merging - forward and backward}

Given a stage-based tree structure $S(s_1, \ldots, s_T)$ and one specific distance function, a straightforward approach to build an optimal tree given a set of pre-sampled scenarios is to start by calculating the distance of each scenario path to each other with the chosen distance. Then proceed in a forward-sweep $t = 1, \ldots, T$ or in a backward-sweep $t = T, \ldots, 1$, and as long as there are more nodes than $s_t$ within the current stage $t$ merge the nodes of the two scenarios, which are closest to each other.

The advantage is that this process is rather straightforward, quite efficient and flexible as the distance calculation is done only once at the beginning, but forward and backward sweeps may lead to different results. Consider the simple simple example shown in Figure \ref{fig:scenmerg}. We like to approximate $4$ sampled scenarios to build a $S(2,3,4)$ tree. In this case there are already differences between using forward and backward-based algorithms, notably not in the terminal stage, but within intermediary stages, which are crucial in any multi-stage setting.

\begin{figure}
\begin{center}
\begin{tabular}{cc}
\scalebox{0.25}{\includegraphics{smf}} &
\scalebox{0.25}{\includegraphics{smb}} \\
\end{tabular}
\end{center}
\caption{Scenario merging example - forward sweep (left), and backward sweep (right).}
\label{fig:scenmerg}
\end{figure}

\subsubsection{Forward-based scenario clustering}

Another method which can be computed efficiently, and can be used for scenario trees with a large amount of stages (e.g. more than 100), is to resemble the node-per-node iterative approximation, but to use the future paths of approximated nodes for subsequent stages.

The implementation is quite tricky, and in the case of a high number of stages non-branching scenarios might occur. The process is schematically shown in Figure \ref{fig:forward}. 

\begin{figure}
\begin{center}
\begin{tabular}{cc}
\scalebox{0.25}{
\includegraphics{forward1}
}
&
\scalebox{0.25}{
\includegraphics{forward2}
}
\\
\scalebox{0.25}{
\includegraphics{forward3}
}
&
\scalebox{0.25}{
\includegraphics{forward4}
}
\end{tabular}
\caption{Forward-based scenario clustering: given pre-sampled paths (top, left), approximate the first stage (top, right), and then use all future paths for each node (bottom, left) to continue the approximations iteratively node-by-node (bottom, right).}
\label{fig:forward}
\end{center}
\end{figure}

\section{Numerical results and comparison}
\label{sec:results}

To obtain numerical results using the multi-stage stochastic financial asset management problem shown in Section \ref{sec:model}, we use two different financial assets. A fixed-income instrument with a deterministic fixed return $r$ at each stage, and a stochastic stock, which will be used for comparing scenario generation techniques. Historical data of the IBM stock has been used, using daily data from January 3rd, 2007 to February 29th, 2008, which is summarized in Figure \ref{fig:ibm}. 

\begin{center}
\begin{figure}
\begin{tabular}{cc}
\scalebox{0.3}{
\includegraphics{ibm1} 
}
&
\scalebox{0.3}{
\includegraphics{ibm2} 
} \\
\end{tabular}
\caption{IBM: historical closing prices (left) and daily returns (right).}
\label{fig:ibm}
\end{figure}
\end{center}

The following additional parameters have been used for solving the asset management model: the initial budget is $B = 10000$ in the first stage, and $0$ for all following stages (t=2,\ldots,T), i.e. budget is only available at the root stage, which does not have economic implications. The risk quantile CVaR $\alpha = 0.9$, and the risk aversion parameter $\kappa = 1$. The fixed-income return is $r = 1.03$ p.a. The return will be replicated for each node of the stock tree (adjusted given the respective stage of the node). 

The node-based optimization model has been implemented using AMPL (implicit formulation of non-anticipativity constraints). All scenario generation methods have been implemented in MatLab 7.6. The optimization models have been solved using Mosek 5 on a MacBook Pro, Mac OS X 10.5.1, Intel, with Intel Core 2 Duo (2.4 GHz), i.e. 1 CPU with 2 Cores, and 2 GB 667 GHz DDR2 SDRAM.

Different aspects have to be analyzed: the run-time of the optimization problem given the specific tree, the convergence of the objective function, and more subtle, the convergence of the optimal decision, i.e. in our specific example the decision on how much budget should be invested into either the stock or the fixed-income product. The visualization of the convergence cannot be done satisfyingly in most cases due to the high dimensionality, such that we will focus on the first two aspects using the stochastic process from Eq. (\ref{eq:process}) above. The financial asset management problem is solved and both the convergence of the objective function (left) as well as the solution time of the optimization solver (right) is shown. The tree structure severely affects the required solution time and the stability of the underlying optimization problem solution. A trade-off has to be taken, which has to be analyzed for each chosen combination of an optimization problem and a scenario generation technique.

The results are shown for two stages in Fig. \ref{fig:conv1}, for three stages using different node structures in Fig. \ref{fig:conv2} and Fig. \ref{fig:conv3}, 
and for four stages in Fig. \ref{fig:conv5}.

\begin{figure}
\begin{center}
\begin{tabular}{cc}
\scalebox{0.25}{\includegraphics{conv1a}} &
\scalebox{0.25}{\includegraphics{conv1b}} \\
\end{tabular}
\end{center}
\caption{Distribution fitting and direct node sampling. $N(k,k) \; \forall k = 2, \ldots, 50$}
\label{fig:conv1}
\end{figure}

\begin{figure}
\begin{center}
\begin{tabular}{cc}
\scalebox{0.2}{\includegraphics{conv2a}} &
\scalebox{0.2}{\includegraphics{conv2b}} \\
\end{tabular}
\end{center}
\caption{Distribution fitting and direct node sampling. $N(k,k,k) \; \forall k = 2, \ldots, 20$}
\label{fig:conv2}
\end{figure}

\begin{figure}
\begin{center}
\begin{tabular}{cc}
\scalebox{0.2}{\includegraphics{conv3a}} &
\scalebox{0.2}{\includegraphics{conv3b}} \\
\end{tabular}
\end{center}
\caption{Distribution fitting and direct node sampling. $N(k,3,3) \; \forall k = 2, \ldots, 40$}
\label{fig:conv3}
\end{figure}

\begin{figure}
\begin{center}
\begin{tabular}{cc}
\scalebox{0.2}{\includegraphics{conv5a}} &
\scalebox{0.2}{\includegraphics{conv5b}} \\
\end{tabular}
\end{center}
\caption{Distribution fitting and direct node sampling. $N(k,3,3,3) \; \forall k = 2, \ldots, 30$}
\label{fig:conv5}
\end{figure}

\section{Conclusion}
\label{sec:conclusion}

In this paper, we presented a unified comparison of different decision process optimization techniques for optimization under uncertainty models. Different generation techniques are described, some of which have not been explicitly described in the literature yet. This description and comparison should create an awareness of choosing some method, and may be used as a guideline for creating and improving new methodologies towards scenario generation, which is the crucial part for solving multi-stage decision problems.

\bibliographystyle{splncs}
\bibliography{adt2009}

\end{document}