A Multilevel Monte Carlo Method for the Valuation of Swing Options

In this study, we propose a novel approach for the valuation of swing options. Swing options are a kind of American options with multiple exercise rights traded in energy markets. Longstaff and Schwartz have suggested a regression-based Monte Carlo method known as the least-squares Monte Carlo (LSMC) method to value American options. In this work, first we introduce the LSMC method for the pricing of swing options. ,en, to achieve a desired accuracy for the price estimation, we combine the idea of LSMC with multilevel Monte Carlo (MLMC) method. Finally, to illustrate the proper behavior of this combination, we conduct numerical results based on the Black–Scholes model. Numerical results illustrate the efficiency of the proposed approach.


Introduction
Swing options are American type contracts with multiple early exercise rights in energy markets, such as gas, oil, and electricity. In these markets, unknown future demands due to weather conditions, political situations, and unstorability of energy lead to price fluctuations. Swing contracts give the option holders multiple exercise rights in long or short positions up to the maturity time to protect themselves from fluctuations in commodity prices. A swing option gives its holder the right to change the volume of periodic deliveries of a particular commodity on certain dates in the future, within a given delivery period, at the prescribed strike price [1].
American option has no closed-form pricing formula, and consequently swing options have no explicit solutions; thus, numerical methods should be employed to price these options approximately. Numerical methods for pricing the swing options can be classified into the following three categories: lattice/tree methods, numerical methods for partial differential equations (PDEs), and Monte Carlo (MC) simulations. ompson [2] was the first to present a paper on the pricing of swing options as take-or-pay options, in which he proposed a binomial tree algorithm to price a version of swing options. Later, Lari-Lavassani et al. developed another binomial tree in which the spot price was a two-factor mean-reverting process [3]. Jaillet et al. applied multi-layered trinomial trees based on a dynamic programming approach for valuing swing options [4]. Indeed, they proposed a one-factor mean-reverting process for energy prices that explicitly incorporates seasonal effects. In the study by Shao and Xiang [5], the pricing problem was formulated as an optimal stochastic control problem, which was solved by the trinomial forest approach while the underlying asset followed a seasonal mean-reverting regimeswitching model. Moreover, Kjaer in [6] implemented a finite difference method to solve PIDEs arising from a dynamic programming problem associated with the pricing of swing options. It is worth mentioning that due to easy implementation of PDE-based methods, they have become popular in financial problems in the last decades [6,7].
Carmona and Touzi viewed swing options as optimal multiple stopping problems and solved the problem using the MC method and Malliavin calculus [8]. ey used the theory of the Snell envelope to determine the optimal exercise boundaries. Zeghal and Mnif in [9] extended the latter method to Lévy processes for valuing swing options. Ben Latifa et al. in [10] valued swing options with the jump in price process. ey transformed optimal multiple stopping time problem into a sequence of ordinary stopping time problems and solved them as the unique viscosity solution of Hamilton-Jacobi-Bellman variational inequality. Moreover, Naouara and Trabelsi [11] proposed a generalization of swing options to the case of an optimal multiple stopping time problem with random number of exercise rights. ey transformed the main problem into a sequence of standard optimal stopping time problems using dynamic programming.
e least-squares Monte Carlo (LSMC) method is a regression-based Monte Carlo simulation for valuing early exercise options. For the first time, Carriere [12] proposed the LSMC method, which was later used by Longstaff and Schwartz [13] for approximating the value of American options. ese studies applied regression to approximate value functions in their backwards dynamic program. Figueroa [14] used an extension of the LSMC method for valuing swing options, as described in [15]. In these papers, regression was used to estimate continuation value (CV), and decisions for exercising rights were made based on the estimated CV. e pricing of American-type options via the MC method is a challenging task, as it requires the backward dynamic programming algorithm, which seems to be incompatible with the forward structure of MC methods [16]. Belomestny et al. in [16] implemented local and global regression methods to estimate the CV. ey combined this idea with multilevel Monte Carlo (MLMC) method to price Bermudan options to reduce the complexity of the LSMC method.
In this paper, we focus on the pricing of swing options using regression-based methods. We combine the LSMC approach with the MLMC method to decrease the computational complexity of the MC method by varying the number of paths on different levels and using the variance reduction method (we refer the reader to [16,17] for a survey), also to obtain more accurate approximation with a lower variance. To illustrate the methodology, we consider a numerical example based on the Black-Scholes model using different values of root mean square error (RMSE). Also, due to large computational times, we deal with a swing option with a small maximum number of exercise rights. is paper is structured as follows: In Section 2, we introduce swing options. We present the LSMC method for the pricing of swing options in Section 3. In Section 4, we give a brief overview of the MLMC method and introduce its algorithm for the pricing of swing options. Section 5 provides a numerical example for pricing swing put options via the MLMC scheme. Also, in this section, we compare our pricing results with those obtained by Carmona and Touzi [8]. Finally, we conclude the paper in Section 6.

