COVID-19: data-driven dynamic asset allocation in times of pandemic

1.   Introduction

Nowadays, people, companies and governments in our fast-developing and changing world are increasingly starting to face more situations and problems where they need to take decisions under uncertainty in a multi-period environment. In the financial sector, the ongoing COVID-19 pandemic has demonstrated the lack of robust multi-stage investment strategies in respect to non-stationary returns with high variability. Reflecting the impact of the COVID-19 crisis on the insurance industry, the Swiss Re Group, one of the world's leading providers of reinsurance and insurance, has reported a net loss of EUR 200 mln. for the first quarter of 2020. According to Willis Towers Watson (WTW), a leading global advisory, broking and solutions company, the coronavirus crisis could result in general insurance losses of up to EUR 70 bln. across key insurance classes in the UK and US. On the global scale, the expected monetary loss of global Gross Domestic Product (GDP) ranges from EUR 70.83 bln. to EUR 320.42 bln. dependent on the severity of the scenario. Moreover, even these losses can be underestimated due to the inability of models to accurately reflect the reality.

In this work, we focus on an asset manager who aims to optimize his/her portfolio allocation in a multi-period environment given highly non-stationary historical data on five assets. The historical data spans from 4th January 1999 to 6th July 2020, is on a daily basis and covers the following indices:

$\textbf{US0003M:}$ 3-month US dollar LIBOR interest rate

$\textbf{SBWGU:}$ World Government Bond index in US dollars

$\textbf{LG30TRUU:}$ Barclays Global High Yield Total Return Index

$\textbf{NDDUWI:}$ MSCI World Total Return (Net) Index

$\textbf{NDUEEGF:}$ MSCI Emerging Markets Index.

We consider a multi-period portfolio allocation problem of constrained maximization of the expected CARA utility and we use the multi-stage stochastic optimization as a well-known mathematical tool for the solution of decision-making problems under uncertainty (Shapiro et al. 2009; Pflug and Römisch 2007). We employ numerical solution methods to overcome the fact that the explicit theoretical solution of a multi-stage optimization program may be difficult or even impossible to obtain due to its functional form Pflug and Pichler (2012) and Keshavarz and Boyd (2012). This makes a numerical approach a challenging, important and, very often, irreplaceable method in the field of research.

Our goal in this work is to develop dynamic programming schemes combined with optimal uncertainty approximation techniques. Dynamic programming arises from pioneering papers (Bellman 1956; Bertsekas 1976; Dreyfus 1965), that expressed the optimal policy in terms of an optimization problem with iteratively evolving value function. More recent works (Bertsekas 2007; Keshavarz and Boyd 2012; Hanasusanto and Kuhn 2013; Powell 2007) are built on the fact that the stage-wise estimation of the future value function is a computationally complex procedure. This stresses the necessity to develop numerically efficient dynamic programming algorithms for the solution of multi-stage stochastic optimization problems. Some efficient algorithms are described in the works of (Bertsekas 2007; Hanasusanto and Kuhn 2013; Keshavarz and Boyd 2012; Powell 2007). However, most of the existing methods use a randomized approach in accounting for the stage-wise available information. This may result in a solution bias, especially during a crisis period.

The recent article (Bertsimas et al. 2020) makes a very progressive step towards data-driven approaches in multi-stage stochastic optimization by taking the time dependency structure explicitly into account in a linear case. We follow this path and, in our applied setup, we develop a data-driven dynamic programming scheme for a case with a concave objective function (CARA utility). In our work, optimal allocation decisions for the future are taken in a multi-period environment, given only past time-series data on assets' returns. Clearly, some assumptions about the future of stock returns are necessary for the solution of the problem. Our assumptions concern returns' marginal and conditional distribution functions, for which we consider three scenarios: pessimistic, optimistic and conservative.

$\textbf{Pessimistic:}$ In the pessimistic scenario, we consider distribution estimates which arise based on the COVID-19 time period data and we allow the next wave to happen in the planning horizon.

$\textbf{Optimistic:}$ In the optimistic scenario, we estimate returns' distribution functions based on the data before the COVID-19 pandemic, thus, accounting for returns on pre-crisis levels in our planning horizon.

$\textbf{Conservative:}$ In the conservative case, we use a theoretical approach, given the type of marginal distributions the in-sample data follows. If marginal distributions are Gaussian, the conditional distributions can be expressed in a closed form in line with the well-known work of (Lipster and Shiryayev 1978).

As a fourth scenario, one could account for signal-based nonparametric conditional distribution estimates using COVID-19 early signals to compute means and covariances of future probability models. This, however, is not within the scope of our work due to the lack of historical data on the possible signals. Instead, we work with stochastic processes defined by their continuous-state probability distributions estimated data-based and we assume that these distributions may change over periods in line with their conditional forms (Fort and Pagés 2002; Graf and Luschgy 2000; Pflug 2001). Working, therefore, in a purely distributional setup, we propose two solution schemes: (i) one scheme approximates the initial continuous-state stochastic process by a finitely valued scenario tree prior to an optimization step, while (ii) another uses a data-driven scenario tree approximation directly inside the dynamic programming procedure. Overall, our solution procedure can be subdivided into the following steps:

$\textbf{Step 1:}$ Obtain data-driven estimates of marginal and conditional distribution functions;

$\textbf{Step 2:}$ Approximate the estimated continuous distributions by optimized discrete ones;

$\textbf{Step 3:}$ Solve the underlying optimization problem dynamically,

where Steps 2 and 3 are combined together in the data-driven scheme and where the initial multi-stage stochastic optimization program is formulated in a continuous form while the approximate problem is finite and discrete. The distance between the problems determines the approximation error and, thus, the approximation accuracy.

Next, if a closed-form solution would be available for the multi-stage stochastic optimization problem of interest, Step 1 would directly lead to the optimal decision. However, due to the variational form of multi-stage problems and to the absence of closed-form solutions in a general case, an estimation of continuous-state distribution functions (i.e., Step 1) does not provide the optimal strategy and further approximations are necessary. For the approximation of uncertainty in returns, we focus on optimal quantization methods of continuous-state stochastic processes by finite and discrete scenario trees using the well-known Kantorovich-Wasserstein distance and its optimization, i.e., Step 2 (Kantorovich 1942; Römisch 2010; Timonina 2013; Timonina 2014; Villani 2003). We minimize the Kantorovich-Wasserstein distances between continuous-state distribution functions and their discrete approximations at each stage of the scenario tree.

On one hand, the use of the Kantorovich-Wasserstein distance is motivated by the fact that the amount of information available on scenario trees is crucial for the correct estimation and minimization of the approximation error. In particular, past time-interdependencies and future probabilistic information must be considered: their incorporation would be noisy in case of random sampling (i.e., Monte-Carlo generation). On the other hand, the fundamental result of (Pflug and Pichler 2012) suggests that the approximation error between the optimization problems can be bounded by the term proportional to the sum of stage-wise Kantorovich-Wasserstein distances. Thus, the minimization of this sum over a finite tree structure would lead to the optimal quantization that could provide a fine approximation of the initial problem and, moreover, bound the minimal distance between corresponding optimal values. For example, if the difference in optimal values in wealth-maximization problems could be bounded by 100, - USD, one would be confident that the developed strategy is worse than the true optimal one by no more than 100 USD.

The stage-wise minimization of the Kantorovich-Wasserstein distance between a continuous and a discrete distribution can be performed via existing numerical methods (Pflug and Pichler 2011; Timonina 2013; Timonina 2014; Villani 2003). The question of computational efficiency is clearly of interest in scenario quantization. It is necessary to understand how many values the scenario process should have at each stage of the scenario tree in order to reduce the computational time and to keep the approximation error small. One would wish for small scenario trees and high approximation quality. Importantly, our scenario quantization approach allows to avoid a high number of Monte-Carlo (or Quasi Monte-Carlo) samples which such methods would require, and to enhance the accuracy of the approximation due to the quantization optimality. Furthermore, our data-driven dynamic programming scheme requires much less distribution estimates than the dynamic programming on scenario trees; in particular, the number of distributions grows linearly compared to the exponentially increasing complexity of scenario tree methods. Due to this, our solution method can be applied in a huge variety of areas. Starting from financial planning and inventory control, possible applications include topics such as energy production and trading, electricity generation planning, pension fund management and similar fields.

From the asset manager perspective, we approximate a continuous non-linear asset allocation problem by a discrete one. Solving the problem numerically via our methods, we construct multi-stage optimal asset allocation strategies. In the final part of our work, we compare the profitability of different investment tactics:

$\textbf{Constant tactic:}$ the investor uses the first-stage decision of the three-stage stochastic optimization problem without changing the strategy from month to month;

