Swing contract pricing: with and without Neural Networks
We propose two parametric approaches to evaluate swing contracts with firm
constraints. Our objective is to define approximations for the optimal control,
which represents the amounts of energy purchased throughout the contract. The
first approach involves approximating the optimal control by means of an
explicit parametric function, where the parameters are determined using
stochastic gradient descent based algorithms. The second approach builds on the
first one, where we replace parameters in the first approach by the output of a
neural network. Our numerical experiments demonstrate that by using Langevin
based algorithms, both parameterizations provide, in a short computation time,
better prices compared to state-of-the-art methods.
Vincent Lemaire, Gilles Pagès, Christian Yeo
Swing Contract Pricing: A Parametric Approach with Adjoint Automatic Differentiation and Neural Networks

###### Abstract

We propose two parametric approaches to price swing contracts with firm constraints. Our objective is to create approximations for the optimal control, which represents the amounts of energy purchased throughout the contract. The first approach involves explicitly defining a parametric function to model the optimal control, and the parameters using stochastic gradient descent-based algorithms. The second approach builds on the first one, replacing the parameters with neural networks. Our numerical experiments demonstrate that by using Langevin-based algorithms, both parameterizations provide, in a short computation time, better prices compared to state-of-the-art methods (like the one given by Longstaff and Schwartz).

_Keywords - Swing contracts, stochastic control, stochastic approximation, neural network, Langevin dynamics, greeks._

## Introduction

With the energy market becoming increasingly deregulated, various derivative products have emerged, offering flexibility in delivery dates and purchased energy amounts. "Swing" contracts, also known as "Take-or-Pay" contracts (see [38] for more details), are among the most widely traded contracts in gas and power markets. These contracts allow their holders to purchase amounts of energy on fixed exercise dates, subject to certain constraints. There exist two types of constraints. In the firm constraints setting, the contract holder is not allowed to violate the constraints. Beside this case, there exists an alternative setting where the holder can violate them but will face a penalty proportional to the default. Our paper focuses on the valuation of swing contracts in the firm constraints case.

The valuation of swing contracts is a more challenging task than that of classic American-style contracts [9, 19, 27, 36], mainly due to the presence of constraints related to time and volume. From a probabilistic standpoint, the valuation of swing contracts is related to a Stochastic Optimal Control (_SOC_) problem. Here, the control represents the vector of amounts of energy that can be purchased at each exercise date. To solve this SOC problem, two primary groups of methods have been developed.

The first group concerns methods that are based on the so called "Backward Dynamic Programming Principle" (BDPP). In this group, the price of the swing contract can be written as a solution of a dynamic programming equation (see [1, 2, 3, 20, 24, 38]). More precisely, at each exercise date, given
the underlying price and the cumulative consumption, the contract value is given by the maximum over the possible consumptions of the immediate cash flow plus the (conditional) expected value of future cash flows. The latter amount (which is in fact a conditional expectation) is often called the "continuation value" and its evaluation is the main difficulty when consider BDPP-based approaches. The mostly used method for computing the continuation value is the Longstaff and Schwartz one (see [27]); where the continuation value is approximated as an orthogonal projection over a subspace spanned by a finite number of squared integrable random variables (see [3, 27]). Another method is based on the so called optimal quantization (see [30]) where the underlying (continuous) random variable is approximated by its discrete version which in turns is defined by a "projection" on nearest neighbors or Voronoi cells (see [2]). Thus the continuation value reduces to a conditional expectation of a discrete random variable. One of the major problems with BDPP's approaches is the following. To compute the value of the contract at each exercise date, the maximum is taken over a geometric interval. In practice the latter interval needs to be discretized leading to a precision/complexity tradeoff. Indeed, achieving a high level of accuracy in pricing the contract requires a finer discretization of the geometric interval, which in turn increases computation time. Additionally, the Longstaff-Schwartz method faces a storage challenge, as regression coefficients for each simulation and admissible cumulative consumption must be stored at each exercise date to compute the continuation value. This often leads to a memory overflow when a large number of simulations is required to obtain a price with a tight confidence interval. Meanwhile, the optimal quantization method suffers from the "curse of dimensionality" because it is well known that the rate of convergence of this method is of order \(\mathcal{O}(N^{-\frac{1}{d}})\) (where \(N\) is the number of simulations and \(d\) is the problem dimension).

An alternative approach to BDPP-based methods is to consider the valuation of swing contracts as a global optimization problem. As mentioned, the valuation of swing contract is equivalent to solve a SOC problem where the objective is to find a vector of purchase amounts that maximizes the expected value of cumulative cash flows up to the expiry of the contract. This leads to a stochastic optimization problem where the control is generally substituted by a parametric function, and optimal parameters are identified using Stochastic Gradient Descent (SGD) algorithms. The use of parametric functions to approximate the optimal control reduces the control problem to a parametric optimization problem. However, the success of such approaches is contingent upon finding a suitable parametric function, which often requires a solid understanding and intuition of the optimal control's behavior. Note that, solving SOC problem with SGD based algorithms have been considered in [17] for general SOC problems and in [3] for the swing contracts case. To the best of our knowledge, there has been limited exploration of such approaches in the literature to price swing options. Thus, in this paper, we propose two parametric methods for global optimization and compare them with the Longstaff-Schwartz method. We optimize both parameterizations using two algorithms: Adaptive Moment Estimation (Adam) and (Preconditioned) Stochastic Gradient Langevin (PSGLD). While Adam is widely used in stochastic approximation, PSGLD has demonstrated effectiveness in recent studies for Bayesian learning. Our results show that using Langevin-type algorithms can accelerate the training of our parameterizations and leads to better prices compared to the state-of-the-art methods.

Our paper is organized as follows. Section 1. We describe swing contracts and recall the pricing framework. Section 2. We present our contribution by proposing two numerical methods for the pricing of swing contracts. Section 3. We present generalities about stochastic approximation and study two optimization algorithms to optimize our numerical solutions. Finally in section 4 we present some additional results about the effectiveness of our methods.
On swing contracts

In this section, we will provide an overview of swing contracts and explain how to build the swing volume grid defining the space of volume constraints. In a firm constraints setting, this space is a critical component to price swing contracts, and this will discussed towards the end of this section.

### Description

Swing option is a commonly encountered derivative on energy (gas and power) markets which allows its holder to buy amounts of energy \(q_{k}\) at times \(t_{k}\), \(k=0,...,n-1\) (called exercise dates) until the contract maturity \(t_{n}=T\). At each exercise date \(t_{k}\), the purchase price (or strike price of the contract) is denoted \(K_{k}\) and can be constant (i.e \(K_{k}=K,k=0,...,n-1\)) or indexed on a formula. In this paper, we focus on the fixed strike price case but the case of indexed strike can be treated likewise. In this latter case the strike at a time \(t_{k}\) is often calculated as the average of the underlying prices over a period which is already completed. Therefore the indexed strike price case only adds a computational complexity rather than a theoretical one.

One important feature on swing option are volume constraints. Swing option gives its holder a flexibility on the amount of energy he is allowed to purchase and this flexibility is subject to some (firm) constraints

* **Local constraints**: at each exercise date \(t_{k}\), the holder of the swing contract is allowed to buy at least \(q_{\min}\) and at most \(q_{\max}\) i.e, \[q_{\min}\leq q_{k}\leq q_{\max},\ \ \ \ 0\leq k\leq n-1.\] (1.1)
* **Global constraints**: the cumulative consumption up to the maturity of the contract must be not lower than \(Q_{\min}\) and not greater than \(Q_{\max}\) i.e, \[Q_{n}=\sum_{k=0}^{n-1}q_{k}\in[Q_{\min},Q_{\max}],\ \ \ \text{with}\ \ Q_{0}=0\ \text{and}\ 0\leq Q_{\min}\leq Q_{\max}\leq+\infty.\] (1.2)

Note that in this paper we only consider the **firm constraints** case which means that the holder of the contract is not allowed to violate the constraints. This case is often released in the literature by adding a penalty term whenever the global constraints are violated (see [2, 3]). In penalty setting, the existence of an optimal Markov consumption had been proved in [3] using a smooth penalty function. A similar result has been established in [1] for the firm constraints case.

Now let us focus on the firm constraints case in which one of the main components is the physical space. That is the area of possible actions taking into account the volume constraints.

### Physical space

At each exercise date, through the firm constraints, the achievable cumulative volumes are delimited by boundaries represented by the functions \(t_{k}\mapsto Q^{down}(t_{k})\) and \(t_{k}\mapsto Q^{up}(t_{k})\)