Swing Options
Swing options, with multiple exercise rights, play a major role in energy markets with limited and costly storage. ese options offer flexibility regarding when and how much of a commodity is acquired and also provide protection against daily price fluctuations from the beginning until the expiration date of the contract. In other words, a swing option provides the opportunity to repeatedly exercise the right of receiving variable amounts of the commodity at certain prices, regularly. e swing option could be considered as a generalization of the Bermudan option [4].
Swing options allow their holders to exercise one right in each predetermined exercise date from the beginning to the expiration date of the contract. Suppose that the lifetime of a swing contract is the interval [0, T], and the possible exercise , the time between two successive exercise dates), and t 0 � 0. Let (Ω, F, P) be a complete probability space, with F � F t |0 ≤ t ≤ T is the filtration with the σ-algebra F t at time t. Suppose that S t , 0 ≤ t ≤ T is the underlying price process defined on the filtered probability space. Consider Φ(S t ), K t , and ℓ ∈ N denote payoff function, strike price of the option at time t, and maximum number of exercise rights of the option, respectively. e strike price may be constant across all t or it may follow a stochastic process. e value function of a swing option with ℓ exercise rights at time 0 is calculated using the following formula: where and r is the risk-free interest rate.

LSMC Method for Pricing Swing Options
e LSMC algorithm for the pricing of swing options is based on stochastic dynamic programming (SDP) which approximates the conditional expectations within the SDP scheme using least-squares estimates. e idea is based on dynamic programming which begins from the maturity to the start time. At each time t k , two quantities must be compared which are the payoff from immediate exercise and the CV. CV is the conditional expectation of the option payoff with respect to the risk-neutral pricing measure Q: where C t (S t ) and ] t (S t ), respectively, are the CV and option values at time t when the price is S t . In order to approximate the CV in (3), the least-squares regression on a finite set of basis functions is utilized. Basis functions can be monomials, Laguerre polynomials, splines, and radial basis functions. By estimating the CV, one can find the optimal exercise strategy along each path on each exercise date.
Recursive equations (2)-(4) for swing options with ℓ remaining exercise rights are rewritten as follows: 2 Mathematical Problems in Engineering To use the LSMC algorithm for valuing swing options, we have to modify the original algorithm proposed in [13] with an extra dimension in cash flow matrices as described by Dörr in [15]. Below, we briefly present the steps of the LSMC algorithm for valuing swing options: Step 1. A certain number of sample paths of the underlying asset is simulated and stored in a matrix S.
Step 2. Payoffs from the immediate exercise at each time step is stored in a matrix P.
Step 3. Cash flow (CF) matrices of the remaining exercise rights is determined at each time step where the number of determined CF matrices at each time step is equal to ℓ. e aforementioned matrices should be updated at each time step by moving backward in time from T. In CF matrices, the number of rows is equal to the number of sample paths, and the number of columns is equal to the number of exercise times. For obtaining CF matrices, following steps are needed: (i) According to equation (5), the value of CV equals zero at the maturity, and the CF matrix of the κ remaining exercise rights, κ � ℓ, ℓ − 1, . . . , 1, is evaluated by initializing the last κ columns of the CF matrix according to the last κ columns of the matrix P. e other columns of the CF matrix remain undefined.
(ii) CF matrices at time t < T in which T − t < ℓ are evaluated as described for the maturity time evaluation. (iii) For t � t n− 1 , . . . , t 1 in which T − t ≥ ℓ, CF matrices should be updated according to the decision making in (7) i.e., Comparing two quantities: the values of CV and payoff from immediate exercise plus CV with one exercise right less than the original one: For estimating the CV at time t with κ remaining exercise rights, first, we choose in-the-money (ITM) paths (i.e., paths at the t-th column of the matrix S with positive corresponding elements in the matrix P). en, we regress discounted values of the CF matrix with κ exercise rights at time t + 1 onto the price vector at time t and estimate the CV vector. After estimating the CV, we compare quantities in (7) and make decisions. If it is optimal to exercise the right at time t and path u, the (u, t) entry of the CF matrix, i.e., the element of the u-th row and t-th column of CF, will be replaced with the (u, t) entry of the matrix P. e other paths' elements at time t are replaced with 0 and the next columns of CF are replaced with corresponding elements of the CF matrix at time t + 1 with one exercise right less. Otherwise, the CF matrix equals the same one at time t + 1.
Step 4. After determining the CF matrices at time t 1 and discounting them to the start time, t 0 , the option value will be evaluated by taking the average across paths.