$\textbf{Reoptimization tactic:}$ the investor changes the weights, reoptimizing the solution from month to month (and, thus, using the first-stage decision each month);

$\textbf{Adaptation tactic:}$ the investor adapts the decision from month to month without reoptimization based on the initial multi-stage solution of the dynamic program.

In general, we consider the problem from the side of an investor who reoptimizes or adapts the decision in each period, as well as from the perspective with no reoptimization (or adaptation) due to high fees or operating costs. Performing our analysis with a link to COVID-19 daily data on new infections, we observe that the reoptimization tactic aids in reducing future wealth uncertainty even in times of crisis (COVID-19). Using the first-stage allocation strategy developed via our method, we observe ca. 10% wealth growth on average out-of-sample with a maximum of ca. 20% and minimum of ca. 5% in three months. Also, our analysis accounts for an increase in volatility of assets, highlighting the relation between the pandemic and the assets' behavior, which is in line with the articles (Hoffmann et al. 2005; Hou 2007; Safvenblad 1997; Lo and Mackinlay 1989; Veronesi 1999), studying the effect of bad news on the stock market and demonstrating a faster diffusion of a negative information (e.g., growth in COVID-19 infections). We show that optimistically offsetted distribution parameters (mean and covariance matrix) lead to a reduction in the out-of-sample wealth.

The rest of the research article is structured as follows. Section 2 describes the mathematical framework, introducing multi-period asset allocation problems. In Section 3 we describe a forward-looking solution procedure on scenario trees and develop dynamic programming schemes with optimal quantizers. Finally, describing our data in Section 4, we propose a three-month optimal allocation strategy adapted to COVID-19 in Section 5.

2.   Problem formulation

Consider the following multi-stage multi-asset portfolio allocation problem maximizing the expected total utility $H({x}, {\xi})$ on a horizon $t = 1, ..., T$

Here, $a_t$ is a time-dependent constant representing the degree of the investor's risk preference. We let $W_0$ be the initial wealth, $W_t$ the investor's realized wealth at period $t = 1, ..., T$, $x = (x_0, ..., x_{T-1})$ the vector of optimal weights and $\xi = (\xi_1, ..., \xi_T)$ the multi-dimensional stochastic process of returns following stage-wise multivariate distributions $\xi_t\sim F_t(\cdot)$. We let $\xi_0$ be the current non-stochastic return and assume it to be equal to zero without loss of generality. We denote process history up to time $t$ by $\xi^t = (\xi_1, ..., \xi_t)$. Further, we introduce constraints ${\xi}\triangleleft \mathcal{F}, \; {x}\triangleleft \mathcal{F}$ incorporating consecutively evolving information. Random returns $\xi_t$ and their optimal weights ${x}_t$ are measurable with respect to $\sigma$-algebra $\mathcal{F}_t\; \forall t$. Notations $\mathbb{0}$ and $\mathbb{1}$ stand for the vectors of zeros and ones correspondingly. The vector $\beta_t$ bounds the percentage allocated to each asset.

Note that both the decision $x_{t}$ and the wealth $W_t$ depend on the history $\xi^t$ of the stochastic process. As its realization is unknown at time $t = 0$, the decision $x_t$ differs for each possible path of $\xi^t$. Therefore, the optimal solution varies in time and is a function of the past uncertainty. To obtain the solution, it is not enough to estimate multivariate distributions $F_t(\xi_t)$. One needs to account for conditional distributions $F_t(\xi_t|\xi^{t-1})$. Given the above, problem (1) is a multi-period decision-making problem under time-interdependent uncertainty (Ermoliev et al. 2004; Shapiro et al. 2009), which can be solved via multi-stage stochastic optimization and approximation techniques.

We use numerical methods for the solution, approximating continuous distribution functions by discrete ones. Our goal is to customize the well-known dynamic programming procedure for the problem (1) enhancing the efficiency of the method via fast data-driven estimates of distribution functions and via an effective approximation of the value function. Overall, a numerical approach is very often an irreplaceable solution method due to the variational form of such problems.

In general terms, problem (1) belongs to the class of multi-stage expectation-maximization stochastic optimization programs given in the form with an objective function $H({x}, {\xi})$ represented as a sum of stage-wise profit/loss functions $h_t(\cdot)$, i.e., $H({x}, {\xi}) = h_0(x_0)+\sum_{t = 1}^{T} h_t(x^{t}, \xi^{t})$ (Pflug and Römisch 2007):

where ${\xi} = (\xi_1, ..., \xi_T)$ is a continuous-state stochastic process defined on the probability space $(\Omega, \mathcal{F}, P)$ and ${\xi}^t = (\xi_1, ..., \xi_t)$ is its history up to time $t$. The random process $\xi$ and the solution $x = (x_0, ..., x_{T-1})$ fulfill the non-anticipativity conditions implying the consecutive evolution of information in time (e.g., exact future returns are not known before the realization occurs). For this, we require the measurability with respect to $\sigma$-algebra $\mathcal{F}_t\; \forall t$, e.g., ${\xi}\triangleleft \mathcal{F}$. Also, $\mathbb{X}$ denotes a general set of constraints on ${x}$.

For the numerical solution of the problem (2), we introduce the approximated problem (3) with the objective function $H({\widetilde{x}}, {\widetilde{\xi}}) = h_0(\widetilde{x}_0)+\sum_{t = 1}^{T} h_t(\widetilde{x}^t, \widetilde{\xi}^{t})$:

Here, the finitely-valued scenario process ${\widetilde{\xi}} = (\widetilde{\xi}_1, ..., \widetilde{\xi}_T)$ approximates the stochastic process $\xi$ (Pflug and Römisch 2007, Pflug 2010, Pflug and Pichler 2011, Pflug and Pichler 2012). The process ${\widetilde{\xi}}$ is defined on a probability space $(\widetilde{\Omega}, \widetilde{\mathcal{F}}, \widetilde{P})$. The distance between problems (2) and (3) determines the approximation error (Pflug and Pichler 2012; Timonina 2013).

For the approximation of the problem (1), we use a scenario process ${\widetilde{\xi}}$ which has $S_t$ samples for each conditional distribution $F_t^d(\widetilde{\xi_t}|\widetilde{\xi}^{t-1})\; \forall d = 1, ..., D_t$ at a given time $t = 1, ..., T$, i.e., $\widetilde{\xi}_t^{s, d}, \; t = 1, ..., T, \; d = 1, ..., D_t, \; s = 1, ..., S_t$. These values are attained with probabilities $p_{t}^{s, d}$ equal to the probability of a distribution multiplied by the probability of a scenario at a particular time. The problem, thus, yields

Note that we could avoid the recursion $W_t^{s, d} = W_{t-1}\widetilde{x}_{t-1}^d\cdot(\mathbb{1}+\widetilde{\xi}^{s, d}_t)$ of the problem (1) replacing it with the constraint $W_t^{s, d} = W_{0}\widetilde{x}_{t-1}^d\cdot(\mathbb{1}+\widetilde{\xi}^{s, d}_t)$ if the stage-wise optimal solution $x_{t-1}$ would not depend on the starting wealth $W_{t-1}$ providing the same fraction for any initial wealth. This is due to the fact that the investor is interested in the optimal allocation but not the optimal value of the problem and, thus, we could keep the investment amount fixed to $W_0$ at each period without influencing the optimal weights. Nevertheless, the constrained maximization of CARA utility does not imply this property being counterintuitive in this respect (a simple example can be found in Appendix 7.1).

3.   Numerical solution method

The approximate problem (4) can be solved via a forward-looking procedure, for which the construction of a scenario tree is necessary before the optimization step (Section 3.1). Scenario trees grow exponentially in the number of nodes, making the method accurate but inefficient. Differently, problem (4) can be solved using the Bellman's principle of optimality starting at the stage $T$ due to the fact that the objective function of the optimization problem (4) is separable in decisions $\widetilde{x}_t\; \forall t$. In this part of the article, we demonstrate how to linearize the growth of scenario trees in the number of distributions given the starting wealth $W_0$ at each period $t$. For this, we use the Bellman's principle of optimality and construct a data-driven scenario tree during the optimization step (Section 3.2). We reduce the approximation error between problems (1) and (4) using the stage-wise minimization of the well-known Kantorovich-Wasserstein distance between continuous and approximate (i.e., discrete) distribution functions. Combining the stage-wise optimal quantization with a backtracking dynamic programming enhances the solution accuracy and keeps the method suitable for high-dimensional cases. In Section 3.3, we consider the solution method on scenario trees with the stage-wise evolution of the starting wealth $W_t$. In the numerical section, we test all proposed strategies for the case when the starting wealth evolves in time and the investor uses the total available amount for the allocation.

