# A Multilevel Monte Carlo Method for the Valuation of Swing Options

Authors: Hakimeh Ghodssi-Ghassemabadi; Gholam Hossein Yari
Journal: Mathematical Problems in Engineering
 (2021)

---

## Abstract

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. Then, 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.


---

## Body

## 1. 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. Thompson [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 regime-switching 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]. They 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. They 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. They transformed the main problem into a sequence of standard optimal stopping time problems using dynamic programming.The 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. These 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. The 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. They 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.This paper is structured as follows: In Section2, 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.


## 2. Swing Options

Swing options, with multiple exercise rights, play a major role in energy markets with limited and costly storage. These 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. The 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 interval0,T, and the possible exercise times are t1,t2,…,tn=T, s.t. ti>0 for i=1,2,…,n, δi=ti−ti−1 is the refraction time (i.e., the time between two successive exercise dates), and t0=0. Let Ω,F,ℙ be a complete probability space, with F=Ft|0≤t≤T is the filtration with the σ-algebra Ft at time t. Suppose that St,0≤t≤T is the underlying price process defined on the filtered probability space. Consider ΦSt, Kt, and ℓ∈ℕ denote payoff function, strike price of the option at time t, and maximum number of exercise rights of the option, respectively. The strike price may be constant across all t or it may follow a stochastic process. The value function of a swing option with ℓ exercise rights at time 0 is calculated using the following formula:(1)νℓ0,s=supT0,TℓE∑k=1ℓe−r⋅tdkΦStdk|St0=s,where T0,Tℓ=td1,td2,…,tdℓ|0<tdi≤Tandtdi−tdi−1≥δi∀i=1,2,…,ℓ, with d1,d2,…,dℓ⊂1,2,…,n, and r is the risk-free interest rate.


## 3. LSMC Method for Pricing Swing Options

The 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. The idea is based on dynamic programming which begins from the maturity to the start time. At each timetk, 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:(2)CTST=0,(3)CtSt=EQνt+1St+1|St,fort=tn−1,…,t1,(4)νtSt=maxΦSt,CtSt,where CtSt and νtSt, respectively, are the CV and option values at time t when the price is St. 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:(5)CTℓST=0,(6)CtℓSt=Eνt+1ℓSt+1|St,fort=tn−1,…,t1,(7)νtℓSt=maxΦSt+Ctℓ−1St,CtℓSt.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 matrixS.Step 2.
Payoffs from the immediate exercise at each time step is stored in a matrixP.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ℓ. The 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. The other columns of the CF matrix remain undefined.(ii)
CF matrices at timet<T in which T−t<ℓ are evaluated as described for the maturity time evaluation.(iii)
Fort=tn−1,…,t1 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 timet 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). Then, 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. The 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 timet1 and discounting them to the start time, t0, the option value will be evaluated by taking the average across paths.


## 4. 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ε−2logε2 to achieve a root mean square error (RMSE) of ϵ [17].The MLMC method applies several levels ofl, where l=0,1,…,L. l=0 and l=L correspond to the coarsest and finest levels, respectively. Each level uses steps of length hl=T/Ml for discretization to simulate paths in order to approximate the SDE, where M is the refinement factor and M≥2. The 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:(8)EP^L=EP^0+∑l=1LEP^l−P^l−1,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 EP^0 and EP^l−P^l−1 according to the MC method, we have(9)Y^0=1N0∑j=1N0P^0j,Y^l=1Nl∑j=1NlP^lj−P^l−1j,where the upper index j denotes the jth realization.It is remarkable that the trajectories of simulating paths in the estimation ofEP^l−P^l−1 are the same. They are chosen from the finest level, though different step sizes hl and hl−1 are applied. Let V0 and Vl denote the variance of P^0 and P^l−P^l−1 for l>1, respectively. Vl 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. The Euler–Maruyama (EM) scheme is one of the discretization methods. Consider an SDE of the usual form as follows:(10)dSt=μStdt+σStdWt,S0=s0,where Wt is a Wiener process. For SDE (10), the EM approximations of St (i.e., S^ti for i∈0,1,…,n) are defined recursively through S^t0=S0 and(11)S^ti+1=S^ti+μS^tiΔti+1+σS^tiΔWi+1,with Δti+1=ti+1−ti, and ΔWi+1=Wti+1−Wti for i=0,1,…,n−1.The MLMC algorithm for pricing swing options is presented in which we chooseNl adaptively, which reduces the computational cost to achieve the desired root mean square error (Algorithm 1).Algorithm 1: MLMC Algorithm for the pricing of swing options.
Input: desired accuracy ε, refinement factor M, expiration date T, number of levels L, and number of exercise rights ℓ.Output: estimation of discounted payoff.(1)
Forl=0,1,…,Ldo:Choose initial number of sample paths,Nlinit.Ifl=0(i) GenerateNlinit sample paths according to a discretization scheme such as (11).(ii) Use the least-squares regression method to estimateP^l, then determine Vl.(iii) Determine the optimal number of sample paths by the formula:Nlopt=⌈2ε−2Vlhl∑l=0LVl/hl⌉.(iv) Simulate extra samples,Nlopt−Nlinit, by the formula in (11).(v) Use the least-squares regression method onNlopt paths to estimate P^l.(vi) EstimateY^0 by (9).else(a) GenerateNlinit 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 estimateP^l−P^l−1, and also determine Vl using the following formula:Vl=1/Nlinit−1∑i=1NlinitP^l−P^l−12−1/Nlinit∑i=1NlinitP^l−P^l−12.(c) Determine the optimal number of sample paths using the equation in (iii).(d) GenerateNlopt−Nlinit sample paths by the EM scheme in (11) and use these paths for both levels l and l−1.(e) Use the least-squares regression method onNlopt paths to estimate P^l−P^l−1.(f) EstimateY^l by (9).(2)
End For(3)
Approximate the value by (8).


