A Machine Learning Approach to Adaptive Robust Utility Maximization and Hedging

Tao Chen Department of Statistics and Applied Probability, University of California Santa Barbara, 1201 Building 479, Santa Barbara, CA 93106-3110       Michael Ludkovski Department of Statistics and Applied Probability, University of California Santa Barbara, 5520 South Hall, Santa Barbara, CA 93106-3110


Keywords: adaptive robust control, Gaussian process surrogates, portfolio optimization

1 Introduction

Stochastic control formulations have been a core tool in financial mathematics for over 40 years. One fundamental challenge in applying stochastic models to practice is the issue of model risk, i.e. calibrating system dynamics to real life. In this article by “model risk” we mean uncertainty about the underlying probability measure representing the mis-specification between the assumed dynamics and the true probabilistic structure. More precisely, the uncertainty concerns some parameters , so that we have the parameterization . Parameter uncertainty arises in any application, since the required model calibration necessarily leaves some residual inference error about the true parameter values. For example, in problems involving risky assets such as stocks, asset return and volatility are the two main parameters driving the investment decisions, but are notoriously difficult to calibrate.

As one remedy, there have been recently multiple proposals on robust extensions for the underlying stochastic control formulation. The basic strategy is to take a worst-case view among a collection of potential ’s. In the financial context, this conservative perspective stemming from unknown parameters is tempered by the idea of learning the dynamics. Historical data about the risky asset evolution can often be used to improve estimates about its drift and/or volatility. Combining these two concepts of robustness and learning, we adopt the adaptive robust framework recently proposed in Bielecki et al. [9]. This paradigm elegantly connects static robust approaches that do not allow belief updating with adaptive control that fully incorporates the latest estimate within a fixed-parameter setup. Because adaptive control ignores dynamic updating, it is time-inconsistent. However, adaptive robust control also presents formidable numerical challenges as it features an expanded, multi-dimensional state space and non-trivial nonlinear optimization to find the optimal feedback control and the value function . As a result, closed form solutions are ruled out, and standard PDE methods are not feasible beyond the most simple settings.

In this paper we propose and develop a novel algorithm for adaptive robust control. Our approach belongs to the class of Regression Monte Carlo (RMC) and Control Randomization (CR) strategies [10, 5, 12, 14, 26]. The key idea is to recursively construct a functional approximation which is then evaluated over a stochastic (non-gridded) mesh. However, in the context of robust control, existing techniques run into the double challenge of: (i) strong path-dependence of optimal on the control, preventing direct application of Regression Monte Carlo (see Section 4); (ii) the nested optimization due to considering multiple ’s is computationally intensive and requires proper approximation architecture to extract . We employ machine learning techniques to overcome the above challenges, by constructing a non-parametric value function approximation [21, 7]. Specifically, we leverage tools from statistical emulation in Stochastic Simulation [11, 6] (and more broadly the Design and Analysis of Computer Experiments field [25]) by recasting the task of solving the Bellman equation as a statistical learning problem of fitting a surrogate (i.e. a statistical model) for [23, 18]. Our framework resembles recent results for Monte Carlo based solvers of Hamilton-Jacobi-Bellman equations [16, 2, 13]. However, unlike the above references that propose to apply Deep Neural Networks (DNN) for the function approximation step, we rather rely on Gaussian Progress (GP) surrogates. GPs is a core machine learning technique which is well-suited for cases where simulation is expensive and therefore efficient surrogate training is needed (in contrast, DNNs are best suited for high-dimensional problems where maximum flexibility is needed and a lot of training data is available). Related use of GPs is in [3, 18].

Utilizing our numerical algorithm, we provide an extensive investigation into the adaptive robust approach to optimal investment and nonlinear hedging. These are two fundamental problems in financial modeling, yet few fully numerical algorithms have ever been available. Our analysis provides new insights into the interplay between robustifying beliefs, learning and investing/hedging. In particular, we investigate the dependence of the conrols and the terminal wealth on risk- and robustness-level parameters.