3.1.   Forward-looking procedure

The forward-looking approach used for the solution of multi-stage stochastic optimization problems is based on the approximation of stochastic process ${\xi} = (\xi_1, ..., \xi_T)$ by a scenario tree. At the stage $t = 1$, there is only one conditional distribution which coincides with the marginal one, due to the fact that $\xi_0 = 0$ is a non-stochastic observation. Approximating this distribution with the discrete one sitting on $S_1$ points, one assumes $S_1$ conditional distributions at the stage $t = 2$. Further, approximating each of these $S_1$ distributions by discrete ones with $S_2$ points, one assumes $S_2\cdot S_1$ conditional distributions at the stage $t = 3$ and so on. Thus, starting with the stage $t = 1$ and proceeding to $t = T$, one estimates conditional distributions based on previously sampled realization paths (Mirkov and Pflug 2007; Mirkov 2008). Next, one approximates these distributions by discrete ones using scenario quantization methods (e.g., Monte-Carlo sampling, quasi Monte-Carlo sampling, optimal quantization) and acquires a finitely-valued scenario tree (Heitsch and Römisch 2009; Pflug and Pichler 2011) directly used for the numerical solution of the multi-stage stochastic optimization problem. One can use multiple quantization methods in order to approximate the stochastic process ${\xi}$ by a scenario tree:

$\textbf{Monte-Carlo (random) generation}$ randomly selects $S_t$ points from the conditional distribution function $F^d_{t}$ and assigns equal probabilities to them (Fishman 1995);

$\textbf{Quasi Monte-Carlo generation}$ replaces random samples of the Monte-Carlo method by deterministic points uniformly distributed in $[0, 1]^r$, where $r$ is the dimensionality of the stochastic process;

$\textbf{Optimal quantization}$ minimizes the well-known Kantorovich-Wasserstein distance $dl_{KA}\big(F_t, \widetilde{F}_t\big)$ (see Appendix 7.2) between the continuous multivariate distribution $F_t$ and its discrete approximation $\widetilde{F}_t$ (Kantorovich 1942, Villani 2003). The optimal supporting points $\widetilde{\xi}^s$, $s = 1, ..., S_t$ are found by the minimization of a functional $\int\min_{s}d(\xi_t, \widetilde{\xi}^s)F_t(d\xi_t)$ over $\widetilde{\xi}^1, ..., \widetilde{\xi}^{S_t}$ with $d(\xi_t, \widetilde{\xi}^s)$ being an $l_1$- or an $l_2$-distance.

Algorithm 1 describes the stage-wise optimal quantization with optimal supporting points and probabilities demonstrated in Figure 1a. The stochastic process is Gaussian and every conditional distribution of $\xi_t$ given the history ${\xi^{t-1}}$ is also a normal distribution with known mean and variance (Lipster and Shiryayev 1978). For comparison, Figure 1b shows a randomly-quantized scenario tree. Note that the Monte-Carlo generation assigns equal probabilties to the sampled points.

Overall, trees with high branchiness give a finer approximation of stochastic processes. However, the exponential complexity of forward-looking procedures is the result of an increase in the number of scenarios and in the size of the tree. For example, a binary tree with $T = 3$ has 7 conditional distributions and, thus, 7 decision vectors. Differently, a ternary tree (Figure 1b) has 13 conditional distributions and, thus, 13 decision vectors. Note that Figure 1 demonstrates the univariate case with a one-dimensional $\widetilde{\xi}_t$. In case of multidimensionality of the process $\widetilde{\xi}$, vectors would correspond to each node of the tree.

In this article, we use optimal quantization for the distributional discretization due to the fact that it provides a better guarantee for the approximation quality. Under the assumption of Lipschitz-continuity with constants $L_1, ...., L_T$, i.e., if $dl_{KA}(F_t(\cdot|u), F_t(\cdot|v))\leq L_t d(u, v), \; \forall t = 1, ..., T$, the sum of stage-wise Kantorovich-Wasserstein distances establishes an upper bound for the approximation error between problems (2) and (3): the proof is provided in (Pflug and Pichler 2012), which implies

where value functions $v$ and $\widetilde{v}$ correspond to optimal solutions of the multi-stage problems (2) and (3).

3.2.   Data-driven dynamic programming

The idea of dynamic programming goes back to pioneering papers (Bellman 1956; Bertsekas 1976; Dreyfus 1965), expressing the optimal policy in terms of an optimization problem with iteratively evolving value function (the optimal cost-to-go function). More recent works (Bertsekas 2007; Keshavarz and Boyd 2012; Hanasusanto and Kuhn 2013; Powell 2007) are built on the fact that the evaluation of optimal cost-to-go functions involving multivariate conditional expectations is computationally complex. This stresses the necessity to develop numerically efficient dynamic programming algorithms.

For a time-separable optimization problem, the objective function can be written in a way, which allows to separate (to partition) current decision $x_t$ from all previous decisions at stages $(0, 1, ..., t-1)$. Such problems can be solved via multiple algorithms (Bertsekas 2007; Hanasusanto and Kuhn 2013; Keshavarz and Boyd 2012; Powell 2007). However, the suboptimality of existing methods due to a randomized approach in accounting for the information available at stages $t = 1, ..., T$ of the problems may result in under- or overestimation of the optimal value, especially in case of stochastic processes with a heavy-tailed distribution functions or in case of a non-stationary data, as during a crisis. We use stage-wise optimal quantizers and propose an efficient data-driven technique for the solution of the problem (4) with the constraint $W_t^{s, d} = W_{0}\widetilde{x}_{t-1}^d\cdot(\mathbb{1}+\widetilde{\xi}^{s, d}_t)$.

$\textbf{Step 1 - Uncertainty approximation:}$ Estimate conditional distributions directly from historical data. For this, subdivide the data into $T$ parts of length $N$ and, for each of them, estimate $D$ distributions fitting randomly selected $n < N$ data points.

In our case, the time horizon we plan for is equal to $T = 3$ months (90 days) with one period length of $N = 30$ days. First of all, we select the most recent 90-day historical data, subdividing it into three equal parts. For each of the parts, we estimate $D = 30$ distribution functions, randomly selecting the starting date and considering multiple 10-day periods of historical returns in the corresponding data. We optimize $S = 10$ quantizers for each of these distributions. By this, we linearize the complexity fixing the number of distributions at stage $t$ at $D$ instead of $\prod_{i = 1}^tS_{i-1}$.

$\textbf{Step 2 - Dynamic programming at the stage $T-1$:}$ Start with the stage $t = T-1$ and solve the following single-stage optimization problem for each of $D$ conditional distribution functions:

$\textbf{Step 3 - Dynamic programming at the stage $t-1 < T-1$:}$ We estimate the optimal value $\widetilde{V}_{t-1}^d$ via the solution of the following problem $\forall d = 1, ..., D$:

where $\mathbb{E}\big[\widetilde{V}_{t+1}\mid\widetilde{\xi}^{s, d}_t\big]$ is the expectation of the value function conditional on the scenario realization $s$ of the $d$-th distribution at stage $t$. We evaluate $\mathbb{E}\big[\widetilde{V}_{t+1}\mid\widetilde{\xi}^{s, d}_t\big]$ as

where $dl_{KA}^{s, d, i} = dl_{KA}(F_{t+1}^{s, d}, \widetilde{F}_{t+1}^i)$, $\forall i = 1, ..., D$ are Kantorovich-Wasserstein distances between theoretical and data-driven distributions for each scenario $s = 1, .., S_t$, $d = 1, ..., D$ (Figure 2). Theoretical conditional distributions are dependent on the scenario $\widetilde{\xi}^{s, d}_t$ and have a closed form in the Gaussian case. Data-driven distributions correspond to the estimates in Step 1.

In Figure 2, $\widetilde{\xi}_{t-1}$ is a non-stochastic observation at stage $t-1$. There are $S$ optimal supporting points quantizing the data-driven conditional distribution at stage $t$. Furthermore, there are $D$ data-driven distributions at stage $t+1$ and $D$ optimal values $\widetilde{V}^d_{t+1}$ for each of them. One needs to evaluate the expectation $\mathbb{E}\big[\widetilde{V}_{t+1}\mid\widetilde{\xi}^{s, d}_t\big]$ corresponding to the theoretical distribution conditional on each of the optimal supporting points $\widetilde{\xi}^{s, d}_t$ at the stage $t$. We compute $\mathbb{E}\big[\widetilde{V}_{t+1}\mid\widetilde{\xi}^{s, d}_t\big]$ by weighting the values $\widetilde{V}^d_{t+1}$ according to Kantorovich-Wasserstein distances in line with (7). By this, we assign lower weights to the values based on distributions further away from theoretical probability models conditional on $\widetilde{\xi}^{s, d}_t$.