\[\left\{\begin{array}{l}Q^{down}(t_{0})=0\\ Q^{down}(t_{k})=\max\left(0,Q_{\min}-(n-k)\cdot q_{\max}\right)\ \,k\in\{1,...,n-1\}\\ Q^{down}(t_{n})=Q_{\min}\end{array}\right. \tag{1.3}\]
\[\left\{\begin{array}{l}Q^{up}(t_{0})=0\\ Q^{up}(t_{k})=\min\big{(}k\cdot q_{\max},Q_{\max}\big{)}\\ Q^{up}(t_{n})=Q_{\max}\end{array}\right.,\,k\in\{1,...,n-1\} \tag{1.4}\]

where \(Q^{down}(t_{k})\) and \(Q^{up}(t_{k})\) denote respectively the lower and the upper bound of the cumulative consumption at time \(t_{k}\). These boundaries lead to the physical space of the swing contract, as drawn in Figure 1 and represents at each exercise date, the range of attainable cumulative consumptions.

In Figure 1 we assumed that \(q_{\min}=0\). This assumption can be made without loss of generality because, as shown in [2], the general case can always be reduced to the case where \(q_{\min}=0\) and \(q_{\max}=1\) (see appendix A for the proof). Besides, if at an exercise date \(t_{k}\) we have bought the amounts \(q_{0},...,q_{k-1}\) leading to a cumulative consumption \(Q_{k}=\sum_{i=0}^{k-1}q_{i}\), then due to local constraints and the cumulative consumption borders, the actual range of attainable cumulative consumption at time \(t_{k+1}\) is the following

\[Q_{k}+q_{k}\in\Big{[}\underbrace{\max(Q^{down}(t_{k+1}),Q_{k}+q_{\min})}_{:=L _{k+1}(Q_{k})},\,\,\underbrace{\min(Q^{up}(k+1),Q_{k}+q_{\max})}_{:=U_{k+1}(Q _{k})}\Big{]}. \tag{1.5}\]

That is, starting with a cumulative consumption \(Q_{k}\), at time \(t_{k}\), the range of admissible volume is

\[q_{k}\in\Big{[}\underbrace{L_{k+1}(Q_{k})-Q_{k}}_{:=A_{k}^{-}(Q_{k})},\,\, \underbrace{U_{k+1}(Q_{k})-Q_{k}}_{:=A_{k}^{+}(Q_{k})}\Big{]}. \tag{1.6}\]

With theses basic blocks, we can now set the theoretical framework for the pricing of swing contracts.

### Pricing and sensitivity calculus

Let \((\Omega,\mathcal{F},\{\mathcal{F}_{t}\},\mathbb{P})\) be a filtered probability space. The price at time \(t\) of the forward contract delivered at the maturity \(T\) is denoted by \(F_{t,T}\). In this paper (as in [2, 3]) we consider contract on the spot price (\(S_{t}=F_{t,t}\)) even if in practice the spot is not a tradable instrument on market. Instead, the most encountered gas contract is the day-ahead forward contract whose price is \(F_{t,t+1}\). But this case can be treated in the same way.

The decision process \((q_{k})_{0\leq k\leq n-1}\) is defined on the same probability space and is supposed \(\mathcal{F}_{t_{k}}^{S}\)- adapted. At each exercise date \(t_{k}\), by buying a certain volume \(q_{k}\), the holder of the contract makes a profit (or loss)
\[\psi_{k}\left(q_{k},S_{t_{k}}\right):=q_{k}\cdot\left(S_{t_{k}}-K\right). \tag{1.7}\]

Let \(Q\in\mathbb{R}_{+}\). Given a cumulative consumption \(Q\), an admissible strategy at time \(t_{k}\) is a vector \((q_{k},\ldots,q_{n-1})\) lying within the following set

\[\mathcal{A}_{k,Q}^{Q_{\min},Q_{\max}}=\left\{(q_{\ell})_{k\leq\ell\leq n-1},\;q _{\ell}:(\Omega,\mathcal{F}_{t_{\ell}}^{S},\mathbb{P})\mapsto\left[q_{\min},q_ {\max}\right]\right.\left.\sum_{\ell=k}^{n-1}q_{\ell}\in\left[(Q_{\min}-Q)_{+}, Q_{\max}-Q\right]\right\}.\]

Note that using equation (1.6), the preceding set reads,

\[\mathcal{A}_{k,Q}^{Q_{\min},Q_{\max}}=\left\{(q_{\ell})_{k\leq\ell\leq n-1},\; q_{\ell}:(\Omega,\mathcal{F}_{t_{\ell}}^{S},\mathbb{P})\mapsto\left[A_{\ell}^{-} (Q_{\ell}),A_{\ell}^{+}(Q_{\ell})\right]\right.\text{where}\;Q_{\ell}=Q+\sum_ {i=k}^{\ell-1}q_{i}\right\}. \tag{1.8}\]

with the convention \(\sum_{i=k}^{k-1}q_{i}=0\). Then for every non negative \(\mathcal{F}_{t_{k-1}}^{S}-\) measurable random variable \(Q\), the price of the swing option at time \(t_{k}\) starting from a cumulative consumption \(Q\) is given by

\[P_{k}\left(S_{t_{k}},Q\right)=\operatorname*{ess\,sup}_{(q_{\ell})_{k\leq\ell \leq n-1}\in\mathcal{A}_{k,Q}^{Q_{\min},Q_{\max}}}\;\mathbb{E}\left(\sum_{\ell =k}^{n-1}e^{-r_{\ell}(t_{\ell}-t_{k})}\psi_{\ell}\left(q_{\ell},S_{t_{\ell}} \right)|\mathcal{F}_{t_{k}}^{S}\right), \tag{1.9}\]

where the expectation is taken under the risk-neutral probability and \(r_{\ell}\) are interest rates over the period \([t_{0},t_{n-1}]\) that we will assume to be zero throughout this paper. Then the price of the swing contract is given by (\(S_{t_{0}}\) is assumed to be deterministic)

\[P_{0}:=P_{0}\big{(}S_{t_{0}},0\big{)}=\sup_{(q_{\ell})_{0\leq\ell\leq n-1}\in \mathcal{A}_{0,0}^{Q_{\min},Q_{\max}}}\;\mathcal{J}\big{(}q_{0},\ldots,q_{n-1} \big{)}. \tag{1.10}\]

where, given an admissible strategy \((q_{0},\ldots,q_{n-1})\), the reward function \(\mathcal{J}\) is defined as the expected value of cumulative future cash flows up to the expiry

\[\mathcal{J}\big{(}q_{0},\ldots,q_{n-1}\big{)}:=\mathbb{E}\left(\sum_{\ell=0}^ {n-1}\psi_{\ell}\big{(}q_{\ell},S_{t_{\ell}}\big{)}\right). \tag{1.11}\]

The latter problem appears to be a constrained stochastic control problem in which the aim is to find an admissible strategy that maximizes the reward function \(\mathcal{J}\). As mentioned in the introduction, there exists in the literature two groups of pricing methods to solve this optimization problem. The first group of methods is based on the "backward dynamic programming principle" and the second on a global optimization approach. In this paper, we propose two numerical solutions to solve problem (1.10) based on this latter group.

It is important to note the "bang-bang" feature, which implies that under certain conditions, at each exercise date, the optimal consumption is given by one of two values: \(A_{k}^{-}\) or \(A_{k}^{+}\) (defined in equation (1.6)). This feature has been proven for the case of firm constraints and when global constraints are whole numbers by Pages et al. [1]. For the penalty case, a proof can be found in Gobet et al. [3]. The "bang-bang" feature is particularly valuable since it allows to significantly reduce the computation time it.
In this paper and unless otherwise stated, we use a one-factor model as in [3, 20]. That is,

\[\frac{dF_{t,T}}{F_{t,T}}=\sigma e^{-\alpha(T-t)}dW_{t},\quad t\leq T \tag{1.12}\]

where \(W\) is a standard Brownian motion. As mentioned before we deal with the spot price whose price is given by a straightforward application of Ito formula

\[S_{t}=F_{0,t}\cdot\exp\big{(}\sigma X_{t}-\frac{1}{2}\lambda_{t}^{2}\big{)}, \quad X_{t}=\int_{0}^{t}e^{-\alpha(t-s)}\,\mathrm{d}W_{s}\ \ \ \text{and}\ \ \lambda_{t}^{2}=\frac{\sigma^{2}}{2\alpha}\big{(}1-e^{-2\alpha t}\big{)}. \tag{1.13}\]

Along with this model, we will use two settings (presented below). In each case we set \(\alpha=4,\sigma=0.7\).

**Case 1**

\[31\text{ exercise dates}\quad\ q_{\min}=0\quad q_{\max}=6\quad Q_{\min}=140\quad Q _{\max}=200.\]

**Case 2** (as in [2, 3])

\[365\text{ exercise dates}\quad\ q_{\min}=0\quad q_{\max}=6\quad Q_{\min}=1300 \quad Q_{\max}=1900.\]

The computer that has been used has the following characteristics: _Processor: Intel(R) Core(TM) i7-1185G7 @ 3.00GHz 3.00 GHz, 32Go of RAM, Microsoft Windows 10 Enterprise_. Deep learning part had been implemented using PyTorch toolbox. The GPU device used is: Nvidia A100-PCIE 40GB.

It is worth noting that for practitioners, the prices of derivatives are closely linked to their sensitivities with respect to market data. These sensitivities are essential for hedging purposes, as they indicate how the price of a derivative product changes with the market. Computation of sensitivities involves calculating derivatives, and in our case, we need to differentiate a price with respect to some parameters, where the price is a solution to the stochastic control problem in equation (1.10). To achieve this, we rely on the so called "envelope theorem."

Let \(f(x,\alpha)\) and \(g_{j}(x,\alpha),j=1,2,\ldots,m\) be real-valued continuously differentiable functions on \(\mathbb{R}^{n+\ell}\), where \(x\in\mathbb{R}^{n}\) are some variables and \(\alpha\in\mathbb{R}^{\ell}\) are parameters, and consider the constrained optimization problem

\[\left\{\begin{array}{l}\max\limits_{x}f(x,\alpha)\\ \text{subject to}\ \ g_{j}(x,\alpha)\geq 0,\quad 1\leq j\leq m\end{array}\right.\]

We introduce the Lagrangian function,

\[\mathcal{L}\left(x,\lambda,\alpha\right)=f(x,\alpha)+\langle\lambda,g(x,\alpha)\rangle\]

where \(\lambda\in\mathbb{R}^{m}\) is the Lagrange multipliers, \(\langle\cdot,\cdot\rangle\) is the Euclidean inner-product and \(g=(g_{1},\ldots,g_{m})^{\top}\). Then we define the value function \(V(\alpha)=f(x^{*}(\alpha),\alpha)\) where \(x^{*}(\alpha)\) is a solution that maximizes the function \(f(\cdot,\alpha)\). The following theorem gives the derivative of the value function \(V\) in case it is differentiable.

**Theorem 1** (Envelope theorem): _Assume that \(V\) and \(\mathcal{L}\) are continuously differentiable. Then,_

\[\frac{\partial V(\alpha)}{\partial\alpha_{k}}=\frac{\partial\mathcal{L}\left( x^{*}(\alpha),\lambda^{*}(\alpha),\alpha\right)}{\partial\alpha_{k}}\quad\quad k =1,\ldots,\ell\]

_where \(\frac{\partial\mathcal{L}}{\partial\alpha_{k}}=\frac{\partial f}{\partial \alpha_{k}}+\langle\lambda,\frac{\partial g}{\partial\alpha_{k}}\rangle\)._
The envelope theorem states that, under some regularity conditions, in a SOC problem, for example our problem (1.10), the derivative of the optimal objective function \(\mathcal{J}\) is given by differentiating the cash flows along with the optimal control. In our model (1.12), we deduce the following proposition which is a corollary application of Theorem 1.

**Proposition 1**: _Let \(\left(q_{k}^{*}\right)_{0\leq k\leq n-1}\) be a solution of the problem (1.10). Assume that the function_

\[\left(F_{0,t_{k}}\right)_{0\leq k\leq n-1}\mapsto\mathbb{E}\left(\sum_{k=0}^{n- 1}q_{k}^{*}\cdot\left(F_{0,t_{k}}\cdot e^{\sigma X_{t_{k}}-\frac{1}{2}\Lambda_{ t_{k}}}-K\right)\right)\]

_is continuously differentiable. Let \(P_{0}\) be the price of the swing contract (1.10). Note that \(P_{0}\) is a function of \(\left(S_{t_{k}}\right)_{0\leq k\leq n-1}\) which in turns depends on \(\left(F_{0,t_{k}}\right)_{0\leq k\leq n-1}\) in the model (1.12). Then for all \(k=0,\ldots,n-1\), the delta (sensitivity of swing price with respect to the initial forward price) is given by_

\[\frac{\partial P_{0}}{\partial F_{0,t_{k}}}=\mathbb{E}\left(q_{k}^{*}\cdot e^ {\sigma X_{t_{k}}-\frac{1}{2}\Lambda_{t_{k}}}\right).\]

Proof We define the following functions:

\[f\left((q_{k})_{0\leq k\leq n-1},(F_{0,t_{k}})_{0\leq k\leq n-1} \right):=\mathbb{E}\left(\sum_{k=0}^{n-1}q_{k}\cdot\left(F_{0,t_{k}}\cdot e^{ \sigma X_{t_{k}}-\frac{1}{2}\Lambda_{t_{k}}}-K\right)\right)\] \[g_{1}\left((q_{k})_{0\leq k\leq n-1},(F_{0,t_{k}})_{0\leq k\leq n -1}\right):=\sum_{k=0}^{n-1}q_{k}-Q_{\min}\] \[g_{2}\left((q_{k})_{0\leq k\leq n-1},(F_{0,t_{k}})_{0\leq k\leq n -1}\right):=Q_{\max}-\sum_{k=0}^{n-1}q_{k}\]

and for all \(k=0,\ldots,n-1\):

\[g_{2k+3}\left((q_{k})_{0\leq k\leq n-1},(F_{0,t_{k}})_{0\leq k\leq n-1} \right):=q_{k}-q_{\min}\]

\[g_{2k+4}\left((q_{k})_{0\leq k\leq n-1},(F_{0,t_{k}})_{0\leq k\leq n-1} \right):=q_{\max}-q_{k}\]

Then it suffices to prove that functions \(f,g_{1},g_{2},\ldots,g_{2n+2}\) are continuously differentiable in order to use the envelope theorem. This holds for functions \(g_{1},g_{2},\ldots,g_{2n+2}\). It remains to prove it for function \(f\). The latter reduces to show that function \(f\) is continuously differentiable in each of its components. For any \(k=0,\ldots,n-1\), the random variable \(q_{k}\cdot\left(F_{0,t_{k}}\cdot e^{\sigma X_{t_{k}}-\frac{1}{2}\Lambda_{t_{k} }}-K\right)\) is integrable since

\[\left|q_{k}\cdot\left(F_{0,t_{k}}\cdot e^{\sigma X_{t_{k}}-\frac{1}{2}\Lambda_ {t_{k}}}-K\right)\right|\leq q_{\max}\cdot\left(F_{0,t_{k}}\cdot e^{\sigma X_ {t_{k}}-\frac{1}{2}\Lambda_{t_{k}}}+K\right)\in\mathbb{L}_{\mathbb{R}}^{1}( \mathbb{P}).\]

Moreover, the function \(F_{0,t_{k}}\mapsto q_{k}\cdot\left(F_{0,t_{k}}\cdot e^{\sigma X_{t_{k}}-\frac {1}{2}\Lambda_{t_{k}}}-K\right)\) is differentiable and its derivative does not depends on \(F_{0,t_{k}}\). Furthermore,

\[\left|\frac{\partial}{\partial F_{0,t_{k}}}\Big{(}q_{k}\cdot\left(F_{0,t_{k}} \cdot e^{\sigma X_{t_{k}}-\frac{1}{2}\Lambda_{t_{k}}}-K\right)\Big{)}\right|= \left|q_{k}\cdot e^{\sigma X_{t_{k}}-\frac{1}{2}\Lambda_{t_{k}}}\right|\leq q _{\max}\cdot e^{\sigma X_{t_{k}}}\in\mathbb{L}_{\mathbb{R}}^{1}(\mathbb{P}).\]
Then thanks to Lebesgue theorem to interchange derivation and integral, the function \(f\) is continuously differentiable in all \(F_{0,t_{k}}\). Likewise one may also show that for all \(k=0,\ldots,n-1\) the function \(f\) is continuously differentiable in \(q_{k}\), and one may use the envelope theorem. For any \(k=0,\ldots,n-1\) we introduce the Lagrangian of the problem (1.10)

\[\mathcal{L}\left(\left(q_{k}\right)_{0\leq k\leq n-1},\lambda, \left(F_{0,t_{k}}\right)_{0\leq k\leq n-1}\right) =\mathbb{E}\left(\sum_{k=0}^{n-1}q_{k}\cdot\left(F_{0,t_{k}}\cdot e ^{\sigma X_{t_{k}}-\frac{1}{2}\Lambda_{t_{k}}}-K\right)\right)+\] \[\langle\lambda,g\big{(}(q_{k})_{0\leq k\leq n-1},(F_{0,t_{k}})_{0 \leq k\leq n-1}\big{)}\rangle,\]

where \(g=\left(g_{1},\ldots,g_{2n+2}\right)^{\top}\) and \(\lambda\in\mathbb{R}^{2n+2}\). Then it follows from envelope theorem that for any \(k=0,\ldots,n-1\):

\[\frac{\partial P_{0}}{\partial F_{0,t_{k}}}=\mathbb{E}\left(\sum_{k=0}^{n-1}q _{k}^{*}\cdot\frac{\partial}{\partial F_{0,t_{k}}}\left(F_{0,t_{k}}\cdot e^{ \sigma X_{t_{k}}-\frac{1}{2}\Lambda_{t_{k}}}-K\right)\right)=\mathbb{E}\left( \sum_{k=0}^{n-1}q_{k}^{*}\cdot e^{\sigma X_{t_{k}}-\frac{1}{2}\Lambda_{t_{k}}} \right).\]

This completes the proof.

We have outlined the theoretical framework for pricing swing contracts and computing their associated sensitivities. To obtain practical numerical solutions, we introduce parametric methods in this paper. Specifically, we replace the optimal control variable \(q_{k}\) at a time \(t_{k}\) with a parametric function \(q_{k}\big{(}I_{k};\theta\big{)}\), where \(I_{k}\) represents information needed to determine the volume to purchase, and \(\theta\) is a parameter belonging to a parameter space \(\Theta\) that depends on the chosen parameterization. Thus, the optimization problem in equation (1.10) can be expressed as,

\[\sup_{\theta\in\Theta}\mathcal{J}\Big{(}q_{0}(I_{0};\theta),\ldots,q_{n-1}(I_ {n-1};\theta)\Big{)}. \tag{1.14}\]

Note that the parameter space has to be designed to guarantee that, almost surely, for \(\theta\in\Theta\), \(\big{(}q_{0}(I_{0};\theta),\ldots,q_{n-1}(I_{n-1};\theta)\big{)}\) is an admissible strategy. That is, it lies within \(\mathcal{A}_{0,0}^{Q_{\min},Q_{\max}}\) which is defined through (1.8) and where the range of admissible volumes at each time \(t_{k}\) depends on purchased volumes up to that time.

## 2 Swing pricing: A global optimization approach

In this paper, we propose two parametric methods for approximating the optimal control (i.e., the volume to purchase at each exercise date) in a swing pricing context. It is worth noting that the use of parametric exercise strategies for pricing swing contracts is not a new approach and has demonstrated its advantages when compared to classical methods such as the Longstaff and Schwartz method (see [3, 27]). For instance, Gobet et al. [3] developed two parametric methods, in the context of swing contracts with penalties, which provide satisfactory prices when compared to their benchmark (obtained using the forest of trees method [24]). The first method is based on a neural network that approximates a \([0,1]\)-valued function \(f\), which is then used to determine the volume to purchase in the range \([q_{\min},q_{\max}]\) (recall that they studied a swing contract with penalties. Thus the range of admissible volume is not constrained as in our case through amounts \(A^{-}\) and \(A^{+}\) (1.6)) using the simple transformation \(q_{\min}+(q_{\max}-q_{\min})f\). In addition to the previous parameterization, the authors proposed another one based on the heuristic that a
higher gas price usually results in a higher consumption at a fixed strike price (though this may not always be true due to global constraints). This heuristic suggests the existence of a threshold beyond which the highest possible volume should be purchased. We studied this heuristic using a neural network to replace the volume \(q_{k}\) with a nonlinear parameterization \(q_{k}:=q(t_{k},S_{t_{k}}-K;\theta)\), where \(\theta\) are the parameters of the neural network trained to solve the swing pricing problem (1.14) (see section 3 for details on the training process). The results showed that when the payoff \(S_{t_{k}}-K\) is below a certain threshold \(\beta_{1}^{k}\), the optimal consumption tends to be \(A_{k}^{-}\) (as defined in equation (1.6)), while for payoffs above a certain threshold \(\beta_{2}^{k}(>\beta_{1}^{k})\), the optimal consumption tends to be \(A_{k}^{+}\) (also defined in equation (1.6)). Therefore, the optimal consumption profile exhibits a similar behavior as that depicted in Figure 2.

Based on this observation, we aim to find a parametric function that reproduces the consumption profile shown in Figure 2, which depends on two thresholds \(\beta_{1}\) and \(\beta_{2}\) (that are assumed to be constant over time in a first step). To achieve this, we define a parametric "decision" function

\[I_{k}\mapsto f_{k}(I_{k};\theta):=f(I_{k};\theta)=\mathbf{1}_{\{I_{k}>\beta_{ 2}\}}+\frac{I_{k}-\beta_{1}}{\beta_{2}-\beta_{1}}\mathbf{1}_{\{\beta_{1}\leq I _{k}\leq\beta_{2}\}} \tag{2.1}\]

where \(\mathbf{1}\) denotes the indicator function, \(\theta=(\beta_{1},\beta_{2})^{\top}\in\mathbb{R}^{2}\) and \(I_{k}=S_{t_{k}}-K\) is the payoff. Note that the thresholds and therefore the parameterization do not depend on time. The preceding "decision" function, once defined, allows to choose an admissible control using the following transformation

\[q_{k}(I_{k};\theta)=A_{k}^{-}(Q_{k}^{\theta})+\left(A_{k}^{+}(Q_{k}^{\theta}) -A_{k}^{-}(Q_{k}^{\theta})\right)\cdot f_{k}(I_{k};\theta) \tag{2.2}\]

with \(Q_{0}^{\theta}=0\) and \(Q_{k}^{\theta}=\sum_{i=0}^{k-1}q_{i}(I_{i};\theta)\) for all \(1\leq k\leq n-1\). Then we aim at solving the parametric optimization problem (1.14) with the preceding choice of parameterization (2.2) and with \(\Theta=\{(x,y)\in\mathbb{R}^{2};x<y\}\). To achieve this, we look for values of the thresholds \(\theta=(\beta_{1},\beta_{2})\) which maximize the resulting reward function \(\mathcal{J}\) (using strategy given by (2.1)). To implement this approach, we simulate \(10^{5}\) independent realizations of the payoff \(I_{k}=S_{t_{k}}-K\) and compute the Monte-Carlo price given by strategy (2.2). The values of \(\beta_{1},\beta_{2}\) used to compute this maximum are the \(50\) evenly distributed within the interval \([-10,10]\), for the contract setting of case \(2\) (as described in section 1.3).

The optimal values for \(\beta_{1}\) and \(\beta_{2}\) that lead to the maximum price of \(2611.7\) are obtained as \(-2.24\) and \(-1.02\), respectively. This simple parametric strategy which has no flexibility (the strategy is the same for all exercise dates since thresholds are constants) and does not depend on the cumulative volume (which is obviously essential when deciding which volume to buy at each exercise dates) gives a satisfying

Figure 2: Optimal consumption profile as a function of the payoff.
price when compared to the more elaborated methods in [3]1. Following the promising results obtained by this straightforward and simple strategy, we introduce two improvements leading to two different parameterizations.

Footnote 1: Note the relative error between this price and the best of theirs is roughly \(2\%\) and they studied the penalty case whose theoretical price is greater than that of the firm constraints case we study in this paper.

* **Explicit Payoff-Volume parameterization**: we suppose the optimal consumption behaves like profile 2 and define a (smooth) parametric function reproducing this profile. The final idea is to use stochastic gradient descent based optimization algorithms to find the best parameters. In this parametrization, we allow the strategy to depend not only on time but also on cumulative consumption. This latter is useful and gives flexibility to our strategy.
* **Neural network parameterization**: this is a variant of the preceding parameterization. Here we replace coefficients which define the first parameterization by the output of a neural network. Once the neural network is trained, we get the parameters and define the same strategy.

To illustrate both points of view, we use the settings of case 1 in the following sections.

### Explicit Payoff-Volume parameterization (_PV strat_)

We aim to find a parameterization that replicates the target profile in Figure 2, and that is also sufficiently regular to allow for optimization using Stochastic Gradient Descent methods. The parameterization needs to be a sufficiently regular function to ensure that the resulting reward function \(\mathcal{J}\) is differentiable. To achieve this, we select the logistic function, denoted by \(\sigma(x):=\frac{1}{1+e^{-x}}\), which is infinitely differentiable on \(\mathbb{R}\) and has a shape close enough to that of the target profile. Additionally, to incorporate flexibility into our exercise strategy, we no longer assume that the parameters are constant over exercise dates. We also need to take into account for the dependence of the strategy on cumulative consumption. Hence, we propose the following parameterization.

\[f_{k}(I_{k};\theta)=\sigma\big{(}\langle\theta_{k},I_{k}\rangle\big{)}\qquad \text{with}\quad I_{k}=\big{(}S_{t_{k}}-K,M(Q_{k}),1\big{)}^{\top}\in\mathbb{R }^{3} \tag{2.3}\]

where for some vector \(\theta=(\theta_{1},\ldots,\theta_{n})\in\Theta=\mathbb{R}^{3n}\), \(\theta_{k}\in\mathbb{R}^{3}\) (\(1\leq k\leq n\)) is the subvector of \(\theta\) made of the \(k^{th},(k+1)^{th},(k+2)^{th}\) components of \(\theta\). Thus vector \(\theta\) embed all parameters used to define the strategy for all exercise dates. In other words, each \(\mathbb{R}^{3}\)-valued component \(\theta_{k}\) represents coefficients that drive the optimal decision of _PV strat_ at the corresponding time \(t_{k}\). The function \(M\) is chosen as the (normalized) margin or remaining purchasing capacity defined by: \(M(Q):=\frac{Q-Q_{\min}}{Q_{\max}-Q_{\min}}\). The importance of this choice is discussed in Remark 1 and 2. In case \(Q_{\min}=Q_{\max}\), one may use \(M(Q)=\frac{Q-Q_{\min}}{Q_{\min}}\). Then the problem to solve is the same as in (1.14) where the parameters space is \(\Theta=\mathbb{R}^{3n}\) and the control \(q_{k}(I_{k};\theta)\) in (2.2) is defined using the strategy (2.3).

**Remark 1** (Importance of function _M_): At first glance, the function \(M\) may seem to be a simple normalization of the cumulative consumption. However, it should be noted that several other normalizations were tested and did not lead to a good strategy. We tested the following normalizations: \(M(Q_{k})=Q_{k}\), \(\frac{Q_{k}-Q_{\min}}{q_{\max}-q_{\min}}\), \(\frac{k\cdot q_{\max}-Q_{k}}{q_{\max}-q_{\min}}\) and \(\frac{Q_{k}-Q^{down}(t_{k+1})}{Q^{up}(t_{k+1})-Q^{down}(t_{k+1})}\).

**Remark 2** (Remaining capacity _Versus_ Cumulative Consumption): We observed that the remaining purchasing capacity is a more crucial factor than the current cumulative consumption when determining
the volume to buy at each exercise date. This contrasts with the neural network approach tested in [3], which uses the current cumulative consumption. When using \(M(Q)=Q\), our algorithm tended to purchase the maximum possible volume at the contract's outset, as long as the payoff was positive. This behavior suggests that the algorithm did not take into account the fact that purchasing capacity diminishes as the contract approaches its expiration, due to global constraints.

Remarks 1 and 2 explain the reason for which function \(M\) yields desirable results. In addition, this choice has a rational interpretation for the optimal control behavior (see appendix D).

### Neural Network parameterization (_NN strat_)

Our second parameterization is based on neural networks. The goal of a neural network is to approximate some mapping \(x\mapsto\Phi^{*}(x)\) by a parametric one (neural network) \(x\mapsto\Phi(x;\theta)\) where \(\theta\) is a parameter (or weights) of the neural network that has to be optimized such in order to give a "good" approximation. A Neural Network can approximate a wide class of complicated (linear or non-linear) relationships (For more details, see Universal Approximation Theorem [18]) between the inputs and outputs by composing linear functions and non-linear thresholds. More precisely, note that a neural network is made of nodes connected to one another where a column of nodes forms a layer (when there are more than one layer in the neural network architecture, we speak of deep neural network). The outermost (see diagram 3) are the input and output layers and all those in between are called the hidden layers. The connection between the input and output layers through hidden layer is made by means of linear functions and activation functions (non-linear ones):

* **Linear functions**: a node in connected to another with an associated weight (\(w\)) and bias (\(b\)) so that it will receive a weighted value (\(wx+b\)) from a node in the previous layer. \(\theta=(w,b)\) forms parameters of the neural network.
* **Non-linear functions**: each node contains an activation function that activates the node (so the input can pass through the node) when the input of the node is above a certain threshold.

Training a neural network is performed through two steps: **forward propagation** and **backward propagation**. The forward propagation leads the predicted output \(\hat{y}\) starting from the input \(x\) layer-by-layer. More formally, the forward propagation comes down to the evaluation of the following function:

\[x\in\mathbb{R}^{d}\mapsto\Phi(x;\theta):=a_{I}^{\theta_{I}}\circ\phi_{q_{I-1}} \circ a_{I-1}^{\theta_{I-1}}\circ\ldots\circ\phi_{q_{1}}\circ a_{1}^{\theta_{1 }}(x)\in\mathbb{R}^{\ell} \tag{2.4}\]

where

\(\rhd I,q_{I},...,q_{1}\) are positive integers specifying the depth of the network and the number of nodes for each hidden layer.

\(\rhd a_{1}^{\theta_{1}}:\mathbb{R}^{d}\rightarrow\mathbb{R}^{q_{1}},\ldots,a_{ I-1}^{\theta_{I-1}}:\mathbb{R}^{q_{I-2}}\rightarrow\mathbb{R}^{q_{I-1}}\) and \(a_{I}^{\theta_{I}}:\mathbb{R}^{q_{I-1}}\rightarrow\mathbb{R}^{\ell}\) are affine functions. As explained above, theses functions are of the form: \(a(x)=Wx+b\) where \(W\) is a matrix of weights and \(b\) is a vector of bias.

\(\rhd\) For \(j\in\mathbb{N}\), \(\phi_{j}:\mathbb{R}^{j}\rightarrow\mathbb{R}^{j}\) are the activation functions. In this paper we only use ReLU (Rectified Linear Unit) function i.e, \(x\mapsto\max(x,0)\).

In our case, the neural network is designed/trained to provide a "decision" function which maximizes the reward function \(\mathcal{J}\). To achieve this, we design a \(\mathbb{R}^{3}\)-valued neural network \(\Phi(\cdot;\theta)\) (thus \(\ell=3\), see
diagram 3) of form (2.4) and which depends on the time \(t_{k}\), the payoff \(S_{t_{k}}-K\) and the margin \(M(Q_{k})\) defined as in (2.3). Thus \(d=3\) and parameter \(\theta\) of the neural network \(\Phi(\cdot;\theta)\) lies within

\[\Theta=\mathbb{R}^{q_{1}}\times\mathbb{R}^{3\times q_{1}}\times \Big{(}\prod_{i=2}^{I-1}\mathbb{R}^{q_{i}}\times\mathbb{R}^{q_{i}\times q_{i-1 }}\Big{)}\times\mathbb{R}^{3}\times\mathbb{R}^{3\times q_{I-1}}. \tag{2.5}\]

As for _PV strat_, we define recursively the "decision" function by \(f_{0}(I_{0};\theta)=\sigma\Big{(}\langle\Phi\big{(}t_{0},S_{t_{0}}-K,M(0); \theta\big{)},I_{0}\rangle\Big{)}\) and for all \(1\leq k\leq n-1\)

\[f_{k}(I_{k};\theta)=\sigma\Big{(}\langle\Phi\big{(}t_{k},S_{t_{k}}-K,M(Q_{k}^{ \theta});\theta\big{)},I_{k}\rangle\Big{)}\quad\text{ with }\,Q_{k}^{\theta}=\sum_{i=0}^{k-1}q_{i}(I_{i};\theta) \tag{2.6}\]

where information \(I_{k}\) is the same as in (2.3).Then the parametric optimization problem to solve is the same as defined in equation (1.14) with the preceding choice of \(\Theta\) given by equation (2.5) and where the control is defined using the strategy (2.6) combined with (2.2).

**Remark 3** (_PV strat Versus NN strat_): Note that in _NN strat_, the parameter \(\theta\in\Theta\) (as defined in (2.5)) replaces all vectors \(\theta_{k}\in\mathbb{R}^{3}\) in equation (2.3), which drives the strategy in _PV strat_. Unlike the latter, in _NN strat_, the same set of parameters defines the strategy for all exercise dates, but the strategy is different according to the exercise date since time is integrated as input. In contrast, the size of the parameter vector in _PV strat_ increases linearly with the number of exercise dates. Additionally, in section 3.4, we will, numerically, demonstrate that even for a high number of exercise dates, a relatively small neural network architecture is still sufficient. This is one of the advantages of _NN strat_ over _PV strat_.

Our neural network parameterization is a more robust one compared to the one performed by Gobet et Al. [3]. Indeed, in their article they had chosen, for the decision function, a non linear \([0,1]\)-valued parameterization. Our neural network parameterization is an improvement of theirs in the sense that we first impose a particular parametric shape. Then the parameters that drive the latter shape are chosen as

Figure 3: Illustration of (deep) neural network architecture.
an output of a neural network. This approach not only provides flexibility to the strategy but also helps the algorithm to "learn" the optimal control by imposing a specific shape.

All the strategies presented above require training to find the best parameters. To achieve this, we rely on stochastic approximation theory.

## 3 Stochastic optimization

In this section we aim at optimizing a \(\mathbb{R}\)-valued function \(h:y\in\mathbb{R}^{q}\mapsto\mathbb{E}(H(y,Z))\) with \(H\) being a \(\mathbb{R}\)-valued function defined on \(\mathbb{R}^{q}\times\mathbb{R}^{d}\) and where \(Z\) is a \(d-\)dimensional random vector. Like a probabilistic extension of the classic Newton-Raphson procedure, stochastic approximation [23, 34]) suggests the following procedure to find zeros of function \(h\):

\[\forall n\in\mathbb{N},\ \ \ y_{n+1}=y_{n}-\gamma_{n+1}\cdot h(y_{n}) \ \ \,0<\gamma_{n}\leq\gamma_{0}. \tag{3.1}\]

When \(h\) is the gradient of another function i.e, \(h=\nabla V\), we speak about **stochastic gradient descent** (SGD) which is an extension of gradient descent (see [37]) to stochastic optimization. Hereafter we consider a SGD setting i.e, \(h=\nabla V\). The sequence \((\gamma_{n})_{n\in\mathbb{N}}\) is called the algorithm step or learning rate and has to be chosen carefully. In fact if the step is too large, then the procedure (3.1) can overshoot and completely miss the global optimum of the function \(V\). If the step is too small, the procedure may take a long time to converge as it will be taking small steps towards the global minimum. Thus to ensure convergence (which is widely discussed in the literature [5, 10, 35]), a classical requirement is the "decreasing step" assumption. This means that the step sequence is non-increasing and satisfy the following conditions

\[\sum_{n\geq 0}\gamma_{n}=+\infty\ \ \text{and}\ \ \sum_{n\geq 0} \gamma_{n}^{2}<+\infty. \tag{3.2}\]

From a practical point of view, to implement iteration (3.1), one uses mini-batch procedure. That is, rather than computing gradient over a single observation, leading to one iteration in procedure (3.1), we use small batches. We divide a sample of size \(M\) into \(B\) samples of size \(L:=M/B\). Then we perform \(B\) iterations using

\[y_{n+1}=y_{n}-\frac{\gamma_{n+1}}{L}\sum_{\ell=1}^{L}H(y_{n},Z_{n+1}^{b,[\ell ]}),\ \ 0\leq n\leq N. \tag{3.3}\]

where \(Z_{n+1}^{b,[\ell]}\) are independent copies of \(Z\). This alternative increases update frequency compared to vanilla SGD which allows a more robust convergence, avoiding local minima/maxima. It should be noted that in the SGD framework, \(H\) often represents the gradient of a multivariate function, which makes it difficult or impossible to implement the SGD procedure by hand. To address this issue, we use a technique called "Adjoint Algorithm Differentiation" (AAD) as described in [4, 32]. This technique involves a program that computes a value and automatically computes derivatives of that value by combining the derivatives of several simple arithmetic expressions.

Vanilla SGD can often perform poorly due to various reasons. One of the well-known reasons is its inability to converge or a slow convergence to a local minima/maxima when the objective function has saddle points. A saddle point is a critical point on the surface of the function's graph that is not a local extremum. Saddle points are often found in regions that look like plateaus, as shown in Figure 4.

It should be noted that saddle points are notorious in high-dimensional spaces, especially in cases involving deep neural networks, which often have non-convex objective functions (see [13, 26]). In fact, as the
dimensionality of the problem increases, the number of saddle points has been shown to increase exponentially compared to the number of local minima/maxima. For both _PV strat_ and _NN strat_, the number of parameters to optimize is high, especially in _NN strat_. Our experiments have shown that vanilla SGD does not perform well in these high-dimensional spaces. Therefore, in this paper, we propose two alternative optimization algorithms, Adaptive Moment Estimation (Adam) and Stochastic Langevin Dynamics (SGLD), and compare their performance.

### Adaptive Moment Estimation (Adam)

Adam [22] is an efficient stochastic optimization method that only requires first-order gradients with little memory requirement. The method computes individual adaptive learning rates for different parameters from estimates of first gradient moments. Practically, this is done by means of a preconditioning matrix \(P\) (see updating (3.4)). Adam can be seen at as a combination of RMSprop 2 and Stochastic Gradient Descent with momentum [33]. The procedure in Adam updating is the following

Footnote 2: Unpublished optimization algorithm designed for neural networks, first proposed by Geoffrey Hinton in a Coursera course [https://www.cs.toronto.edu/~tijmen/csc321/slides/lecture_slides_lec6.pdf](https://www.cs.toronto.edu/~tijmen/csc321/slides/lecture_slides_lec6.pdf)

\[y_{n+1}=y_{n}-\gamma_{n+1}P_{n+1}\cdot\widehat{\nabla V}(y_{n}) \tag{3.4}\]

where \(\widehat{\nabla V}\) is a (scaled) estimation of the gradient (see Algorithm 1) computed as a moving average of the estimated gradient in the spirit of momentum technique. \((P_{n})_{n}\) is a sequence of preconditioning matrix. In Adam updating the latter matrix is a diagonal matrix whose diagonal components depend on the squared (elementwise) of the gradient \(g_{n+1}\) estimated after iteration \(n\). More precisely, given a gradient \(g=\big{(}g_{1},\ldots,g_{q}\big{)}\in\mathbb{R}^{q}\), the square elementwise of this gradient denoted \(g\odot g\) is given by

\[g\odot g=\big{(}g_{1}^{2},\ldots,g_{q}^{2}\big{)}. \tag{3.5}\]

In the same spirit, for some matrix \(A=\big{(}a_{i,j}\big{)}_{1\leq i,j\leq q},B=\big{(}b_{i,j}\big{)}_{1\leq i,j \leq q}\) where \(b_{i,j}\neq 0\) for all \(1\leq i,j\leq q\) we define elementwise division as follows

\[A\oslash B=\big{(}c_{i,j}\big{)}_{1\leq i,j\leq q}\quad\text{ where }\;c_{i,j}=a_{i,j}/b_{i,j}. \tag{3.6}\]

Adam algorithm is described in 1.

Figure 4: Saddle points illustration. \((x,y)\mapsto x^{3}-3xy^{2}\) (on left) and \((x,y)\mapsto x^{2}-y^{2}\) (on right). The point \((0,0)\) is a saddle point.
**Algorithm 1** Adam updating. \(\lambda\) is a small correction term to avoid division by zero. \(\odot\) denotes element-wise multiplication and \(\oslash\) denotes elementwise division defined in (3.5) and (3.6). \(Id_{q}\) represents the \(q\times q\) identity matrix.

**Require:**

\(\rhd\)\(\gamma\): Step size

\(\rhd\)\(\mu_{1},\mu_{2}\in[0,1)\): Exponential decay rates for the moment estimates

\(\rhd\)\(V(y)\): Stochastic objective function with parameters \(y\). Note that \(V\) is computed by its Monte-Carlo version using (3.3)

\(\rhd\)\(y_{0}\): Initial parameter vector

\(M_{0}=0\) (Initialize \(1^{st}\) moment vector)

\(MS_{0}=0\) (Initialize \(2^{nd}\) moment vector)

```
1:while\(y_{n}\) not converged do
2:\(g_{n+1}\leftarrow\nabla_{y}V(y_{n})\)\(\triangleright\) Compute gradient w.r.t to parameters
3:\(M_{n+1}\leftarrow\mu_{1}\cdot M_{n}+(1-\mu_{1})\cdot g_{n+1}\)\(\triangleright\) Update biased first moment estimate
4:\(MS_{n+1}\leftarrow\mu_{2}\cdot MS_{n}+(1-\mu_{2})\cdot g_{n+1}\odot g_{n+1}\)\(\triangleright\) Update biased second raw moment estimate
5:\(\widehat{M}_{n+1}\gets M_{n+1}/(1-\mu_{1}^{n+1})\)\(\triangleright\) Compute bias-corrected first moment estimate
6:\(\widehat{MS}_{n+1}\gets MS_{n+1}/(1-\mu_{1}^{n+1})\)\(\triangleright\) Compute bias-corrected second raw moment estimate
7:\(P_{n+1}=\operatorname{diag}\left(Id_{q}\oslash(\lambda\cdot Id_{q}+\sqrt{ \widehat{MS}_{n+1}})\right)\)\(\triangleright\) Compute preconditioning matrix
8:\(y_{n+1}\gets y_{n}-\gamma\cdot P_{n+1}\cdot\widehat{M}_{n+1}\)\(\triangleright\) Update parameters
9:endwhile
```

**Algorithm 2** Adam updating. \(\lambda\) is a small correction term to avoid division by zero. \(\odot\) denotes elementwise multiplication and \(\oslash\) denotes elementwise division defined in (3.5) and (3.6). \(Id_{q}\) represents the \(q\times q\) identity matrix.

Adam uses the squared gradients to scale the learning rate like RMSprop and it takes advantage of momentum by using the moving average (driven by \(\mu_{1},\mu_{2}\)) of the gradient instead of the gradient itself like SGD with momentum. The scaling of the learning rate by the squared gradient allows to solve gradient magnitude problem. Otherwise, sometimes the gradients may be huge and sometimes small which intricate the choice of a learning rate. The exponential moving average of the gradient allows to "de-noise" the estimation of the gradient. Indeed, it is recommended to choose \(\mu_{1},\mu_{2}\approx 1\) so that the contribution of the current estimation of the gradient is small since we give a higher weight to previous data points and lower one to current data points. This allows to smooth subsequent estimation of gradient.

Adam has become the default algorithm for stochastic optimization and especially for neural networks training. Its theoretical convergence had been widely studied in the literature. First, the authors of the algorithm analyzed the convergence in a convex framework with bounded gradient assumption. For non-convex cases, one may refers to [14, 40].

In this paper, we do not limit our study to Adam updating. We propose to explore other algorithm based on Langevin dynamics and which shares, as Adam, good properties when dealing with saddle points.

### Preconditioned Stochastic Gradient Langevin Dynamics (PSGLD)

Stochastic Gradient Langevin Dynamics (SGLD) has been introduced for bayesian learning [39] to approximate a posterior distribution given some data. SGLD combines Robbins-Monro type algorithm with Langevin dynamics which injects noise into the parameter updates such that the trajectory of the parameters converges to the posterior distribution. To be more precise, it has been shown that, under mild conditions, the latter posterior distribution is the unique invariant probability measure of the Langevin
Stochastic Differential Equation (SDE):

\[dy_{t}=-\nabla V(y_{t})dt+\sqrt{2}dW_{t}, \tag{3.7}\]

where \(W\) is a \(q\)-dimensional Brownian motion. Practically, the posterior distribution is approximated using an Euler discretization of the Langevin SDE

\[y_{n+1}=y_{n}-\gamma\nabla V(y_{n})+\sqrt{2\gamma}Z_{n+1}, \tag{3.8}\]

where \((Z_{n})_{n}\) is a sequence of i.i.d. standard \(q\)-dimensional Gaussian vectors.

From equations (3.1) (with \(h=\nabla V\)) and (3.8), SGLD appears to be an extension of vanilla SGD by adding an exogenous white noise to the gradient descent. This allows to regularize the problem and escape from traps (see [28]). The use of Langevin dynamics based algorithms for optimization is justified by the fact that recent studies (see [11, 12]) showed that sampling from a distribution which concentrates around the global minimum of \(V\) is a similar task as minimizing \(V\) via certain optimization algorithms. But, under some conditions, it can be shown (see [29]) that if \((Y_{t}^{y})_{t\geq 0}\) is a \(\mathbb{R}^{d}\)-valued random vector solution of the following SDE

\[dY_{t}=b(Y_{t})dt+\sigma(Y_{t})dW_{t},\ \ Y_{0}=y \tag{3.9}\]

with \(W\) a standard \(q\)-dimensional Brownian motion and a drift term of the form

\[b:=-\frac{1}{2}\Big{(}(\sigma\sigma^{\top})\nabla V-\Big{[}\sum_{j=1}^{d} \partial_{y_{j}}(\sigma\sigma^{\top})_{ij}\Big{]}_{i=1:d}\Big{)}. \tag{3.10}\]

Then, the distribution

\[\nu_{V}(dy)=C_{V}e^{-V(y)}\cdot\lambda_{d}(dy) \tag{3.11}\]

is the unique invariant distribution of the SDE (3.9). In particular if \(\sigma=\sqrt{2}\cdot Id_{d}\) then get the langevin SDE (3.7). Which shows that the latter SDE converges to its stationary distribution, namely the Gibbs measure \(\propto\exp(-V(y))\) which, in fact, concentrates on the global minimum of \(V\). All this motivates the use of Langevin dynamics based algorithms for stochastic optimization (see [6, 8, 21]). In this paper, as in [8], we consider the application of Langevin algorithms in a non-bayesian setting and we will consider the preconditioned version of SGLD which leads to Preconditioned SGLD (PSGLD, see [25]). The latter choice is due to the fact that in its standard version, SGLD (without noise) acts like SGD by updating all parameters with the same step size leading to slow mixing when the components of \(y\) have different curvature. Given a sequence of preconditioning matrix \((P_{n})_{n}\), the updating procedure (3.8), reads in PSGLD setting

\[y_{n+1}=y_{n}-\gamma_{n+1}P_{n+1}\cdot\nabla V(y_{n})+\sigma_{n+1}\sqrt{ \gamma_{n+1}}\mathcal{N}(0,P_{n+1}) \tag{3.12}\]

where \((\sigma_{n})_{n}\) is a constant or non-increasing sequence controlling the amount of injected noise. As in Adam updating the preconditioning matrix sequence \((P_{n})_{n}\) is a diagonal matrix whose components depend on the squared (elementwise) of the gradient estimation. The general procedure is presented in Algorithm 2.
The convergence of Langevin algorithm has been widely studied. One may refers to [15, 16] when the sequence \((\sigma_{n})_{n}\) is constant and to [7, 29] when the sequence is decreasing.

In the next section we implement _PV strat_ and _NN strat_ with settings of case 1 and optimize them using both Adam and PSGLD presented above.

### Adam versus PSGLD

We implement _PV strat_ and _NN strat_ using either Adam or PSGLD optimization algorithm with case 1 setting. For this case the Longstaff-Schwartz method, with canonical polynomial functions of degree 3 provides a price of 65.14 with a 95 % confidence interval: \([65.08,65.21]\) and in 4 minutes and 40 second in average. Note that with case 1 (even for case 2) setting, regardless the chosen methods, the resulting Monte-Carlo estimator of the swing price exhibits high variance due to the high variance of the underlying. Thus, in the following table and just for this section, prices are computed by averaging _100 replications_ of prices. Each of the 100 prices is obtained using \(B=1\) batch of size \(L=M=2^{14}\) and a learning rate \(\gamma=0.1\). In what follows we denote by \(N\) the number of iterations in the stochastic procedure (see (3.3)), by \(\widehat{P_{0}}\) the (raw) price resulting from _PV strat_ or _NN strat_. And \(\widehat{P_{0}}^{BG}\) denotes the price given by forcing the decision function (either (2.3) for _PV strat_ or (2.6) for _NN strat_) to be strictly bang-bang. Results are recorded in Tables 1 and 2 for _PV strat_ and in Tables 2 and 4 for _NN strat_.
Tables 1 and 2 suggest that Adam requires several iterations and, as a result, increases computation time when using a learning rate of 0.1. Whereas with only 1000 iterations, using PSGLD optimization algorithm, leads to a better price. This fact is also observable in _NN strat_ where PSGLD optimization algorithm allows to achieve a better price than Adam.

Tables 3 and 3 show that prices obtained by _NN strat_ are better than that of _PV strat_ and, as already mentioned, PSGLD optmization algorithm provides a better exercise strategy. This means that adding some noise in the optimization procedure helps to converge quickly. Following the performance of PSGLD updating, we consider this algorithm in the remainder along with the following configuration: \(\sigma=1\cdot e^{-6},\beta=0.8,\lambda=1\cdot e^{-10}\). With the latter setting, one may compute sensitivity of the swing price with respect to the initial forward price (see Figure 5).

\begin{table}
\begin{tabular}{c c c c c} \hline \hline \(\sigma\) & \(\beta\) & \(\widehat{P_{0}}\) & \(\widehat{P_{0}}^{BG}\) & Time (s) \\ \hline \hline \(1\cdot e^{-6}\) & 0.9 & 65.23 ([65.16, 65.30]) & 65.21 ([65.14, 65.28]) & 53.4 \\ \hline \(1\cdot e^{-5}\) & 0.9 & 65.23 ([65.14, 65.31]) & 65.22 ([65.13, 65.30]) & 52.3 \\ \hline \(1\cdot e^{-6}\) & 0.8 & 65.27 ([65.20, 65.35]) & 65.26 ([65.18, 65.33]) & 51.3 \\ \hline \end{tabular}
\end{table}
Table 4: Results for _NN strat_ using PSGLD optimization algorithm. For the neural network architecture, we used 2 hidden layers (\(I=2\)) with 10 units per layer (\(q_{1}=10,q_{2}=10\)). Values in brackets are confidence intervals (95%). Columns “time” includes the training and the valuation times. We used \(N=1000\) iterations and \(\lambda=1\cdot e^{-10}\).

\begin{table}
\begin{tabular}{c c c c c} \hline \hline \(\sigma\) & \(\beta\) & \(\widehat{P_{0}}\) & \(\widehat{P_{0}}^{BG}\) & Time (s) \\ \hline \hline \(1\cdot e^{-6}\) & 0.8 & 65.19 ([65.12, 65.26]) & 65.21 ([65.13, 65.28]) & 21.3 \\ \hline \(1\cdot e^{-6}\) & 0.9 & 65.21 ([65.15, 65.28]) & 65.23 ([65.17, 65.29]) & 21.3 \\ \hline \(1\cdot e^{-5}\) & 0.9 & 65.22 ([65.16, 65.28]) & 65.23 ([65.18, 65.29]) & 21.6 \\ \hline \end{tabular}
\end{table}
Table 2: Results for _PV strat_ using PSGLD optimization algorithm. Values in brackets are confidence intervals (95%). Column “time” includes both training and valuation times. We used \(N=1000\) iterations and \(\lambda=1\cdot e^{-10}\).

\begin{table}
\begin{tabular}{c c c c c} \hline \hline \(N\) & \(\widehat{P_{0}}\) & \(\widehat{P_{0}}^{BG}\) & Time (s) \\ \hline \hline
3000 & 65.22 ([65.16, 65.28]) & 65.23 ([65.17, 65.29]) & 150.4 \\ \hline
1000 & 65.13 ([65.05, 65.20]) & 65.22 ([65.15, 65.29]) & 53.2 \\ \hline \end{tabular}
\end{table}
Table 3: Results for _NN strat_ using Adam optimization algorithm. For the neural network architecture, we used 2 hidden layers (\(I=2\)) with 10 units per layer (\(q_{1}=10,q_{2}=10\)). Values in brackets are confidence intervals (95%). Columns “time” includes the training and the valuation times.
Note that for both optimization algorithms (Adam and PSGLD) different other hyperparameters had been tested and results are reported in Appendix E.

### Practitioner's corner: Transfer learning

As mentioned, in practice, training both parameterizations we proposed in this paper is allowed by AAD. However, the latter can be time consuming especially in the case of a swing contract having several exercise dates and using more than one batch per iteration (see Appendix E). For instance, let us consider a swing contract with maturity one year and daily exercise (365 exercise dates). Alongside we use the diffusion model (1.12) with the setting of case 2. The results are reported in the table below. As observed in [3] the Longstaff-Schwartz method generates numerical instabilities. Moreover, note that we get better prices than optimal quantization method like in [2].

From Table 5 one may notice that the computation time is not negligible even if, compared to other methods (like Longstaff-Schwartz) it remains reasonable. To reduce it, we propose a method based on the so called **transfer learning** (for details see [31]). Transfer learning refers to a machine learning method where a model developed for a task is reused as the starting point for a model on a second and similar task. This method may accelerates the evaluation of swing contract with several exercise dates. In our case, we perform "transfer learning" as follows. For a swing contract with several exercise dates, we first consider an aggregated version. That means that we consider a contract with less exercise dates and where the local constraints are aggregated. For instance for a swing contract over one year with daily exercise and with \(q_{\min}=0,q_{\max}=6\), we may consider a contract with one exercise per month (the middle of each month) and where for example, for months with 30 days, the local constraints become \(q_{\min}=0\times 30=0,q_{\max}=30\times 6=180\). Then we run few iterations to optimize the pricing of the

Figure 5: Delta forward for setting of case 1 (left) and case 2 (right) using _PV strat_ and _NN strat_.

\begin{table}
\begin{tabular}{l l l l} \hline \hline  & \(\widehat{P_{0}}\) & \(\widehat{P_{0}}^{BG}\) & \((T_{train},T_{eval})\) \\ \hline \hline _PV strat_ & 2690.10 ([2689.20, 2691.00]) & 2692.59 ([2691.68, 2693.49]) & (274.1, 74.7) \\ \hline _NN strat_ & 2693.74 ([2692.83, 2694.64]) & 2694.12 ([2693.21, 2695.02]) & (680.9, 398.3) \\ \hline \hline \end{tabular}
\end{table}
Table 5: Results for a one-year swing contract. Values in brackets are confidence intervals (95%). The valuation was performed with a sample of size \(1\cdot e^{8}\). For the training we used \(N=1000\) iterations. \(T_{train}\) denotes the training time and \(T_{eval}\) the valuation time.
aggregated contract. The resulting parameters are used as an initial guess for the actual pricing problem. To extend the parameters obtained in the problem with 12 exercise dates to the problem with 365 exercise dates, we assume that all days within the same month behave the same. Therefore the initial guess parameters for all days within the same month are those obtained in the aggregated problem for this month. It turns out that this allows to reach more quickly a descent optimum.

We can notice from Table 6 that by means of transfer learning the computation time is drastically reduced. With our adaptation of transfer learning, the training time is reduced by a factor of 3 without degrading the accuracy. Besides, always with in mind the aim of reducing the computation time, transfer learning can be used differently. In a few words, by means of transfer learning, when the market data change or when the contracts settings change we do not need to start all over again. To illustrate this point let us consider below three cases: M1, M2, M3, representing some shifts of market/contract settings. Recall that the baseline settings are those consider above (case 2 setting).

In Table 7 we consider three different market move scenarios described as follows. In the first case we bumped the initial forward curve from 20 to 22. In the second case we change the global constraints from \(Q_{\min}=1300\) and \(Q_{\max}=1900\) to \(Q_{\min}=1400\) and \(Q_{\max}=2000\). In the final case we reduced the strike price from 20 to 18. We could have changed the local constraints but as we already pointed out, swing pricing can always reduce to the case where \(q_{\min}=0,q_{\min}=1\) (see Appendix A). Now the question is the following: Could our baseline model give accurate prices without training again? Or could it help us to get accurate prices quickly? The answers to theses questions are recorded in the following tables.

\begin{table}
\begin{tabular}{l l l l} \hline \hline  & \(\widehat{P_{0}}\) & \(T_{agg}\) & \((T_{train},T_{eval})\) \\  & \(\widehat{P_{0}}^{BG}\) & & \\ \hline \hline _PV strat_ & 2688.73 ([2687.82, 2689.63]) & 4.8 & (83.2, 76.3) \\  & 2692.81 ([2691.90, 2693.71]) & & \\ \hline _NN strat_ & 2693.22 ([2692.31, 2694.13]) & 6.9 & (206.5, 398.1) \\  & 2694.12 ([2693.21, 2695.02]) & & \\ \hline \end{tabular}
\end{table}
Table 6: Results for a one-year swing contract using transfer learning. We used 500 iterations for the aggregated contract and 300 iterations for the actual one-year contract. For the valuation, we used a sample of size \(1\cdot e^{8}\). \(T_{agg}\) denotes the training time for the aggregated contract, \(T_{train}\) and \(T_{eval}\) denote respectively the training time and the valuation time for the actual contract.

\begin{table}
\begin{tabular}{l l l} \hline \hline
**Case M1** & **Case M2** & **Case M3** \\ \hline \hline \(F_{0,t_{k}}=22\) & \((Q_{\min},Q_{\max})=(1400,2000)\) & \(K=18\) \\ \hline \end{tabular}
\end{table}
Table 7: Market data move scenarios
Tables 8 and 9 demonstrate that in the case where global constraints change (M2), using the baseline parameters without additional training gives prices that are comparable to those obtained by training with a few iterations starting from the baseline parameters. Furthermore, it is important to note that increasing the number of iterations beyond 300 (we tested up to 1000 iterations in our experiments) does not lead to better prices than those obtained after 300 iterations (see columns labeled "Re-train"). This suggests that reusing the baseline parameters can speed up convergence. This observation remains true regardless of the case (M1, M2, or M3).

## 4 Numerical experiments

We now perform additional simulations to demonstrate the effectiveness of the two parameterizations that we proposed in this paper.

### Three factor model

We first consider a three factor model whose dynamics is given by

\begin{table}
\begin{tabular}{l l l l} \hline \hline Cases & Re-use & Re-train & \((T_{train},T_{eval})\) \\ \hline \hline
**Case M1** & 5942.65 ([5941.64, 5943.65]) & 6005.75 ([6004.70, 6006.79]) & (85.4, 77.2) \\  & 5944.27 ([5943.26, 5945.27]) & 6008.54 ([6007.50, 6009.58]) & \\ \hline
**Case M2** & 2514.86 ([2513.90, 2515.81]) & 2516.99 ([2516.04, 2517.93]) & (85.0, 76.9) \\  & 2517.26 ([2516.30, 2518.21]) & 2518.73 ([2517.78, 2519.67]) & \\ \hline
**Case M3** & 5655.98 ([5655.06, 5656.89]) & 5744.91 ([5743.95, 5745.86]) & (85.9, 77.05) \\  & 5657.67 ([5656.75, 5658.58]) & 5747.39 ([5746.43, 5748.34]) & \\ \hline \end{tabular}
\end{table}
Table 8: Results for _PV strat._ Column “Re-use” provides results when the baseline model parameters are reused as is. Column “Re-train” gives results when we re-train with only 300 iterations with the baseline model parameters used as starting values. \(T_{train}\) denotes the training time and \(T_{eval}\) the valuation time. The testing computation time is the same when we re-use the model as when we re-train.

\begin{table}
\begin{tabular}{l l l l} \hline \hline Cases & Re-use & Re-train & \((T_{train},T_{eval})\) \\ \hline \hline
**Case M1** & 5981.17 ([5980.14, 5982.19]) & 6011.08 ([6010.03, 6012.12]) & (203.1, 397.8) \\  & 5982.57 ([5981.54, 5983.59]) & 6011.62 ([6010.57, 6012.66]) & \\ \hline
**Case M2** & 2509.97 ([2509.01, 2510.90]) & 2516.92 ([2515.96, 2517.87]) & (203.7, 397.5) \\  & 2511.08 ([2510.12, 2512.03]) & 2517.58 ([2516.62, 2518.53]) & \\ \hline
**Case M3** & 5713.80 ([5712.85, 5714.72]) & 5749.97 ([5749.01, 5750.92]) & (203.4, 397.6) \\  & 5715.24 ([5714.34, 5716.17]) & 5750.51 ([5749.55, 5751.46]) & \\ \hline \end{tabular}
\end{table}
Table 9: Results for _NN strat_. Column “Re-use” provides results when the baseline model parameters are reused as is. Column “Re-train” gives results when we re-train with only 300 iterations with the baseline model parameters used as starting values. \(T_{train}\) denotes the training time and \(T_{eval}\) the valuation time. The testing computation time is the same when we re-use the model as when we re-train.
\[\frac{dF_{t,T}}{F_{t,\,T}}=\sigma_{1}e^{-\alpha_{1}(T-t)}dW_{t}^{1}+\sigma_{2}e^{- \alpha_{2}(T-t)}dW_{t}^{2}+\sigma_{3}e^{-\alpha_{3}(T-t)}dW_{t}^{3}, \tag{4.1}\]

where for all \(1\leq i,j\leq 3\), the instantaneous correlation is given by

\[\langle dW_{\cdot}^{i},dW_{\cdot}^{j}\rangle_{t}=\left\{\begin{array}{ll}dt& \text{ if }\,i=j\\ \rho_{i,j}\cdot dt&\text{ if }\,i\neq j\end{array}\right.\]

In this model, the spot price is given by

\[S_{t}=F_{0,t}\cdot\exp\Big{(}\langle\sigma,X_{t}\rangle-\frac{1}{2}\lambda_{t }^{2}\Big{)},\]

where \(\sigma=\big{(}\sigma_{1},\sigma_{2},\sigma_{3}\big{)}^{\top}\), \(X_{t}=\big{(}X_{t}^{1},X_{t}^{2},X_{t}^{3}\big{)}^{\top}\) and for all \(1\leq i\leq 3\),

\[X_{t}^{i}=\int_{0}^{t}e^{-\alpha_{i}(t-s)}\,\mathrm{d}W_{s}^{i}\ \ \text{and}\ \,\lambda_{t}^{2}=\frac{1}{2}\sum_{i=1}^{3}\frac{\sigma_{i}^{2}}{\alpha_{i}} \big{(}1-e^{-2\alpha_{i}t}\big{)}+\sum_{i\neq j}\rho_{i,j}\frac{\sigma_{i} \sigma_{j}}{\alpha_{i}+\alpha_{j}}\Big{(}1-e^{-(\alpha_{i}+\alpha_{j})t}\Big{)}.\]

We fix the following configuration \(\sigma_{i}=\sigma=0.7,\alpha_{i}=\alpha=1.5,\rho_{i,j}=\rho\in[-1,1]\). The swing contract setting corresponds to case 1. We use \(N=1000\) iterations and a learning rate \(\gamma=0.1\). Results are recorded in Tables 10 and 11.

It should be noted that there is a different approach to implement _NN strat_, and taking into account the structure of forward prices. Specifically, in the three-factor framework (4.1), the forward price depends on state variables that are components of the \(\mathbb{R}^{3}\)-valued random vector \(X_{t}\). In this context, one can include

\begin{table}
\begin{tabular}{c l l} \hline \hline \(\rho\) & \(\widehat{P_{0}}\) & \(\widehat{P_{0}}^{BG}\) \\ \hline \hline
0.6 & 172.71 ([172.51, 172.90]) & 172.72 ([172.52, 172.92]) \\ \hline
0.3 & 147.79 ([147.63, 147.95]) & 147.78 ([147.62, 147.95]) \\ \hline -0.2 & 91.02 ([90.92, 91.12]) & 91.02 ([90.92, 91.12]) \\ \hline \end{tabular}
\end{table}
Table 10: Results using _PV strat_. Values in brackets are confidence intervals (95%). The valuation had been performed with a sample of size \(1\cdot e^{8}\). For each result, the execution time (training plus testing) is roughly equal to 22s.

\begin{table}
\begin{tabular}{c l l} \hline \hline \(\rho\) & \(\widehat{P_{0}}\) & \(\widehat{P_{0}}^{BG}\) \\ \hline \hline
0.6 & 173.05 ([172.85, 173.24]) & 172.98 ([178.78, 173.17]) \\ \hline
0.3 & 147.97 ([147.81, 148.14]) & 147.92 ([147.75, 148.08]) \\ \hline -0.2 & 91.18 ([91.08, 91.28]) & 91.13 ([91.03, 91.23]) \\ \hline \end{tabular}
\end{table}
Table 11: Results using _NN strat_. Values in brackets are confidence intervals (95%). The valuation had been performed with a sample of size \(1\cdot e^{8}\). For the neural network architecture we used \(I=2\) layers with \(q_{1}=q_{2}=10\) units. For each result, the execution time (training plus testing) is roughly equal to 45s.
as an additional input to the neural network described in (2.4), in addition to time \(t_{k}\), payoff \(S_{t_{k}}-K\), and volume normalization \(M(Q_{k})\). That is, using the above notation, \(I_{k}=\left(t_{k},S_{t_{k}}-K,X_{t_{k}},M(Q_{k})\right)\in\mathbb{R}^{d+3}\), where \(d\) is the dimension of \(X_{t_{k}}\) (in the three-factor model, \(d=3\)). This approach can help to capture the correlation structure in a multi-factor framework. The prices resulting from this alternative approach are presented in Table 12.

It appears that, regardless of whether the model is multi-factor or not, both strategies (_PV strat_ and _NN strat_) can determine the optimal volume to purchase at each exercise date, based only on the payoff and the cumulative volume. In the following section, we will evaluate the performance of our strategies on a more complex diffusion model.

### Multi-curve forward diffusion

We finally consider a multi-curve model. At a certain valuation date \(t\) we consider \(p\) risk factors. Each risk factor \(i\in\{1,\ldots,p\}\), denoted by \(F_{t,T_{i}}\), is a forward contract observed at date \(t\) and expiring at date \(T_{i}\). It is modeled through an instantaneous volatility function \(\sigma_{t}(T_{i})\) with a dynamics given by,

\[\frac{dF_{t,T_{i}}}{F_{t,T_{i}}}=\sigma_{t}(T_{i})dW_{t}(T_{i}) \tag{4.2}\]

where \((W_{t}(T_{i}))_{t\geq 0}\) is a standard Brownian motion. The instantaneous correlation between Brownian motions are given by

\[\langle dW_{\cdot}(T_{i}),dW_{\cdot}(T_{j})\rangle_{t}=\rho_{i,j}dt\]

The model (4.2) implies to diffuse one curve per risk factor; \(F_{t,T_{i}}\) defined by its maturity date \(T_{i}\). Likewise, the spot at a given date generates its own risk factor. Thus if we consider a swing contract on the day-ahead where we can exercise on all working days of a year, then we have about 365 delivery dates. Therefore 365 curves to diffuse. But due to high correlation (at any date, two risk factors are very correlated when maturities are very close), one can see that only a few risk factors really impact the price. In practice forward contracts delivering in the same month are highly correlated. Thus we will consider there are as many factor as months. For example, for a contract on a whole year \(p=12\); therefore 12 (correlated) Brownian motions per exercise date. In this model we consider a valuation date given by \(17^{th}\) march 2021 and a swing contract delivering on the period January 2022-December 2022. The swing constraints are: \(q_{\min}=0,q_{\max}=1,Q_{\min}=180,Q_{\max}=270\). For the diffusion model, the correlation matrix is (_in percentage_)

\begin{table}
\begin{tabular}{c l l} \hline \hline \(\rho\) & \(\widehat{P_{0}}\) & \(\widehat{P_{0}}^{BG}\) \\ \hline \hline
0.6 & 172.76 ([172.56, 172.95]) & 172.70 ([172.50, 172.89]) \\ \hline
0.3 & 147.95 ([147.78, 148.11]) & 147.91 ([147.75, 148.07]) \\ \hline -0.2 & 91.11 ([91.01, 91.20]) & 91.07 ([90.97, 91.17]) \\ \hline \end{tabular}
\end{table}
Table 12: Results using _NN strat_ including state variables. Values in brackets are confidence intervals (95%). The valuation had been performed with a sample of size \(1\cdot e^{8}\). For the neural network architecture we used \(I=2\) layers with \(q_{1}=q_{2}=10\) units. For each result, the execution time (training plus testing) is roughly equal to 50s.
Volatilities and initial forward prices are assumed to be constant per month. For months in 2022, the square of the instantaneous volatility \(\sigma_{t}^{2}(T_{i})\) is given by Figure 6 and the initial forward curve is given by Figure 7.

Finally, the strike price corresponds to the average of all initial forward prices; which gives 17.865. For this diffusion model, _PV strat_ does not perform well. This is probably due to the multi-curve framework

Figure 6: Square of instantaneous volatility per risk factor observed on 17-March-2021. Values are: \(\left(45.09\%,44.88\%,42.63\%,40.57\%,36.99\%,34.37\%,31.18\%,29.59\%, 30.22\%,30.6\%,31.05\%,30.25\%\right)\)

Figure 7: Forward prices per risk factor observed on 17-March-2021. Prices per risk factor are: \(\left(20.07,20,19.6,17.4,16.75,16.5,16.56,16.53,16.71,17.31,18.31,18.64\right)\)
which is more complex than a multi factor model. However modifying _NN strat_ as in the preceding section by adding, in neural network inputs, Brownian motions \(\big{(}W_{t_{k}}(T_{i})\big{)}_{1\leq i\leq p}\) at each exercise date \(t_{k}\) allows to achieve very good prices. We compare the latter implementation of _NN strat_ with a method based on a static replication of swing contracts by means of spread options. The latter method aims at approximating swing options prices by selecting a combination of spread options (_SO_) where the underlying is driven by the same dynamics (4.2). The selection is done through a linear programming method. Results are recorded in Table 13.

## Conclusion

We introduced two parametric approaches for pricing swing contracts. The first method involves a direct parameterization of the optimal control based on some heuristics, while the second method improves upon the first by using neural networks. We conducted numerical experiments to compare two optimization algorithms (Adam and PSGLD). Our results demonstrate that using Langevin-based optimization algorithms allows us to achieve better prices with short computation time. We also found that our methods outperform the state-of-the-art methods in terms of accuracy. Additionally, we tested our neural network parameterization in various complex diffusion models and demonstrating its robustness. We can hence conclude that our neural network-based method is well-suited to the pricing swing options.

## Acknowledgments

The first author is grateful for the financial support provided by Engie Global Markets via CIFRE agreement and would like to thank Asma Meziou and Frederic Muller for throughout insights on swing contracts.

## Appendix A Swing contract decomposition

We aim to prove that any swing contract can be reduced to a normalized contract, namely a swing contract with local constraints \(q_{\min}=0,q_{\max}=1\). The swing contract price is given by

\[P_{0}=\operatorname*{ess\,sup}_{(q_{\ell})_{0\leq\ell\leq n-1} \in\mathcal{A}^{\mathcal{Q}_{\min},\,Q_{\max}}_{0,0}}\ \mathbb{E}\left(\sum_{\ell=0}^{n-1}q_{\ell}\times(S_{t_{\ell}}-K)\right)\]

where

\[\mathcal{A}^{\mathcal{Q}_{\min},\,Q_{\max}}_{0,0}=\left\{(q_{ \ell})_{0\leq\ell\leq n-1},\ q_{\ell}:(\Omega,\mathcal{F}^{S}_{t_{\ell}}, \mathbb{P})\mapsto[q_{\min},q_{\max}],\ \sum_{\ell=0}^{n-1}q_{\ell}\in\,[Q_{\min},Q_{\max}]\right\}.\]

It follows from the linearity of the expectation that,

\begin{table}
\begin{tabular}{c c c} \hline \hline \(\widehat{P_{0}}\) & \(\widehat{P_{0}}^{BG}\) & \(SO\) \\ \hline \hline
708.46 ([705.76, 711.16]) & 711.32 ([708.61, 714.03]) & 645.3 \\ \hline \end{tabular}
\end{table}
Table 13: Results using _NN strat_. Values in brackets are confidence intervals (95%). The valuation had been performed with a sample of size \(5\times 1\cdot e^{6}\). For the neural network architecture we used \(I=2\) layers with \(q_{1}=q_{2}=50\) units.
\[P_{0}=q_{\min}\times\mathbb{E}\left(\sum_{\ell=0}^{n-1}(S_{t_{\ell}}-K)\right)+(q_{ \max}-q_{\min})\times\underset{(q_{\ell})_{0\leq\ell\leq n-1}\in\mathcal{A}_{0,0 }^{Q_{\min},\,Q_{\max}}}{\text{ess}\sup}\ \mathbb{E}\left(\sum_{\ell=0}^{n-1}\frac{q_{\ell}-q_{\min}}{q_{\max}-q_{\min} }\times(S_{t_{\ell}}-K)\right).\]

The first term in the last equality is given by

\[\mathbb{E}\left(\sum_{\ell=0}^{n-1}(S_{t_{\ell}}-K)\right)=\sum_{\ell=0}^{n-1} \mathbb{E}\left(S_{t_{\ell}}\right)-n\cdot K\]

and can be easily computed using either closed formula (depending on the underlying diffusion model) or Monte-Carlo method. Let us consider the second term. Let \((q_{\ell})_{0\leq\ell\leq n-1}\in\mathcal{A}_{0,0}^{Q_{\min},Q_{\max}}\) and define \(\tilde{q}_{\ell}=\frac{q_{\ell}-q_{\min}}{q_{\max}-q_{\min}}\). Note that \((\tilde{q}_{\ell})_{0\leq\ell\leq n-1}\in\mathcal{A}_{0,0}^{\tilde{Q}_{\min},\tilde{Q}_{\max}}\) where

\[\mathcal{A}_{0,0}^{\tilde{Q}_{\min},\tilde{Q}_{\max}}=\left\{(q_{\ell})_{0 \leq\ell\leq n-1},\ q_{\ell}:(\Omega,\mathcal{F}_{t_{\ell}},\mathbb{P}) \mapsto[0,1],\ \sum_{\ell=0}^{n-1}q_{\ell}\in[\tilde{Q}_{\min},\tilde{Q}_{\max}]\right\}.\]

and

\[\tilde{Q}_{\min}=\frac{(Q_{\min}-n\cdot q_{\min})_{+}}{q_{\max}-q_{\min}} \qquad\tilde{Q}_{\max}=\frac{(Q_{\max}-n\cdot q_{\min})_{+}}{q_{\max}-q_{\min}}.\]

Thus,

\[\mathbb{E}\left(\sum_{\ell=0}^{n-1}\frac{q_{\ell}-q_{\min}}{q_{ \max}-q_{\min}}\times(S_{t_{\ell}}-K)\right) =\mathbb{E}\left(\sum_{\ell=0}^{n-1}\tilde{q}_{\ell}\times(S_{t_ {\ell}}-K)\right)\] \[\leq\underset{(q_{\ell})_{0\leq\ell\leq n-1}\in\mathcal{A}_{0,0} ^{\tilde{Q}_{\min},\tilde{Q}_{\max}}}{\text{ess}\sup}\ \mathbb{E}\left(\sum_{\ell=0}^{n-1}q_{\ell}\times(S_{t_{\ell}}-K)\right).\]

Therefore taking the supremum yields,

\[\underset{(q_{\ell})_{0\leq\ell\leq n-1}\in\mathcal{A}_{0,0}^{Q_{\min},Q_{ \max}}}{\text{ess}\sup}\ \mathbb{E}\left(\sum_{\ell=0}^{n-1}\frac{q_{\ell}-q_{\min}}{q_{\max}-q_{\min} }\times(S_{t_{\ell}}-K)\right)\leq\underset{(q_{\ell})_{0\leq\ell\leq n-1}\in \mathcal{A}_{0,0}^{Q_{\min},\tilde{Q}_{\max}}}{\text{ess}\sup}\ \mathbb{E}\left(\sum_{\ell=0}^{n-1}q_{\ell}\times(S_{t_{\ell}}-K)\right).\]

Conversely let \((q_{\ell})_{0\leq\ell\leq n-1}\in\mathcal{A}_{0,0}^{\tilde{Q}_{\min},\tilde{Q} _{\max}}\) and define \(\tilde{q}_{\ell}=q_{\min}+(q_{\max}-q_{\min})\cdot q_{\ell}\in[q_{\min},q_{\max}]\). It follows \(\sum_{\ell=0}^{n-1}q_{\ell}\in[Q_{\min},Q_{\max}]\) so that \(\left(\tilde{q}_{\ell}\right)_{0\leq\ell\leq n-1}\in\mathcal{A}_{0,0}^{Q_{ \min},Q_{\max}}\). Thus,

\[\mathbb{E}\left(\sum_{\ell=0}^{n-1}q_{\ell}\times(S_{t_{\ell}}-K)\right) =\mathbb{E}\left(\sum_{\ell=0}^{n-1}\frac{\tilde{q}_{\ell}-q_{\min }}{q_{\max}-q_{\min}}\times(S_{t_{\ell}}-K)\right)\] \[\leq\underset{(q_{\ell})_{0\leq\ell\leq n-1}\in\mathcal{A}_{0,0}^ {Q_{\min},\,Q_{\max}}}{\text{ess}\sup}\ \mathbb{E}\left(\sum_{\ell=0}^{n-1}\frac{q_{\ell}-q_{\min}}{q_{\max}-q_{\min} }\times(S_{t_{\ell}}-K)\right).\]

Taking the supremum, we get,

\[\underset{(q_{\ell})_{0\leq\ell\leq n-1}\in\mathcal{A}_{0,0}^{\tilde{Q}_{\min},\,\tilde{Q}_{\max}}}{\text{ess}\sup}\ \mathbb{E}\left(\sum_{\ell=0}^{n-1}q_{\ell}\times(S_{t_{\ell}}-K)\right)\leq \underset{(q_{\ell})_{0\leq\ell\leq n-1}\in\mathcal{A}_{0,0}^{Q_{\min},\,Q_{ \max}}}{\text{ess}\sup}\ \mathbb{E}\left(\sum_{\ell=0}^{n-1}\frac{q_{\ell}-q_{\min}}{q_{\max}-q_{\min} }\times(S_{t_{\ell}}-K)\right).\]

Therefore,

\[P_{0}=q_{\min}\times\mathbb{E}\left(\sum_{\ell=0}^{n-1}(S_{t_{\ell}}-K)\right)+( q_{\max}-q_{\min})\times\underset{(q_{\ell})_{0\leq\ell\leq n-1}\in\mathcal{A}_{0,0}^ {Q_{\min},\,Q_{\max}}}{\text{ess}\sup}\ \mathbb{E}\left(\sum_{\ell=0}^{n-1}q_{\ell}\times(S_{t_{\ell}}-K)\right).\]
On the rate of convergence

This section aims to estimate, at least numerically, the rate of convergence. We denote by \(U_{n}\) the price obtained after the \(n^{th}\) iteration in the training step. We make the assumption that for some constant \(C>0\)

\[U_{n}=U_{\infty}+\frac{C}{n^{\alpha}}\]

where \(U_{\infty}\) is the limit of the stochastic procedure that we assume to exist. Thus,

\[\log(|U_{2n}-U_{n}|)=K_{\alpha}+\alpha\ \log(\frac{1}{n}),\ \ \ \ \text{ where }K_{\alpha}=\log(|C(2^{-\alpha}-1)|).\]

Therefore the coefficient \(\alpha\) (representing the rate of convergence following the assumption) appears to be the slope in the \(\log-\log\) regression of \(|U_{2n}-U_{n}|\) against \(\frac{1}{n}\).

From figures 8 and 9 it suggests that both parameterizations give a rate of convergence of order \(\mathcal{O}(\frac{1}{N})\) (where \(N\) is the number of iterations). This convergence rate is much faster than Monte-Carlo, and suggests that our parameterizations are efficient alternatives to Longstaff-Schwartz method for this problem.

## Appendix C Estimator variance and computation time

In this section we illustrate the high variance phenomenon which may appear when pricing swing option. To this end, we compute 100 realizations of swing contract price with case 1 setting for the three methods: _PV strat_, _NN strat_ and Longstaff-Schwartz. The distributions of prices are represented in Figures 10, 11, 12.
Figures 10, 11, and 12 demonstrate that prices can significantly fluctuate, regardless of the pricing method. Using \(1000000\) simulations is insufficient to obtain a reliable price estimate. Storage limitations prevent the Longstaff-Schwartz method from exceeding this number of simulations, and this limit is even lower when the number of exercise dates is high. However, our proposed methods allow to increase the number of simulations as needed once the strategy trained. This is possible because we can evaluate our strategies on mini-batches sequentially instead of using the entire test set. For instance, to compute a price with a sample size of \(10^{8}\), we estimate \(100\) prices over \(100\) mini-batches, each with a sample size of \(10^{6}\). The final price is then the average of the \(100\) estimated prices.

Hereafter (see Figure 13) we present CPU time for _PV strat_ and _NN strat_. We use one batch with size of \(2^{14}\) and the same PSGLD setting.

Figure 11: Distribution of swing prices using _NN strat_. From left to right we used successively \(10^{5},10^{6},10^{7},5\times 10^{8}\) simulations.

Figure 12: Distribution of swing prices using Longstaff-Schwartz method. From left to right we used successively \(10^{5},10^{6}\) simulations. Higher number of simulations leads to memory overflow.

Figure 10: Distribution of swing prices using _PV strat_. From left to right we used successively \(10^{5},10^{6},10^{7},5\times 10^{8}\) simulations.
## Appendix D _PV strat_ coefficients

In this appendix we present the coefficients obtained with _PV strat_ and either Adam or PSGLD updating. For that purpose, we used \(N=1000\) iterations, a learning rate of \(0.1\) and the following configurations. For Adam, we use \(B=4\) batches of size \(2^{12}\). For the PSGLD updating, we use \(\sigma=1\cdot e^{-6},\beta=0.8,\lambda=1\cdot e^{-10}\) and \(B=1\) batch of size \(2^{14}\). In the following graphics, coefficients \(a_{k},b_{k}\) represent the coefficients which multiply respectively the payoff \(S_{t_{k}}-K\), the margin of cumulative consumption \(M(Q_{k})\). The coefficient \(c_{k}\) is the constant (see (2.3)). For both optimization algorithms, the strategy is estimated with the same sample.

In general, it can be observed that the coefficients estimated using the PSGLD algorithm are smoother across the exercise dates compared to those estimated with the Adam updating method. Additionally, _PV strat_ results in interpretable coefficients. Specifically, the coefficient \(a_{k}\) tends to increase globally, with a rapid growth over the last exercise dates. This is because, as we move forward in exercise dates, the minimal global constraint is gradually fulfilled, giving us more flexibility on the control or allowed volumes to purchase. In this particular example, the maximum global constraint is always fulfilled, so towards the last exercise dates, we can decide to buy the minimum possible amount if the payoff is negative and buy the maximum possible amount if the payoff is positive (i.e., the parameterization (2.3) approaches 0 or 1).

On the other hand, the coefficient \(b_{k}\) generally decreases at the beginning of the contract and starts increasing towards zero over the last exercise dates. This indicates that the remaining exercise capacity has less influence on the strategy towards the end of the contract. This is consistent with the claim made earlier, which stated that towards the end of the contract, the optimal strategy depends mainly on the payoff and less on the volume constraints.

Figure 14: Coefficients of _PV strat_ using PSGLD updating.

Figure 13: CPU time (in seconds) as a function of the number of iterations.
## Appendix E Summary tables: Adam versus PSGLD

Hereafter are recorded all the results obtained using Adam and PSGLD updatings with different hyperparameters combinations. Prices are computed by averaging 100 replications of prices; each obtained using a sample size of \(10^{6}\). We consider a swing contract with case 1 setting. Results are recorded in tables 14, 15, 16, 17.
\begin{table}
\begin{tabular}{c c c c c} \hline \hline \(\sigma\) & \(\beta\) & \(\widehat{P_{0}}\) & \(\widehat{P_{0}}^{BG}\) & Time (s) \\ \hline \hline \(1e^{-5}\) & 0.8 & 65.15 ([65.09, 65.22]) & 65.17 ([65.10, 65.23]) & 22.6 \\ \hline \(1e^{-5}\) & 0.7 & 65.13 ([65.05, 65.21]) & 65.15 ([65.07, 65.22]) & 22.1 \\ \hline \(1e^{-5}\) & 0.9 & 65.22 ([65.16, 65.28]) & 65.23 ([65.18, 65.29]) & 21.6 \\ \hline \(1e^{-6}\) & 0.8 & 65.19 ([65.12, 65.26]) & 65.21 ([65.13, 65.28]) & 21.3 \\ \hline \(1e^{-6}\) & 0.9 & 65.21 ([65.15, 65.28]) & 65.23 ([65.17, 65.29]) & 21.4 \\ \hline \end{tabular}
\end{table}
Table 16: Summary table for _PV strat_ using PSGLD. The values in brackets are confidence intervals (95%). The column “time” includes both training and valuation time. We used a learning rate equal to 0.1, one batch of size \(2^{14}\), \(\lambda=1\cdot e^{-10}\) and 1000 iterations.

\begin{table}
\begin{tabular}{c c c c c} \hline \hline \(L\times B\) & \(\gamma\) & \(N\) & \(\widehat{P_{0}}\) & Time (s) \\ \hline \hline \(2^{14}\times 1\) & 0.1 & 3000 & 65.22 ([65.16, 65.28]) & 65.23 ([65.17, 65.29]) & 150.4 \\ \hline \(2^{14}\times 1\) & 0.1 & 1000 & 65.13 ([65.05, 65.20]) & 65.22 ([65.15, 65.29]) & 53.2 \\ \hline \(2^{12}\times 4\) & 0.1 & 3000 & 65.26 ([65.21, 65.32]) & 65.23 ([65.17, 65.29]) & 603.4 \\ \hline \(2^{12}\times 4\) & 0.1 & 1000 & 65.20 ([65.14, 65.26]) & 65.20 ([65.14, 65.26]) & 133.5 \\ \hline \(2^{14}\times 1\) & 0.01 & 3000 & 64.78 ([62.65, 64.90]) & 65.06 ([64.98, 65.14]) & 157.1 \\ \hline \(2^{14}\times 1\) & 0.01 & 1000 & 64.81 ([64.74, 64.88]) & 65.05 ([64.98, 65.12]) & 55.1 \\ \hline \(2^{12}\times 4\) & 0.01 & 3000 & 65.20 ([65.15, 65.26]) & 65.20 ([65.14, 65.26]) & 607.3 \\ \hline \(2^{12}\times 4\) & 0.01 & 1000 & 65.02 ([64.92, 65.11]) & 65.17 ([65.09, 65.24]) & 135.4 \\ \hline \end{tabular}
\end{table}
Table 15: Summary table for _NN strat_ using Adam. We used a neural network architecture as follows: 2 hidden layers (\(I=2\)) and 10 units per layer (\(q_{1}=10,q_{2}=10\)). The values in brackets are confidence intervals (95%). The time includes the training and the valuation time.

\begin{table}
\begin{tabular}{c c c c c} \hline \hline \(\sigma\) & \(\beta\) & \(\widehat{P_{0}}\) & \(\widehat{P_{0}}^{BG}\) & Time (s) \\ \hline \hline \(1e^{-5}\) & 0.8 & 65.15 ([65.09, 65.22]) & 65.17 ([65.10, 65.23]) & 22.6 \\ \hline \(1e^{-5}\) & 0.7 & 65.13 ([65.05, 65.21]) & 65.15 ([65.07, 65.22]) & 22.1 \\ \hline \(1e^{-5}\) & 0.9 & 65.22 ([65.16, 65.28]) & 65.23 ([65.18, 65.29]) & 21.6 \\ \hline \(1e^{-6}\) & 0.8 & 65.19 ([65.12, 65.26]) & 65.21 ([65.13, 65.28]) & 21.3 \\ \hline \(1e^{-6}\) & 0.9 & 65.21 ([65.15, 65.28]) & 65.23 ([65.17, 65.29]) & 21.4 \\ \hline \end{tabular}
\end{table}
Table 16: Summary table for _PV strat_ using PSGLD. The values in brackets are confidence intervals (95%). The column “time” includes both training and valuation time. We used a learning rate equal to 0.1, one batch of size \(2^{14}\), \(\lambda=1\cdot e^{-10}\) and 1000 iterations.