The rest of the paper is organized as follows. The rest of Section 1 sets up the adaptive robust formulation and lays out the associated discrete-time Bellman recursions. Section 2 develops our numerical methodology and the Gaussian Process surrogate-based algorithm. The second half of the paper illustrates the algorithm on an Optimal Investment (Section 3) and Hedging (Section 4) case studies. Finally, Section 5 provides discussion on the numerical aspects and concludes with a list of additional enhancements.

1.1 The Adaptive Robust Bellman Equation

Let be a measurable space and a fixed time horizon. Let be the discrete time set and be a non-empty set, which plays the role of the global parameter space, so that for each there is a probability measure on . We write to denote the expectation operator corresponding to . We consider a controlled stochastic Markov state process taking values in the state space . We postulate that this process is observed, and we denote by its natural filtration. The optimization problem involves the family of admissible feedback control processes, which are –adapted processes defined on with taking values in a measurable space and being the feedback control applied at time .

Consider the objective of maximizing a reward functional

with running rewards and terminal reward . In the two financial motivations we have , but since our algorithm is of independent interest, we continue with this slightly more general setup. In the classical setup, the dynamics of are described through a probability measure and dynamic optimization of is reduced to the classical dynamic programming Bellman equation


with and .

The pervasive challenge of prescribing -dynamics implies model mis-specification for determining . Robust control injects consideration of multiple models directly into the optimization step (1) via a generic min-max formulation. Therefore, the single is replaced with many ’s. In the parametric setup we consider, model uncertainty is indexed via whose effect is through the dynamics


The mapping captures the joint effect of the initial condition , the control applied at epoch , the parameterized measure , the stochastic shock and the time-dependence . Note that (2) postulates that we may re-write the evolution of in terms of external, model-independent stochastic shocks; in the common setting we consider, are independent and identically distributed random variables (without loss of generality taken to be standard Gaussian in some ).

The robust control setup then includes an intermediate worst-case optimization over a collection , of measures:


where is the density of the random variable . The expression (3) can be seen as a min-max game between nature that picks and the controller who counteracts with her control . In information-theoretic language, (3) is a game against the adversary who selects a probability measure among the Knightian uncertainty set [15]. Note that the adversarial parameter set is understood to be dynamic, i.e.  is adapted to the information available by period and in particular can depend on .

The adaptive robust control as developed by Bielecki et al. [9], couples the Bayesian updating procedure to the set of measures in (3). Specifically, we take

re-interpreted as the robustified belief set defined by the latest estimate . As a canonical example, is the posterior credible interval at some level (say the 95% CI) centered around the unbiased point estimate . This interpolates dynamic adaptive control, where (no uncertainty but learning), with the robust versions that do not allow for learning, independent of . (We should also mention the myopic adaptive formulation where one plugs-in into the static model that treats as a parameter. This is in fact the most common applied usage which myopically separates learning from control.)

The learning algorithm is meant to reduce the uncertainty about the true probabilistic structure driving . If one assumes that there exists such that for all and the learning and estimation are asymptotically perfect, as then it is also reasonable to expect and the adaptive robust method brings a consistently conservative extension of Bayesian adaptive control. However, the framework can also handle more general situations, such as when never converges to a singleton, or when no exists. Thus, the adaptive robust can handle non-stationary models where the dynamics change over time (as might be imagined happens in real markets, leading to a process ); in such settings represents the best current guess about . So the size of need not go to zero.

Returning to the control objective, since now affects the inner optimization in (4), it is augmented to the system state

The dynamics of consist of the autonomous dynamics of , as well as the filtering equations that govern the evolution of . We summarize them by the mapping

which aggregates the previous dynamics of the state process with the filtering equations for . Note that those joint dynamics are degenerate since a single noise term simultaneously drives both and . As a result, there is a strong dependence between and , and more generally , which is to be contrasted with the static robust framework where and are independent. The resulting adaptive robust Bellman equation which is the main object of analysis in this article is then