## 5. 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. The price process St follows a geometric Brownian motion given as follows:(12)dSt=rStdt+σStdWt,where Wtt≥0 is a Wiener process, r and σ are constant risk-free rate of interest and volatility, respectively. The parameter will be set as follows: S0=100, K=100, σ=0.3, and r=0.05. The life of the contract is in the interval 0,1, and the payoff function is ΦSt=maxK−St,0. Furthermore, we consider n=5.To approximate the CV in (6), we choose a set of functions 1,x,x2 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 hl=T/Ml 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:(13)EP^L=EP^1+∑l=2LEP^l−P^l−1.In the MLMC method, we considerM=5 and Nlinit=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 N1init=100 sample paths via the EM method and store them in a matrix S, after that, we use the least-squares regression method presented in Section 3 to estimate the values of P^1 (there are N1init values of P^1). We evaluate the variance of this estimation as V1 and utilize this value to determine the required optimal number of sample paths, N1opt, using the equation given in (iii) in Algorithm 1. By determining N1opt, we simulate N1opt−N1init new sample paths via the EM algorithm and use the total number of N1opt sample paths for estimating P^1. Lastly, we use the following formula in order to approximate Y^1:(14)Y^1=1N1opt∑j=1N1optP^1j.(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 approximateY^2. Therefore, we should use N2init and N2opt number of sample paths for estimating both values of P^1 and P^2 and finally(15)Y^2=1N2opt∑j=1N2optP^2j−P^1j.As a benchmark solution, we use the results of the MC-based method proposed by Carmona and Touzi in [8]. These swing option prices are ν01=9.85, ν02=19.26, ν03=28.80, ν04=38.48, and ν05=48.32. Swing option values using the MLMC method considering ε∈1e−01,5e−02,1e−1.4, L=4, and absolute errors are reported in Table 1.Table 1
Values of swing put option withℓ=5, using the MLMC method L=4.
ν0ℓε=1e−01Errorε=5e−02Errorε=1e−1.4Errorν019.78(0.07)9.63(0.22)9.48(0.37)ν0219.48(0.22)19.26(0.00)19.03(0.23)ν0329.31(0.51)28.87(0.07)28.68(0.12)ν0439.28(0.80)38.57(0.09)38.45(0.03)ν0549.29(0.97)48.60(0.28)48.38(0.06)CPU time (s)545.83085618.09103232.4965In Table2, we compare the results of pricing the swing options using the MLMC method with L=3andε=1e−1.7 and the LSMC method using 500,000 sample paths. The 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.Table 2
Value of swing put option withℓ=5, using the MLMC method (L=3, ε=1e−1.7) and LSMC method.
ν0ℓResults of the MLMCErrorResults of the LSMCErrorν019.93(0.08)10.13(0.28)ν0219.68(0.42)19.76(0.50)ν0329.31(0.51)28.87(1.07)ν0438.98(0.50)37.41(0.08)ν0548.92(0.60)45.36(2.96)CPU time (s)1119.09621165.0012Figure1 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. The 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 Vl. 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 hl=T/Ml.Figure 1
Optimal number of sample paths per levell, i.e., Nlopt, in MLMC method for ε∈1e−01,5e−02,1e−1.4 with L=4.Figure 2
Level log variances: variance of multilevel correction,Pl−Pl−1, using the MLMC method for ε∈1e−01,5e−02,1e−1.4, L=4, and variance of Pl implementing the LSMC method, M=5.


## 6. 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.



---

*Source: 8407324-2021-01-19.xml*