In case of Gaussian-type distributions with mean vector $\mu = (\mu_1, \dots, \mu_T)$ and a non-singular covariance matrix $C = (c_{t_1, t_2})_{t_1, t_2 = 1, ..., T}$, the theoretical distribution conditional on the current state has a closed form and is equal to the distribution (8) (Lipster and Shiryayev 1978):

where $\mu^{t-1} = \big(\mu_1, ..., \mu_{t-1}\big)$ is the unconditional mean process up to time $t-1$ and the submatrix

Algorithm 2 describes the overall dynamic optimization procedure. The total number of data-driven distribution estimates in Algorithm 2 is $D\cdot T$.

Going backwards in time $t = T-1, ..., 1$ in the dynamic programming procedure, one needs to estimate the expected value function at the stage $t+1$. As shown in equation (7), we weight the value function estimates dependent on the distance between two probability distributions. The first distribution is the data-driven approximation $\xi_{t}\mid\xi^{t-1} \sim \mathcal{N}\big(\mu^d_t, \Sigma^d_t)$, while the second one should be the theoretical (true) one and, in general, should coincide with (8) for the Gaussian case. However, to evaluate (8), one requires the complete history of stochastic process quantizers, which is unavailable if not following the forward-looking procedure. Thus, we approximate this distribution using the Markovian structure of means and covariances and assuming that unconditional distribution at the stage $t+1$ coincides with the conditional estimate at the stage $t$ (i.e., $\mathcal{N}\big(\mu^d_t, \Sigma^d_t)$). Thus, we employ the following structure:

where the parameters change their location dependent on the current quantizer and the conditional mean estimate. Here, $\text{cov}(r_{t+1}, r_{t})$ is a covariance estimate for the returns in periods $t+1$ and $t$, while $\hat{\xi}_t^s$, $s = 1, ..., S$ is an optimal quantizer of the Gaussian distribution with mean $\mu^d_t$ and covariance $\Sigma^d_t$.

This sheds light on our approach to estimate distances $dl^{t, i}_{KA}$ during the dynamic programming procedure, starting with the stage $T$ and going backwards in time without explicit information on the historical mean process. We evaluate the Kantorovich-Wasserstein metric between distribution (10) and the data-driven estimate $\mathcal{N}\big(\mu^d_{t+1}, \Sigma^d_{t+1})$.

In Section 3.3 we provide a general dynamic programming scheme, accounting for state variables and incorporating scenario trees. Clearly, this makes the procedure more accurate but much less efficient than the one in Section 3.2. In Section 4, we assess the historical data on several asset classes and analyze it in relation with WHO data on COVID-19 cases, searching for a reaction of returns to the announcement on new infection cases.

3.3.   Dynamic programming on a scenario tree

As stated in Section 3.2, the objective function of a time-separable optimization problem can be written in a way which allows to separate current decision $x_t$ from all previous decisions. If such a partition is not straightforward in time (e.g., wealth $W_{t-1}$ directly depends on $x_{t-2}$, which, thus, influences stage $t$ in the problem (1)), the separation can be done artificially by the introduction of state variables $w_t = \begin{pmatrix}x^{t-1}, \; \xi^{t-1}\end{pmatrix}'$, accumulating all available information on previous decisions and on random component realizations: see (Shapiro et al. 2009) for model state equations for linear optimization. In a general form, optimization problems (2) and (3) can be written as follows:

Note that the state variable $w_t$ may grow in time as $w_t = \begin{pmatrix} w_{t-1}, \; {x}_{t-1}, \; {\xi}_{t-1}\end{pmatrix}'$ $\forall t = 1, ..., T$, describing the accumulation of information (Timonina and Pflug 2017). However, in many practical cases part of the information becomes irrelevant, allowing to neglect some dimensions of the vector $w_t$. If the stochastic process has the Markovian structure, the next value of the process depends on its current value only, being conditionally independent of all the previous values of the stochastic process. Furthermore, some non-Markovian processes can still be represented as Markov chains by expanding the state space. Thus, in line with the recursive wealth constraint in our problem (1), we assume that the dimension of the endogenous variable $w_t$ does not change in time and that the variable obeys the recursion $w_{t+1} = g_t(w_t, x_t, \xi_{t+1})$ $\forall t = 0, ..., T-1$ with the given initial state $w_0$ at $t = 0$. Further, we use Algorithm 1 for the construction of the optimally-quantized scenario tree and combine it with a dynamic programming scheme accounting for the state variable $w_t$:

$\textbf{Step 1 - Uncertainty approximation:}$ Fix the scenario tree structure and quantize conditional distributions optimally as in Algorithm 1. One acquires optimal supporting points sitting at each node of the tree and the corresponding conditional probabilities. Use a grid for the endogenous variable $w_t$, $\forall t = 1, ..., T$. Let us denote points in the grid as $\{\widetilde{w}_t^k\}_{t = 1}^{T}, \; \forall k = 1, ..., K$. Differently, one can use random trajectories for the endogenous state variable or, as in the work (Hanasusanto and Kuhn 2013), one can employ the historical data paths for $w_t$, $\forall t = 1, ..., T$.

$\textbf{Step 2 - Dynamic programming at the stage $T$:}$ Use the scenario tree discretization, as well as the grid for the endogenous variable at stage $t = T$. Let $\widetilde{V}^{k, s, d}_{T}$ be the optimal value computed at the point $(\widetilde{w}^k_{T}, \widetilde{\xi}^{s, d}_{T})$ of the scenario tree. The quantizer $\widetilde{\xi}^{s, d}_{T}$ corresponds to the $d$-th theoretical distribution at stage $T$. We estimate the value $\widetilde{V}_{T}^{k, s, d}$ via the solution of the problem (13) $\forall k, s, d$:

$\textbf{Step 3 - Dynamic programming at the stage $t$:}$ Solving the optimization problem at stage $t+1$, we proceed to stage $t$. Let $\widetilde{V}^{k, s, d}_{t}$ be the optimal value evaluated at the point $(\widetilde{w}^k_{t}, \widetilde{\xi}^{s, d}_{t})$ of the scenario tree. We estimate $\widetilde{V}_{t}^{k, s, d}$ via the following problem $\forall k = 1, ..., K, \; \forall s = 1, ..., S_{t}, \; \forall d = 1, ..., D_{t}$:

where $\widetilde{V}_{t+1}^{k, i, j} = \widetilde{V}_{t+1}(w_{t+1}, \widetilde{\xi}^{i, j}_{t+1})$ and $j$ is the index of the conditional distribution at stage $t+1$. Note that the index $j$ is a function of the realized scenario $s$ and the distribution $d$.

Importantly, in order to solve the optimization problem (14), one needs to evaluate the optimal value $\widetilde{V}_{t+1}^{k, i, j}$ at the point $(g_t(\widetilde{w}^{k}_t, x_t, \widetilde{\xi}^{i, j}_{t+1}), \widetilde{\xi}^{i, j}_{t+1})$. A variety of value function approximation methods are presented in works of (Boyd 2009; Johnson et al. 1993; Nadaraya 1964; Rosen and Marcia 2004; Watson and Geoffrey 1964). However, one needs to guarantee the global solution of a multi-stage optimization problem of interest. In this article, we approximate the value function continuously in $w_{t+1}$ under assumptions about concavity and monotonicity of functions $h_t(w_t, x_t, \xi_{t})$, $g_t(w_t, x_t, \xi_{t+1})$ and $V_{t+1}(w_{t+1}, \xi_{t+1})$ (see Theorems 1 and 2 in Appendix). If convexity (resp. concavity) and monotonicity conditions of Theorems 1 or 2 hold for functions $h_t(w_t, x_t, \xi_{t})$, $g_t(w_t, x_t, \xi_{t+1})$ and $V_{t+1}(w_{t+1}, \xi_{t+1})$ in the dynamic program (14), we can guarantee that the function $V_{t}(w_{t}, \xi_{t})$ is also convex (resp. concave) and monotone. Moreover, these properties stay recursive $\forall t = 1, ..., T$, due to Theorems 1 and 2.

For dynamic programs (13) and (14), Theorems 1 and 2 give the possibility to approximate the optimal value function $V_{t+1}(w_{t+1}, \xi_{t+1})$ by a concave and monotone interpolation in $w_{t+1}$ prior to the solution of the corresponding optimization problem and, therefore, to evaluate the value function at any point $w_{t+1}$, which does not necessarily coincide with grid points $\{\widetilde{w}^{k}_t\}, \; k = 1, ..., K$.

We use the following value functions for our interpolations:

$\textbf{Quadratic:}$ requires monotonicity and concavity constraints and is suitable for general problems;