Note that now simply denotes expectation over which is the sole stochastic, model-independent component above; the parameter dependence and model risk are encoded in . The solution entails applying backward induction over to compute .

1.2 Solving the Bellman Equation

To implement the adaptive robust framework, we must solve (5) which practically reduces to (4) with the state-dependent uncertainty set as in (13). If we define the operators via

then the discrete time Bellman recursion is . is the propagation operator that computes the continuation value based on taking conditional expectations of step-ahead value function. is the optimization operator that computes the optimal control at , recording the optimizers . is the worst-case parameters (i.e. beliefs) for the model dynamics given and is the robustified feedback control in state . The difference vs is precisely the precaution taken by the controller against model mis-specification.

Three challenges must be handled to solve (5):

  • Continuous state space: because the state space of is continuous and multi-dimensional, some discretization is required to handle it computationally. This discretization however implies that (i.e. ) will not be computed for all potential , but only for some finite subset. Therefore, space discretization must be accompanied by interpolation in order to be able to evaluate for arbitrary in subsequent steps.

  • Integration: the integral over in (4) requires approximation given that the integrand is not analytically available. This implies that the true conditional operator is to be replaced with an approximation , i.e.  is replaced with .

  • Optimization: finding the optimizers and is generally not possible analytically, so another approximation is required to carry out these optimizations numerically. This implies that we replace with an approximation .

In low dimensions, the standard technique for Bellman equations are PDE-based solvers that discretize . However, the augmentation of beliefs to the -state necessarily leads to multi-dimensional setups, essentially ruling out PDE-based strategies that are practical only in dimension . At the same time, simulation-based strategies from conventional regression Monte Carlo are also difficult to apply to (5):

  • The control intrinsically affects the evolution of the state and hence prevents direct simulation of the overall as is done in the standard RMC paradigm.

  • Parametric representation of (e.g. in terms of polynomials in ) is challenging in dimension and brings the concern of overfitting/underfitting;

  • Because we face the additional step of inner optimization over , the representation of is critical to yield financially reasonable/accurate estimates of , which is another limitation of classical parametric approximations;

  • One usually learns by constructing a grid of -values (akin to finite-difference approaches for PDEs); however such gridding is extremely inefficient for and essentially impossible for .

As a result, fully numerical approaches to robust stochastic control have been traditionally viewed as intractable and there remains a large gap in actual use of robust models for financial modeling. In this article we resolve these challenges, making contributions along two directions. On the algorithmic side, we propose a new, machine-learning-inspired algorithm for (4). The key concept is to employ a non-parametric value function approximation strategy [21, 7], namely a Gaussian process surrogate. In addition, we borrow concepts from the Design and Analysis of Computer Experiments field [25]) to construct the underlying stochastic meshes. This methodology is quite general, and could also be easily modified to tackle other robust control formulations (e.g strong robust), and other contexts (e.g robust optimal stopping). Moreover, it mitigates the issue of scalability to higher dimensions. Related extensions would be treated in separate sequels, and we see a lot of further potential for such blending of RMC and machine learning. On the finance side, we present two detailed case studies of using adaptive robust paradigm in optimal investment and loss-based hedging. In particular, the latter problem (33) was essentially out of reach until now. Employing our algorithm, we are able to give a comprehensive investigation of the resulting strategy and its sensitivity to different model parameters. We also demonstrate that the adaptive robust framework outperforms the common (simpler) alternatives of static robust and adaptive control.

1.3 Motivation: Adaptive Robust Optimal Investment

The problem of optimal investment in risky financial instruments dates back to Markowitz and Merton and has been extensively studied over the past 50 years. Asset return and volatility are the two main parameters driving the investment decisions, but are notoriously difficult to calibrate. Therefore, there is strong interest in methods that explicitly take into account parameter uncertainty, while allowing one to partially learn asset dynamics.