Multilevel Monte Carlo Scheme for Pricing Swing Options
MLMC method was first introduced by Heinrich [18] for parametric integration. Giles proposed it for the estimation of an expected value arising from stochastic differential equations (SDEs) in finance [19]. Applying the MLMC method with Euler-Maruyama (EM) discretization reduces the computational complexity of MC simulations for payoff estimation from O(ε − 3 ) to O(ε − 2 (log ε) 2 ) to achieve a root mean square error (RMSE) of ϵ [17]. e MLMC method applies several levels of l, where l � 0, 1, . . . , L. l � 0 and l � L correspond to the coarsest and finest levels, respectively. Each level uses steps of length h l � (T/M l ) for discretization to simulate paths in order to approximate the SDE, where M is the refinement factor and M ≥ 2. e main idea of the MLMC method is to evaluate the expectation in the finest level as a telescopic sum and to apply the standard MC method to estimate each expectation as follows: where P l denotes the approximation of the payoff function on level l. Considering that Y 0 and Y l , respectively, are the estimates for E[P 0 ] and E[P l − P l− 1 ] according to the MC method, we have where the upper index j denotes the jth realization. It is remarkable that the trajectories of simulating paths in the estimation of E[P l − P l− 1 ] are the same. ey are chosen from the finest level, though different step sizes h l and h l− 1 are applied. Let V 0 and V l denote the variance of P 0 and (P l − P l− 1 ) for l > 1, respectively. V l decreases with l, and after decreasing in variance, very few samples will be required on finer levels to estimate the expected values accurately [20].
Applying the MLMC method requires the discretization of SDE. e Euler-Maruyama (EM) scheme is one of the discretization methods. Consider an SDE of the usual form as follows: Mathematical Problems in Engineering 3 where W t is a Wiener process. For SDE (10), the EM approximations of S t (i.e., S t i for i ∈ 0, 1, . . . , n { }) are defined recursively through S t 0 � S 0 and with e MLMC algorithm for pricing swing options is presented in which we choose N l adaptively, which reduces the computational cost to achieve the desired root mean square error (Algorithm 1).