$\textbf{Linear:}$ does not require constraints to guarantee the global solution, but can be over-simplistic;

$\textbf{Exponential:}$ does not require monotonicity constraints and is especially suitable for the problem (1) due to its form with CARA utility.

Example 1. Let $S = S_t\; \forall t$ and suppose that the optimization problem (14) is evaluated at all points $(\widetilde{w}^k_{t}, \widetilde{\xi}^{s, d}_{t})$ of stage $t$: Table 1 schematically demonstrates these optimal values of the cost-to-go function with respect to subtrees outgoing from the first two nodes at the stage $t-1$. For each scenario tree node $\widetilde{\xi}_t^{s, d}, \; \forall s = 1, ..., S_t, \; d = 1, ..., D_t$, one approximates the function $\widetilde{V}_t(w_t, \widetilde{\xi}^{s, d}_t)$ of $w_t$ by fitting a concave and monotone (under conditions of Theorem 1 or 2) function $\widetilde{V}_t(w_t, \widetilde{\xi}^{s, d}_t)$ to the data in the corresponding column of Table 1. Therefore, one obtains $S_tD_t$ concave and monotone function estimates and, hence, one can evaluate optimal value $\widetilde{V}_{t}({w}_{t}, \widetilde{\xi}^{s, d}_{t})$ at any point $w_{t}$.

The quadratic approximation of the function $\widetilde{V}_{t}(w_{t}, \widetilde{\xi}^{i, j}_{t})$, $\forall i$ is described below:

where $A_i$, $b_i$ and $c_i$ are to be estimated by fitting concave and monotone function $\widetilde{V}_{t}({w}_{t}, \widetilde{\xi}_{t}^{i, j})$ to the points $\widetilde{V}_{t}(\widetilde{w}^k_{t}, \widetilde{\xi}_{t}^{i, j})$, $\forall k = 1, ..., K$.

Importantly, if conditions of Theorem 1 hold, the estimates are obtained via the sum of squares minimization under the constraint implying monotonicity in the sense $w_1\succeq w_2 \Rightarrow \widetilde{V}_t(w_1, \widetilde{\xi}_{t}^{i, j})\geq \widetilde{V}_t(w_2, \widetilde{\xi}_{t}^{i, j})$, i.e.,

where $(w_t)_m$ is the $m$-th coordinate of the vector $w_t$. Differently, if conditions of Theorem 2 hold, the opposite constraint should be used, i.e.,

which implies monotonicity in the sense $w_1\succeq w_2 \Rightarrow \widetilde{V}_t(w_1, \widetilde{\xi}_{t}^{i, j})\leq \widetilde{V}_t(w_2, \widetilde{\xi}_{t}^{i, j})$.

The quadratic function (15) can be computed efficiently by solving the following semidefinite program $\forall i = 1, ..., S_tD_t$:

where $\mathbb{S}^m$ is the set of symmetric matrices and $l$ is the corresponding dimensionality (in our case of CARA-maximization, $l = 1$ as future value functions and future total wealth are directly dependent).

Further, as the wealth at stage $t$ is exponentially related to the value function at stage $t$, we can also use an exponential approximation for the relationship between the wealth at stage $t$ and the value function at stage $t+1$ in the problem (1): this would automatically imply monotonicity if

Also, one does not require monotonicity conditions in case of linear programming (i.e., if functions $h_t(w_t, x_t, \xi_{t})$, $g_t(w_t, x_t, \xi_{t+1})$ and $V_{t+1}(w_{t+1}, \xi_{t+1})$ are linear in $w_t$ and $x_t$). Indeed, linearity conditions are a special case of the requirements of Lemma 2 and they are recursively preserved in the dynamic programming (see Corollary 1). In case of linearity of the program (i.e., in case conditions of Corollary 1 are satisfied), we implement the well-known linear interpolation at the next stage, i.e.,

Algorithm 3 describes the overall dynamic optimization procedure. Also, Figure 3 demonstrates linear and quadratic approximations of the value function in problem (1).

Next, comparing the complexity of Algorithms 2 and 3 one can compute the cumulative number of necessary distribution estimates in both of them. In a logarithmic scale, one observes that the number grows exponentially for the dynamic programming on scenario trees (similar to forward-looking algorithms) while the number is linear in the data-driven scheme.

4.   Data analysis

We employ historical time series data on five assets extracted from Bloomberg: 3-month US dollar LIBOR interest rate (US0003M), World Government Bond index in US dollars (SBWGU), Barclays Global High Yield Total Return Index (LG30TRUU), MSCI World Total Return (Net) Index (NDDUWI) and MSCI Emerging Markets Index (NDUEEGF). The data spans from 4th January 1999 to 6th July 2020 and is on daily basis. Assuming the planning horizon of $T = 3$ months (90 days) to the future, we subdivide the data into in-sample and out-of-sample parts, in order to construct and to test our multi-stage optimal allocation strategies. The in-sample data includes data points on a 3-month horizon until 27th March 2020, while the out-of-sample tests are performed on the basis of data between 28th of March and 6th July 2020. The in-sample data is demonstrated in Figure 5a.

Our in-sample data includes the beginning of the COVID-19 pandemic. Thus, in Figure 5, we observe a strong increase in the volatility of assets, pointing out the relation between the pandemic and the assets' behavior. This goes in line with the works (Hoffmann et al. 2005; Hou 2007; Safvenblad 1997; Lo and Mackinlay 1989; Veronesi 1999), that study the effect of bad news on the stock market and demonstrate a faster diffusion of a negative information (e.g., growth in COVID-19 cases). Denoting daily returns by $r_t$ and computing monthly cumulative returns as $R_t = \prod_{i = 1}^{30}(1+r_{t-i})-1$, we plot monthly values w.r.t. COVID-19 cases in the US in Figure 5b, observing a dramatic drop of returns coinciding in time with the beginning of the COVID-19 world-wide expansion (March 2020). The out-of-sample data does not suggest stabilization of returns to the level before the pandemic. However, to some extent, the variance decreases after April 2020. In Figure 5, we observe a strong interdependency in assets' behavior as negative news arrives.

Clearly, one cannot assume stationarity of the stochastic process based on the assets' behavior in Figure 5. Thus, we allow means, variance and other elements of the covariance matrix to vary in time. To demonstrate the changes in these parameters, we estimate the Gaussian mixture distribution in the rolling interval of 30 days for daily and monthly returns in Figure 6. As expected, we observe a sharp decrease in returns and a strong increase of the covariance matrix elements as soon as negative news about COVID-19 cases arrives. The asset with the smallest variance during the COVID-19 wave is the 3-month LIBOR (US0003M). The increased interdependency between assets and the increased volatility in bad times are align with the studies on diffusion of information (Hoffmann et al. 2005; Hou 2007; Safvenblad 1997; Lo and Mackinlay 1989, Veronesi 1999).

To model the volatility peaks in the dynamic programming procedure, we estimate 30 Gaussian mixture distributions $\mathcal{N}(\mu_t, C_t)$ for each of the months $t = 1, 2, 3$ of our daily in-sample data ($N = 30$ days periods). We model monthly return distributions as normal with means $N\mu_t$ and covariance matrix $\sqrt{N}C_t$ (by this, we assume independence in daily data within each month). However, we account for a correlation structure of returns across different months, assuming that the effects are cumulative in time. The process of modeling monthly correlations is in line with equations (9), (10).

In the next section, we construct single- and multi-stage optimal portfolios accounting for the changes in conditional distributions in our planning horizon ($T = 3$ months). We consider the following three cases for the construction of the multi-stage optimal portfolio:

$\textbf{Pessimistic:}$ Distributions of returns in our planning horizon are equal to the time-dependent estimates provided by our in-sample data. Thus, one expects stabilization of returns in April and May, as well as the second wave of COVID-19 starting in June. As one can observe in the out-of-sample data, this does not happen in reality, suggesting a pessimistic scenario. In general, one could employ epidemiological models using growth in infections as a signal for prediction.

$\textbf{Optimistic:}$ Differently from the first case, one does not expect the second wave of COVID-19 in the planning horizon, assuming the situation to be resolved already after the first wave. This case is optimistic, as it underestimates the variance.

$\textbf{Conservative:}$ Here, we assume only marginal distributions to follow the in-sample data, while conditional distributions are evaluated based on the equation (8).

5.   Optimal portfolio allocation