Let be the constant risk-free interest rate. We assume a fixed time grid . The excess log-return of a risky asset between and is Gaussian with mean and variance . Both and are treated as unknown, yielding the framework of (3) with . Thus, the dynamics of the controlled wealth process under are given by


where is the proportion of portfolio wealth invested in the risky asset. Above we restrict the investor from shorting or leveraging her risky position. The objective is to maximize expected utility of terminal wealth

for a utility function . The augmented state is three-dimensional: and the Bellman equation for is


where the expectation is over the standard Gaussian increment .

To specify the transition map , we recall the learning of the asset dynamics. Given past observations of , the respective MLE estimators are denoted by and under satisfy for the recursions [8]


Note that the above are time-dependent, capturing the idea that learning slows over time as the information set grows. As increases, the weight of the prior rises and the role of the latest shock declines.

In sum, the dynamics of under are prescribed by (6)-(9)-(10), all driven by the 1-D exogenous factor . They are summarized by the map


where are treated as external parameters.

Assuming there is a true , the uncertainty set is described through the (1-)-confidence region for for some confidence level . In the case where only the drift is uncertain and is given, is an interval centered on ,


where denotes the -quantile of a standard normal distribution. In the 2-D case where both and are uncertain, it is shown in [8] that the -confidence region for at time is an ellipsoid given by


where is the -quantile of the distribution with two degrees of freedom.

When the utility function is of the CRRA power type, , then the problem (7)-(8) is separable with respect to the current wealth : with the dimension-reduced recursion , and for


We will solve (14) in Section 3. Another case study which also relies on is considered in Section 4.

2 Methodology

Our algorithm is based on value function approximation. This is achieved by evaluating the right hand side of (5) at a collection of inputs and then constructing a statistical surrogate for . Thus, we use a machine learning tool, namely Gaussian Process surrogate [22], to learn based on training data , where , see below. To compute the expectation in (14) is approximated with a quantized sum, and the inner saddle-point is computed locally using a nonlinear optimization algorithm. More precisely we use a numerical quadrature representation of the integral over and a nested fminbnd call to get which is interpreted as the estimated optimal control at . Picking the collection is known as experimental design.

The overall solver consists of a fitpredictoptimize loop over time based on the recursion


where predict is needed to evaluate and fit is to construct . For the latter, we evaluate the right hand side of (15) at which yields the pointwise . Because the latter is based on multiple approximations, including approximating the true , approximating the true saddle-point, approximating the true , we treat them as noisy versions of and thus the surrogate fitting step includes smoothing, rather than solely interpolation. Similarly, we view as pointwise samples from the feedback control map that is also obtained via (an independent) surrogate fitting.

Our algorithm returns the fitted surrogates and for each time step. The interpretation of is as the value function of the underlying adaptive robust control problem. However, for practical purposes we typically wish to know the strategy and resulting utility based on a fixed measure , while (3) assumes that the parameters are dynamically and adversarially drawn from at each time-step. To this end, we concentrate on the outputted feedback control surrogate . To evaluate the resulting utility, we rely on forward Monte Carlo, i.e. we generate out-of-sample forward paths by drawing -realizations based on that induces the controlled trajectories and then evaluate the Monte Carlo estimate

Observe that can also be written as an expectation, but under the path-dependent, “pasted” measure which is defined in terms of the worst-case beliefs at each time step and does not admit a simple interpretation. Indeed, the expectation of is carried out with respect to which depends on . For this reason, the representation is generally not practically relevant.

Remark 2.1.

For an arbitrary , we distinguish between the actual minimizer that is obtained by calling predict to evaluate and then optimizing, and the surrogate-based that is obtained from a statistical prediction. Because optimization is expensive, it is much faster to utilize . We also stress that both and are defined based on and hence will necessarily differ from the true optimal feedback due to error back-propagation.

Remark 2.2.

When the optimization over is over a compact space such as , an alternative to generic fminbnd or optim gradient-free solvers is to employ direct search over a discrete candidate set such as . This is faster but of course introduces a fixed discretization error relative to the true optimum. Conversely, the surrogate can output the gradient that could be combined with the chain rule to enable gradient-based optimization.