Numerical Result
In this section, we use the combination of the LSMC and MLMC methods to value a swing put option. We consider a swing put contract with maximum number of exercise rights ℓ � 5. e price process S t follows a geometric Brownian motion given as follows: where W t t≥0 is a Wiener process, r and σ are constant riskfree rate of interest and volatility, respectively. e parameter will be set as follows: S 0 � 100, K � 100, σ � 0.3, and r � 0.05. e life of the contract is in the interval [0, 1], and the payoff function is Φ(S t ) � max(K − S(t), 0). Furthermore, we consider n � 5.
To approximate the CV in (6), we choose a set of functions 1, x, x 2 as basis functions for regression. All simulations are carried out on a computer with Intel (R) Core (TM) i5-4310U CPU @ 2.60 GHz processors with 16.0 GB RAM. To approximate SDE in formula (12), we use the EM discretization defined in formula (11) on the uniform time step h l � (T/M l ) on each level l in order to simulate the price paths. Here, we determine the value of M according to the number of exercise rights and exercise times considered in the problem. Due to the number of exercise rights in the problem, we have to choose M ≥ 5 and also start the multilevel estimators in formula (8) from level 1: In the MLMC method, we consider M � 5 and N init l � 100. For more details of implementing the MLMC method, we give a brief description of the MLMC algorithm for two levels: (i) For the first level, l � 1, we start with simulating N init 1 � 100 sample paths via the EM method and store them in a matrix S, after that, we use the leastsquares regression method presented in Section 3 to estimate the values of P 1 (there are N init 1 values of P 1 ).
We evaluate the variance of this estimation as V 1 and utilize this value to determine the required optimal number of sample paths, N opt 1 , using the equation given in (iii) in Algorithm 1. By determining N opt 1 , we simulate (N opt 1 − N init 1 ) new sample paths via the EM algorithm and use the total number of N opt 1 sample paths for estimating P 1 . Lastly, we use the following formula in order to approximate Y 1 : (ii) For the second level, we implement the same method as described for the first level, just the difference is that, here, we need to approximate Y 2 . erefore, we should use N init 2 and N opt 2 number of sample paths for estimating both values of P 1 and P 2 and finally As a benchmark solution, we use the results of the MCbased method proposed by Carmona and Touzi in [8] Table 1.
In Table 2, we compare the results of pricing the swing options using the MLMC method with L � 3 and ε � 1e − 1.7 and the LSMC method using 500,000 sample paths. e number of sample paths in levels 1 to 3 for obtaining results of the MLMC method in Table 2 is 16085592, 5454601, and 3325175, respectively. Comparing the results of the MLMC and LSMC methods in these tables with the benchmark values reveals that the MLMC method is more accurate than the LSMC method for valuing swing options. It is remarkable that we obtain the results of the MLMC method in Table 2 for L � 3 (because of the lack of memory). Due to implementing regression, the run-time is high in both methods. Figure 1 shows the optimal number of paths for each level according to the formula given in (iii) in Algorithm 1, for different values of ε � 1e − 01, 5e − 02, and 1e − 1.4. e optimal number of paths decreases in each level l in order to achieve the desired accuracy. Also, the required number of sample paths increase for obtaining more accurate approximation. In Figure 2, we report the variance reduction of the multilevel correction (P l − P l− 1 ), which we denote by V l . As shown in this figure, an increase in the number of levels leads to a reduction in the values of variance in the MLMC method. We compare this reduction with the variance of P l , which is the approximation of the option price for the LSMC method using the EM discretization with the time step Input: desired accuracy ε, refinement factor M, expiration date T, number of levels L, and number of exercise rights ℓ. Output: estimation of discounted payoff.
(1) For l � 0, 1, . . . , L do: Choose initial number of sample paths, N init l . If l � 0 (i) Generate N init l sample paths according to a discretization scheme such as (11). (ii) Use the least-squares regression method to estimate P l , then determine V l . (iii) Determine the optimal number of sample paths by the formula: , by the formula in (11). (v) Use the least-squares regression method on N opt l paths to estimate P l . (vi) Estimate Y 0 by (9). else (a) Generate N init l sample paths according to the EM scheme in (11), and use these paths for both levels l and l − 1 (coarsen level l to the level l − 1), we refer the interested readers to [21] for more details. (b) Use the least-squares regression method to estimate P l − P l− 1 , and also determine V l using the following formula:    h l � (T/M l ).

Conclusions
In this paper, we have studied the numerical pricing of the swing put option under a geometric Brownian motion model. To obtain a numerical solution to the problem, we have employed the least-squares regression method along with the multilevel Monte Carlo method. We calculated the optimal number of samples for each level adaptively and ran the LSMC algorithm at each level. By increasing the number of levels, we observed a reduction in the optimal number of sample paths and the variance of multilevel corrections, which led to the accurate valuation of swing options with less computation. However, running the algorithm is hard-demanding due to the backward characteristic of the problem. Indeed, decreasing the optimal number of sample paths at each level is the key benefit of using the MLMC method, which leads to lower cost. To show the MLMC method's efficiency, we conducted numerical results with different targeted values of RMSE and compared the obtained results with the previous study and results from the LSMC method. Several further extensions to this work could be (i) using different discretization methods with the higher-order convergence rate for simulation and investigating their effects on the approximated solutions and run-time of the algorithm and (ii) using different regression methods as represented in [16] and proposing the other types of basis functions in regression.

Data Availability
No data were used to support this study.