Consider the optimization problem (1) with $T = 1$. Using the in-sample data on five assets listed in Section 4, we approximate the uncertainty about future daily returns via a Gaussian mixture distribution, finding mean $\mu^1$ and covariance matrix $C^1$. We sample future scenarios for monthly returns quantizing the distribution with mean $N\mu^1$ and covariance matrix $\sqrt{N}C^1$ ($N = 30$ days) by the Kantorovich-Wasserstein distance minimization. We demonstrate the wealth in 100 uncertainty scenarios dependent on the risk-aversion parameter $\alpha_t = \alpha_1$ in Figure 7a. With the initial wealth $W_0 = 1000$, we define the probability of success of our investment strategy as the probability of the total wealth $W_1$ exceeding the initial wealth $W_0$, i.e., as $P(W_1>W_0)$. We plot this probability dependent on the risk-aversion parameter $\alpha_1$ in Figure 7b. We observe low variation in the resulting wealth dependent on the parameter $\alpha_1$. Due to this and as the scenarios have different probabilities with a Gaussian form, we fix $\alpha_t = 0.5\; \forall t$ for our further computations.

The optimal allocation obtained via a single-stage problem is based on distribution functions estimated using the in-sample data on returns. The result is shown in Figure 8a, where the maximal weight is assigned to the 3-month US dollar LIBOR interest rate remaining stable despite COVID-19.

In Section 5.1, we consider the optimization problem (1) with $T = 3$, i.e., a three-stage stochastic optimization problem, which we solve via Algorithms 2 and 3.

5.1.   Data-driven asset allocation

Using Algorithm 2, we approximate the uncertainty about future daily and monthly returns as described in Section 4. We incorporate the dependency between months according to equations (9)–(10) and estimate multiple conditional distributions at each stage of the planning horizon. The optimal weights are shown in Figure 8b, where the variability strongly increases from stage $t = 1$ to $t = 3$ of the planning horizon in line with a natural increase of uncertainty in time. We observe that the amount allocated to the emerging market index is minimal among all other assets and decreases from stage $t = 1$ to $t = 3$. This corresponds to the fact that the first outbreaks of the pandemic took place in emerging countries. Furthermore, the weights at $t = 3$ are evocative of the allocation of the single-stage problem (Figure 8a). This corresponds to the increased volatility of assets during the COVID-19 wave. Differently, the first and the second-stage decisions assign the maximal weight to the World Total Return index leaving only a small fraction to the 3-month LIBOR. This results in ca. 2.5% increase of the total wealth compared to the single-stage strategy tested on the out-of-sample data (Figure 9a).

Maximal achieved wealth of the three-stage problem in Figure 9 is equal to ca. $W_3 = 1200$, while the minimal wealth is ca. $W_3 = 1050$. The difference in possible wealths comes from the fact that we work with a five-dimensional Gaussian mixture distribution whose optimal quantization is not unique. Due to the symmetricity of the distribution, this leads to multiple possible allocations based on the expected wealth maximization. In our analysis, we estimate 30 conditional distributions for each period, quantize each of them with 10 points and find 30 optimal allocations corresponding to the different quantizations, testing the performance of each of them in Figure 9a.

Furthermore, in Figure 9a, we use the first-stage decision of the three-stage stochastic optimization problem and compare its performance to the behavior of the single-stage problem decision without changing the strategy from month to month. We observe that the decision corresponding to the single-stage problem underperforms the first-stage decision of our dynamic program. Thus, our initial recommendation to employ dynamic programming for optimal allocation can be confirmed unless the operating costs for taking such a decision are too high (e.g., salaries to additional analysts).

Further, we consider several possibilities for adapting the optimal weights via our dynamic programming procedure:

i) we can change the weights, reoptimizing the solution from month to month (and, thus, using the first-stage decision each month), or

ii) we can adapt the decision from month to month without reoptimization, based on the initial solution of our dynamic program.

Clearly, case (i) requires the in-sample data growth, while case (ii) is based on the initial in-sample data and does not account for new data realizations.

In general, changing the strategy allows to collect some part of the profit at each period but comes with processing fees. Importantly, any of our estimated strategies can appear to be more profitable in the long run, because of some uncovered and unpredictiable properties of the out-of-sample data. For example, the first-stage dynamic programming strategy demonstrates higher profitability than the strategy (i) in the third month in Figure 9b. Moreover, strategy (ii) outperforms both the first-stage decision and the strategy (i) in this month. Nevertheless, strategy (i) demonstrates a better performance in the second month of our planning horizon. Obviously, Figure 9b corresponds to only one possible allocation, though there are multiple optimal quantizations and different strategies corresponding to them due to multidimensionality of the asset returns process.

Next, studying the performance of 30 estimated strategies, we observe a clear reduction of uncertainty in wealth due to the use of monthly reoptimization (i); see Figure 10. In Figure 10a, each period starting investment is assumed to be the mean among resulting wealths one period before; this allows to test the average performance of our decisions. In Figure 10b, the reinvested amount is the maximal resulting wealth one period before.

We observe that the uncertainty reduction comes at a cost of maximal possible wealth, which, on average, drops from ca. $W_3 = 1200$ to ca. $W_3 = 1150$ (the minimal possible wealth increases from ca. $W_3 = 1050$ to ca. $W_3 = 1100$). In the best case scenario, i.e., if the reinvested amount takes its maximal possible value in each period, the maximal wealth of the reoptimization strategy is comparable with the no-reoptimization case (ca. 0.5% overperformance of no-reoptimization). To take the decision about the employment of these strategies, the asset manager needs to account for the risk tolerance level in the case of a wealth reduction, for the corresponding probability and for the implementation fees.

As mentioned before, the reoptimization strategy allows dividends to be paid (as one does not necessarily invests the total available amount), but is subject to reinvestment fees. Clearly, one also needs to account for operating costs of the reoptimization while making the decision about the strategy implementation. If the cost of working hours or machinery is too high, one can avoid reoptimization, adapting the strategy based on the second-stage decision of the initial dynamic program. In Figure 11, we compare the reoptimization strategy (i) to the strategy (ii), which utilizes adaptation instead of reoptimization. This results in a small increase of uncertainty in wealth while neglecting reoptimization.

Based on the out-of-sample data, we can observe that in our pessimistic scenario, the second wave of COVID-19 was assumed to be earlier than what happened in reality. This is due to the fact that our in-sample data suggested the COVID-19 wave in the third month of the planning horizon. Clearly, in the absence of fine epidemiological models a prediction of COVID-19 signals is a complex task by itself. Nevertheless, from the asset-management perspective, our pessimistic strategy (expecting COVID-19 second wave) appears to be profitable on the out-of-sample data.

Next, we test the offsetted means and covariances of assets according to the optimistic scenario which expects the resolution of the pandemic in the planning horizon. For this, we assume that the parameters stay constant for the whole planning horizon and are estimated based on data in January 2020, i.e., before the first-wave of COVID-19. Solving the problem (1), we observe ca. 1% maximal wealth reduction in case of offsetted estimations (see Figure 12a). Further, neglecting monthly reoptimization can become beneficial in case of the optimistically offsetted parameters as seen in Figure 12b, which suggest ca. 1.3% increase in maximal wealth with optimistic mean and covariance.

5.2.   Asset allocation on scenario trees

Further, we solve problem (1) using Algorithm 3, which is more general than Algorithm 2 but arrives at a cost of high computational efficiency. To perform this method, we construct a five-dimensional scenario tree, which considers theoretical rather than data-driven conditional distributions (Figure 13).

On one hand, such theory-driven scenario trees allow to cover possible variations in data more broadly, avoiding data-driven estimates of conditional distributions and keeping the marginal estimates only. On the other hand, this method can be over-pessimistic and can result in offsetted decisions due to the fact that one would be willing to reduce complexity by decreasing the size of the scenario tree, which, in turn, would also result in the accuracy reduction. Indeed, we observe that the first-stage decision in Figures 14a, c are closely related to the decision of a single-stage optimization problem whose solution is presented in Figure 8. Furthermore, the performance of an out-of-sample wealth is even more uncertain than the performance of wealth in a single-stage problem (Figure 14b, d).

Finally, we consider multiple optimal discretizations and make the second- and the third-stage decisions vary, dependent on both the quantization and the uncertainty realization. We present the allocation in Figure 15 (the first-stage decision varies dependent on the quantization only). Here, we logically observe similar patterns but a higher variance of the solution driven by scenario trees.

6.   Conclusion

In our work, we focus on an asset manager aiming to optimize his/her portfolio allocation in a multi-period environment given highly non-stationary historical data on five assets. To address this goal, we formulate and solve a multi-period portfolio allocation problem of constrained maximization of the expected CARA utility. Using the multi-stage stochastic optimization as a well-known mathematical tool for the solution of decision-making problems under uncertainty, we develop a data-driven dynamic programming scheme for a case with a concave objective function (CARA utility).