Before proceeding to discuss the details of operationalizing the above method, we summarize the original proposal of Bielecki et al. [9] that relied on a Control Randomization (CR) + Regression approach. Bielecki et al. started with a Monte Carlo-based paradigm that works with trajectories of and exploited the fact that the underlying state is an autonomous process. In their approach, one begins by generating trajectories of . In addition, one fixes a measure and uses the above trajectories to construct the corresponding filtered posterior estimates as in (9)-(10). Note that is picked arbitrarily and serves as a “baseline” probability measure like in the CR method [1, 17]. Next, the Bellman equation (3) is solved pathwise along the above as follows. Given , one solves for by replacing the integral with a weighted sum


where the Gaussian quadrature recipe specifies the appropriate knots and corresponding quadrature weights . The latter are optimized to minimize a certain global criterion [4, 20]. Since the noise distribution is fixed throughout, the optimal quadrature recipe is pre-computed and stored offline. In the second sub-step, the one-step-ahead values are approximated—since they will not fall on the mesh —through linear interpolation of the training collection .

The above recipe has several limitations that our new proposal overcomes. First, the use of a pre-specified leads to a non-adaptive experimental design, which puts potentially very far from the likely values of (depending on whether is “good” or not). By breaking up the path-based construction, we are able to adaptively build better designs, which in particular leverage the structure of . Second, linear interpolation of ’s brings three restrictions: (i) it poorly scales in dimension , because the interpolation requires sorting of all the sites to identify nearest neighbors; (ii) it implies that is piecewise linear which a restrictive approximation architecture (in particular non-smooth which may cause problems when solving for ); (iii) it takes the estimated values of are taken as exact, providing no ability to smooth previous approximations. Third, the procedure only outputs the pointwise estimates ; there is no simple way to obtain for arbitrary , essentially ruling out out-of-sample estimation under any other measure than . Fourth, the approach has no way to avoid extrapolation when using (16); linear extrapolation is prone to significant errors that are likely to back-propagate. Fifth, pre-simulation under requires to use the same mesh size across all time steps.

2.1 Gaussian Process Surrogates

Gaussian process models is a popular choice for flexible statistical models. They allow a consistent treatment of interpolation and regression with very few tunable hyperparameters. The idea of GPs is to view the function to be approximated as a realization of a Gaussian random field with covariance kernel . Training the model then reduces to applying the Gaussian conditional equations, i.e. evaluating the distribution of the random field given the data. The Gaussian structure implies that the conditional process is still Gaussian, so its distribution is summarized through the posterior mean and posterior (co-)variance , interpreted as the model-predicted approximation and the corresponding posterior uncertainty or “standard error”.

Intuitively, the GP offers a methodology to put a smooth surface through the data labelled generically as . For numerical stability, we introduce the “nugget” , so that rather than exactly interpolating the given ’s, the surrogate also smoothes. Namely, the nugget is equivalent to assuming

In our context, adding ’s is justified via the errors coming from quantization and numerical optimization and we use . The predicted at any input is given by

with corresponding posterior covariance

with the vector and matrix defined by


Above is a identity matrix and is the noise matrix.

To construct a GP surrogate requires selecting a kernel family and the corresponding hyper-parameters. The latter step is typically done through Maximum Likelihood Estimation, i.e. a nonlinear optimization problem involving the respective likelihood of observations. A common kernel family we use is the Matern-5/2,


where the lengthscales and the process variance are estimated via maximum likelihood. Specifically, we use the GPML Matlab package to obtain the MLE for the hyperparameters . Another kernel is the squared-exponential