Firstly, we show the benefits of the dynamic programming compared to the use of a single-stage optimization problem, which suggests lower out-of-sample returns on average and does not account for interdependencies of returns in time. Further, we consider pessimistic, optimistic and conservative scenarios to account for the time-dependent evolution of probability distribution parameters in time, necessary in multi-stage optimization. Solving the problem numerically via our methods, we construct multi-stage optimal asset allocation strategies and compare the profitability of investment tactics given constant, reoptimized and adapted decision strategies. Thus, we consider the problem from the side of an investor who reoptimizes or adapts the decision in each period, as well as from the perspective of no reoptimization (or adaptation) due to high fees or operating costs. Using the first-stage allocation strategy developed via our method for the whole planning horizon, we observe ca. 10% wealth growth on average out-of-sample with a maximum of ca. 20% and minimum of ca. 5% in three months. We observe that the reoptimization tactic aids in reducing future wealth uncertainty even in times of crisis (COVID-19). Also, we demonstrate that the adaptation tactic, which avoids frequent reoptimization and, thus, allows to save operating costs such as cost of machinery or additional working hours, leads only to a small increase in uncertainty on the resulting wealth. Performing our analysis with a link to COVID-19 daily data on new infections, we account for an increase in volatility of assets pointing out the relation between the pandemic and the assets' behavior. This goes in line with works that study the effect of bad news on the stock market, showing a faster diffusion of a negative information (e.g., growth in COVID-19 infections). To demonstrate the effects of a scenario that does not account for future possible crises, we optimistically offset distribution parameters (mean and covariance matrix) and show that this results in a reduction in the out-of-sample wealth. In our future research, we would like to implement signal-based nonparametric conditional distributions linking COVID-19 early signals with distribution estimates by way of their means and covariances. This will allow better predictions of the evolution of distributional parameters in time, requiring, however, the employment of epidemiological models for COVID-19, not yet present in the literature.

Acknowledgments

We are very grateful to the Austrian Science Fund (FWF Der Wissenschaftsfonds), project No. J3674-N26 for funding our research on data-driven dynamic programming. Furthermore, we would like to thank Prof. Norman Schürhoff at the University of Lausanne and the asset manager OLZ AG (Bern, Switzerland) for very valuable remarks and support of the project's financial application.

Conflict of interest

The author declares no conflicts of interest in this paper.

7.   Appendix

The following statements hold true for a general minimization problem. Similar statements are valid for maximization problems of interest.

Lemma 1. If function $h(w, x)$ is jointly convex in $(w, x)$, then the function $\underset{x}{\operatorname{min}}\; h(w, x)$ is convex in $w$.

Proof 1. For any $w_1$ and $w_2$, let $x_1: = \underset{x}{\operatorname{argmin}}\; h(w_1, x)$ and $x_2: = \underset{x}{\operatorname{argmin}}\; h(w_2, x)$. As the function $h(w, x)$ is jointly convex, the following holds by definition and is in line with our notations:

Therefore, the statement of Lemma 1 follows:

Lemma 2. The following holds true:

1. If $g:\mathbb{R}^n\rightarrow \mathbb{R}^n$ is componentwise convex and $V: \mathbb{R}^n\rightarrow \mathbb{R}^1$ is convex and monotonically increasing as $y_1\succeq y_2 \Rightarrow V(y_1)\geq V(y_2)$, then $V\circ g$ is convex;

2. If $g:\mathbb{R}^n\rightarrow \mathbb{R}^n$ is componentwise concave and $V: \mathbb{R}^n\rightarrow \mathbb{R}^1$ is convex and monotonically decreasing as $y_1\succeq y_2 \Rightarrow V(y_1)\leq V(y_2)$, then $V\circ g$ is convex;

3. If $g:\mathbb{R}^n\rightarrow \mathbb{R}^n$ is linear and $V: \mathbb{R}^n\rightarrow \mathbb{R}^1$ is convex and monotone (i.e. either increasing as $y_1\succeq y_2 \Rightarrow V(y_1)\geq V(y_2)$ or decreasing as $y_1\succeq y_2 \Rightarrow V(y_1)\leq V(y_2)$), then $V\circ g$ is convex.

Proof 2. 1. Proof of the statement 1 of Lemma 2:

As $g:\mathbb{R}^n\rightarrow \mathbb{R}^n$ is componentwise convex, the following holds by definition:

Further, as the function $V$ is monotonically increasing and convex, we claim:

which implies the statement 1 of Lemma 2.

2. Proof of the statement 2 of Lemma 2:

As $g:\mathbb{R}^n\rightarrow \mathbb{R}^n$ is componentwise concave, the following holds by definition:

Further, as the function $V$ is monotonically decreasing and convex, we claim:

which implies the statement 2 of Lemma 2.

3. Statement 3 directly follows from the statements 1 and 2 of Lemma 2.

Theorem 1. Let $\xi_0$ and $\xi_{1}$ be two dependent random variables defined on some probability space $(\Omega, \mathcal{F}, P)$ and let the following hold:

1. Function $h(w, x, \xi_0): \mathbb{R}^n\rightarrow \mathbb{R}^1$ is jointly convex in $(w, x)$ and is monotonically increasing as $w_1\succeq w_2 \Rightarrow h(w_1, x, \xi_0)\geq h(w_2, x, \xi_0), \; \forall x, \xi_0$;

2. Function $g(w, x, \xi_1): \mathbb{R}^n\rightarrow \mathbb{R}^n$ is componentwise convex in $(w, x)$ and is componentwise increasing in $w$;

3. Function $V_{1}(y, \xi_1): \mathbb{R}^n\rightarrow \mathbb{R}^1$ is convex in $y$ and is monotonically increasing as $y_1\succeq y_2 \Rightarrow V_1(y_1, \xi_1)\geq V_1(y_2, \xi_1), \; \forall \xi_1$.

Then the function $V_0(w, \xi_0): = \underset{x}{\operatorname{min}}\bigg\{h(w, x, \xi_0)+\mathbb{E}\big[V_1(g(w, x, \xi_1), \xi_1)\mid \xi_0\big]\bigg\}$ is convex in $w$ and is monotone in the sense of $w_1\succeq w_2 \Rightarrow V_0(w_1, \xi_0)\geq V_0(w_2, \xi_0), \; \forall \xi_0$.

Proof 3. By statement 1 of Lemma 2, function $V_1(g(w, x, \xi_1), \xi_1)$ is jointly convex in $(w, x)$. Further, the function $\big(h(w, x, \xi_0)+\mathbb{E}\big[V_1(g(w, x, \xi_1), \xi_1)\mid \xi_0\big]\big)$ is jointly convex in $(w, x)$ as the sum of two jointly convex functions. Therefore, the convexity result of Theorem 1 follows by Lemma 1.

In order to prove the monotonicity result of Theorem 1, notice, that the following holds:

that implies $w_1\succeq w_2 \Rightarrow V_1(g(w_1, x, \xi_1), \xi_1)\geq V_1(g(w_2, x, \xi_1), \xi_1)$, preserving monotonicity of the function $V_1(g(w, x, \xi_1), \xi_1)$ in $w$. As the minimized sum of two monotone functions, the function $V_0(w, \xi_0)$ is also monotone in $w$ in the sense of $w_1\succeq w_2 \Rightarrow V_0(w_1, \xi_0)\geq V_0(w_2, \xi_0), \; \forall \xi_0$.

Theorem 2. Let $\xi_0$ and $\xi_{1}$ be two dependent random variables defined on some probability space $(\Omega, \mathcal{F}, P)$ and let the following hold:

1. Function $h(w, x, \xi_0): \mathbb{R}^n\rightarrow \mathbb{R}^1$ is jointly convex in $(w, x)$ and is monotonically decreasing as $w_1\succeq w_2 \Rightarrow h(w_1, x, \xi_0)\leq h(w_2, x, \xi_0), \; \forall x, \xi_0$;

2. Function $g(w, x, \xi_1): \mathbb{R}^n\rightarrow \mathbb{R}^n$ is componentwise concave in $(w, x)$ and is componentwise increasing in $w$;

3. Function $V_{1}(y, \xi_1): \mathbb{R}^n\rightarrow \mathbb{R}^1$ is convex in $y$ and is monotonically decreasing as $y_1\succeq y_2 \Rightarrow V_1(y_1, \xi_1)\leq V_1(y_2, \xi_1), \; \forall \xi_1$.

Then the function $V_0(w, \xi_0): = \underset{x}{\operatorname{min}}\bigg\{h(w, x, \xi_0)+\mathbb{E}\big[V_1(g(w, x, \xi_1), \xi_1)\mid \xi_0\big]\bigg\}$ is convex in $w$ and is monotone in the sense of $w_1\succeq w_2 \Rightarrow V_0(w_1, \xi_0)\leq V_0(w_2, \xi_0), \; \forall \xi_0$.

Proof 4. By statement 2 of Lemma 2, function $V_1(g(w, x, \xi_1), \xi_1)$ is jointly convex in $(w, x)$. Further, the function $\big(h(w, x, \xi_0)+\mathbb{E}\big[V_1(g(w, x, \xi_1), \xi_1)\mid \xi_0\big]\big)$ is jointly convex in $(w, x)$ as the sum of two jointly convex functions. Therefore, the convexity result of Theorem 2 follows by Lemma 1.

In order to prove the monotonicity result of Theorem 2, notice that the following holds:

that implies $w_1\succeq w_2 \Rightarrow V_1(g(w_1, x, \xi_1), \xi_1)\leq V_1(g(w_2, x, \xi_1), \xi_1)$, preserving monotonicity of the function $V_1(g(w, x, \xi_1), \xi_1)$ in $s$. As the minimized sum of two monotone functions, the function $V_0(w, \xi_0)$ is also monotone in $w$ in the sense of $w_1\succeq w_2 \Rightarrow V_0(w_1, \xi_0)\leq V_0(w_2, \xi_0), \; \forall \xi_0$.

Corollary 1. Let $\xi_0$ and $\xi_{1}$ be two dependent random variables defined on some probability space $(\Omega, \mathcal{F}, P)$ and let the following hold:

1. Functions $h(w, x, \xi_0): \mathbb{R}^n\rightarrow \mathbb{R}^1$ and $g(w, x, \xi_1): \mathbb{R}^n\rightarrow \mathbb{R}^n$ are linear in $(w, x)$;

2. Function $V_{1}(y, \xi_1): \mathbb{R}^n\rightarrow \mathbb{R}^1$ is linear in $y$.

Then the function $V_0(w, \xi_0) = \bigg(h(w, x, \xi_0)+\mathbb{E}\big[V_1(g(w, x, \xi_1), \xi_1)\mid \xi_0\big]\bigg)$ is linear in $(w, x)$.

Proof 5. Directly follows from Lemma 2.

7.1.   Counterexample

Consider the following single-stage optimization problem, which has a larger feasible region that the problem (1):

This optmization problem has a closed-form solution if a random variable $\xi_T$ is Gaussian, i.e., if $\xi_T\sim\mathcal{N}(\mu, \Sigma)$. Indeed, due to monotonicity of the objective function, the Lagrangian can be written as $\mathcal{L}(x_{T-1}, \lambda) = x_{T-1}^T\mu + 1 - \frac{a_TW_{T-1}}{2}x_{T-1}^T\Sigma x_{T-1} +\lambda(\mathbb{1}'x_{T-1}-1)$, where $\lambda\in \mathbb{R}$ is the dual variable corresponding to the constraint $\mathbb{1}'x_{T-1} = 1$. The optimal solution yields

Both the optimal solution and the optimal value depend on the starting wealth $W_{T-1}$. Importantly, the optimal value has the following form

where $f(W_{T-1})$ is a non-linear function of the starting wealth. Thus, in a multi-stage setup, the optimal solution $x_{T-2}$ at stage $T-2$ would influence both $\exp{(a_TW_{T-1})}$ and $f(W_{T-1})$.

7.2.   Distance between probability measures

To get an intuitive understanding of the distance between probability measures let us start with a description of the well-known Monge Transportation Problem, which was initially formulated in the following way:

Split two equally large volumes into small particles and then associate them with each other so that the sum of products of the paths along which the particles are transported (i.e., the distances between associated particles) to the volume of the transported particle is least. Along what paths must the particles be transported and what is the smallest transportation cost?

Denote the initial and final volumes by $\Omega$ and $\widetilde{\Omega}$ respectively and notice that in general it is not necessary to assume them to be of equal volumes. Instead, they are bodies with equal masses though with not necessarily uniform densities. Let $\pi(\cdot, \cdot)$ be the probability measure on $\Omega\times\widetilde{\Omega}$ describing the shipping plan. Then the fraction of volume of $\Omega$ that is transported from $w$ into $\widetilde{w}$ is defined by $\pi(w, \widetilde{w})$ for any $w\subset \Omega$ and $\widetilde{w}\subset\widetilde{\Omega}$. According to the constraints of the problem, $\pi(w, \widetilde{\Omega}) = \frac{\text{Volume}(w)}{\text{Volume}(\Omega)}$ and $\pi(\Omega, \widetilde{w}) = \frac{\text{Volume}(\widetilde{w})}{\text{Volume}(\widetilde{\Omega})}$. The projections of $\pi(\cdot, \cdot)$ on the first and the second coordinates, denoted by $F(\cdot)$ and $\widetilde{F}(\cdot)$, describe masses of $\Omega$ and $\widetilde{\Omega}$ respectively: $\pi(\Omega, \cdot) = \widetilde{F}(\cdot)$ and $\pi(\cdot, \widetilde{\Omega}) = F(\cdot)$.

If we denote the cost of the transportation from $w$ to $\widetilde{w}$ by $d(w, \widetilde{w})$ then the total cost is $\int_{\Omega\times\widetilde{\Omega}}d(w, \widetilde{w})\pi(dw, d\widetilde{w})$, where $\pi(dw, d\widetilde{w})$ is the mass shipped from the neighborhood of $w$ to the neighborhood of $\widetilde{w}$. Hence, we obtain the so-called Kantorovich formulation of the Monge problem:

Suppose that $F$ and $\widetilde{F}$ are two Borel probability measures given on separable metric spaces $(\Omega, d_{\Omega})$ and $(\widetilde{\Omega}, d_{\widetilde{\Omega}})$ respectively and $\mathcal{P}^{(F, \widetilde{F})}$ is the space of all Borel probability measures $\pi$ on $\Omega\times\widetilde{\Omega}$ with fixed marginals $F(\cdot) = \pi(\cdot, \widetilde{\Omega})$ and $\widetilde{F}(\cdot) = \pi(\Omega, \cdot)$. Evaluate the functional $\mathcal{D}(F, \widetilde{F}) = \inf_{\pi}\bigg\{\int_{\Omega\times\widetilde{\Omega}}d(w, \widetilde{w})\pi(dw, d\widetilde{w}): \pi\in \mathcal{P}^{(F, \widetilde{F})}\bigg\}$, where $d(w, \widetilde{w})$ is the cost function for the transportation of $w\in\Omega$ to $\widetilde{w}\in\widetilde{\Omega}$.

Definition 1 (Kantorovich (multivariate) distance). The Kantorovich distance (Kantorovich 1942) between measures can be defined in the following way:

The Wasserstein distance between measures is a generalization of the Kantorovich distance for order $q\geq1$, i.e.,

Assume now that we search for the approximation of the continuous probability measure $F$ by a discrete one that sits only on $S$ points, i.e., the discrete measure $\widetilde{F}$. The following theorem addresses the question of how the approximation depends on the number of points $S$.

Theorem 3 (Graf and Luschgy 2000). Suppose that an $r$-dimensional distribution with probability measure $F$ has a density $f$ with $\int|u|^{1+\delta}f(u)du < \infty$ for some $\delta>0$ and suppose that $\bigg[\int f(x)^{\frac{r}{r+1}}dx\bigg]^{\frac{r+1}{r}}\leq c$, where $c$ is some constant.

Then $\exists \widetilde{F}$, such that $dl_{KA}(F, \widetilde{F})\leq cS^{-\frac{1}{r}}$, where $\widetilde{F}$ is the discrete approximation of $F$ that sits on $S$ points, and $dl_{KA}(F, \widetilde{F})\rightarrow 0$ for $S\rightarrow\infty$.

Proof. The result of the theorem follows from the Zador-Gersho formula (Graf and Luschgy 2000).

Therefore, increasing the number of points in the approximation, one can guarantee that the Kantorovich-Wasserstein distance converges to zero.

In this article, we approximate five-dimensional distribution functions by discrete distributions, using "Algorithm 1D, $l_1$-distance" from the article (Timonina 2013). For a fixed number of optimal supporting points $S$, we continue iterations of the "Algorithm 1D, $l_1$-distance" until the Kantorovich-Wasserstein distance between two subsequent approximations reaches $10^{-6}$ (see Figure 16 for an example of a lognormal case).

The discrete probability measure $\widetilde{F}$ sits on $S$ optimal supporting points denoted by $\widetilde{{z}} = (\widetilde{z}_1, \widetilde{z}_2, ..., \widetilde{z}_S)^T$, which have corresponding probabilities $\widetilde{{p}} = (\widetilde{p}_1, \widetilde{p}_2, ..., \widetilde{p}_S)^T$, i.e., $\widetilde{F} = \sum_{i = 1}^S \widetilde{p}_i \delta_{\widetilde{z}_i}$.