The GP predictive surface is akin to kernel regression in the sense that the prediction (17) at is always a weighted average of the ’s, with the weights driven by the spatial covariance structure encoded in . In particular, if the spatial correlation decays quickly, the surface will tend to be more “bumpy”, while for very strong spatial correlation the surface will be nearly flat/linear. The differentiability of the predictive surface is driven by the properties of . Under the Matern-5/2 choice (20), the resulting is twice-differentiable, while for (21) it is . In both families, the lengthscales determine the spatial “wiggliness” of the fitted surface in the respective coordinate. Note that different lengthscales (anisotropy) in different coordinates allow for to be, say, more flat in , but more flexible in .

For the propagation operator showing in the Bellman equation we proceed as in (16), approximating the integral over a standard normal distribution by a -points optimal quantizer , i.e. a discrete sum over quadrature points .

Remark 2.3.

In CR-RMC [1, 17], solving the Bellman equation

is handled by defining the q-value across , building a surrogate and then finding as the analytic maximizer of . In contrast, in our approach, we first fix and maximize at to obtain , then fit a surrogate , effectively reversing the order of optimization and emulation.

2.2 Experimental Design

The quality of the regression estimates is critically linked to the choice of the simulation design . Specifically, accuracy of the GP-based at some given input is directly related to the density of the design around . This is the localization property of the GP and matches the intuition of constructing an interpolant for the training data : the interpolated prediction would be good in the neighborhoods of ’s, and progressively worse far away. In particular, as gets far from the training inputs , the GP prediction is driven by the chosen asymptotes of reflecting extreme extrapolation. Specifically, for GPs we have that reverts to its prior mean.

Based on the above discussion, the design should concentrate on the region of interest. In our context, the latter is driven by the pairs that are likely to be encountered by the controller at step . In addition, since the key output we seek is the investment strategy , the hardest learning task is to identify the optimal control when the constraints are not binding. Consequently, the region of interest is .

Another problem that requires attention is extrapolation in estimating the surrogate . Extrapolation tends to lead to large estimation errors in approximating the true conditional expectation operator . With this in mind, we need to select the design so that extrapolation is minimized. In order to compute , we must consider at any within . This, in turn, means that we need good estimation of , , where recall that comes from Gaussian quantization. Figure 1 displays the relationship between and the set of all next-step locations . These are the locations where we must evaluate and hence highlights to what extent the method requires extrapolation. On the right panel we compare all needed predictive sites against the used for training .

Given a region of interest , there are two main concepts to generate a training set that targets . The first idea is a “density-based” approach that aims to make mimic the distribution of under . Recall that includes the exogenous and the path-dependent . One method to handle this path-dependency is based on Control Randomization [1, 17] which was also the motivation for generating the path-simulations in Bielecki et al. [9]. We first pick some and generate paths under . In other words, we pick some “true” model parameters (e.g. using the prior mean estimates ) that are used to generate . A collection of such independent “pilot” paths can be then used to set .

The second concept is to apply space-filling schemes which aim to yield a “uniform” sample from the region of interest. This approach is similar to classical maximin/max-entropy experimental designs in statistics. Space filling is achieved by specifying an input domain (usually some polyhedral or rectangular bounding box) and then a space-filling sequence. A popular choice are Quasi Monte Carlo (QMC) sequences that are used as a variance reduction tool to ensure that does not contain any “holes” that might degrade estimation in their neighborhoods. A QMC sequence offers a deterministic way to sequentially fill the unit hypercube and can be straightforwardly scaled and clipped to fill any polyhedral domain. For example, Sobol QMC sequences offer a discrete of any size .

In Figure 1, we show an adaptive “mixture” design that we utilize for the portfolio optimization problem with unknown . It is obtained in three steps. In the first step, we simulate forward paths under a pre-specified that yield . We then utilize these “pilot” paths to obtain a convex hull . In the second step, we utilize the Sobol QMC sequence to fill the respective with sites. Finally, in the third step, we augment with another sites based on , i.e. add new sites where the control was non-trivial in the previous time-step.

Figure 1: Constructed simulation design for the portfolio optimization case study. There are a total of 256 sites: Sobol QMC based (200 blue); adaptive based on