Valuation of electricity storage contracts using the COS method
VValuation of electricity storage contracts using the COS method
Boris C. Boonstra a, ∗ , Cornelis W. Oosterlee a,b a CWI - National Research Institute for Mathematics and Computer Science, Amsterdam, the Netherlands b Delft Institute of Applied Mathematics, Delft University of Technology, Delft, the Netherlands
Abstract
Storage of electricity has become increasingly important, due to the gradual replacement offossil fuels by more variable and uncertain renewable energy sources. In this paper, we providedetails on how to mathematically formalize a corresponding electricity storage contract, takinginto account the physical limitations of a storage facility and the operational constraints ofthe electricity grid. We give details of a valuation technique to price these contracts, wherethe electricity prices follow a structural model based on a stochastic polynomial process. Inparticular, we show that the Fourier-based COS method can be used to price the contractsaccurately and efficiently.
Keywords:
Renewable energy, Valuation, Electricity storage, COS method, Fourier-cosineexpansions, Electricity prices, Energy markets, Polynomial processes
1. Introduction
One of the main options for the reduction of greenhouse gasses is the use of renewableenergy [28], like wind and solar energy. The current, highly variable, output of these energysources however creates a great challenge in maintaining a balance between demand and supplyand ensuring a reliable and stable electricity network [34]. Among all solutions for the highlyvariable output [22], electricity storage is acknowledged as a solution with promising potential[7]. There are many different technologies for large-scale electricity storage systems and eachtechnology has its own technical characteristics (see [34], [27], [38]). In addition, new techniquesand concepts are being developed that can be used for electricity storage (e.g. Car as PowerPlant [31], [37]).The rapid technological improvements of electricity storage are also increasingly interestingfrom a financial point of view, i.e., storing electricity when there is a lot of supply (and thereforea low price) and selling when the demand is high (and therefore a high price). The business-economic consequences, profitability analyses, technological developments and applications ofelectricity storages have been thoroughly researched, [40], [32], [21], [8], [13],[12], [7]. In thispaper, quantitative research is conducted into the valuation of contracts for storing electricalenergy by trading on the electricity market. We formalize the contracts mathematically, takinginto account the physical limitations of the electricity storage and the operational constraintsof the electricity grid, while allowing the contracts to be valued efficiently by a dynamic pricingalgorithm.Prior to the liberalization of the electricity markets, price fluctuations were minimal andoften regulated. Nowadays, electricity and electricity derivatives, such as forwards and options,are traded over-the-counter (OTC) or on electricity exchanges (e.g. APX Group, Nord PoolAS, NYMEX). Electricity prices behave differently from other commodities because electricitycannot yet be stored on a large scale (technological progress could change this) causing the ∗ Corresponding author at Delft Institute of Applied Mathematics, TU Delft, Delft, the Netherlands.
Email addresses: [email protected] (Boris C. Boonstra),
[email protected] (CornelisW. Oosterlee)The views expressed in this paper are the personal views of the authors and do not necessarily reflect theviews or policies of their current or past employers.
Preprint submitted to Elsevier January 11, 2021 a r X i v : . [ q -f i n . M F ] J a n rices to respond more directly to the relationship between supply and demand [20]. Thisfeature makes the prices exhibit unique characteristics, such as seasonality, high volatility,mean reversion and price spikes.The early commonly used stochastic models for electricity prices were developed by Schwartzand Smith (2000) [39] and Lucia and Schwartz (2002) [30], where a two-factor model with short-term mean reverting effects was presented. Further, Barlow (2002) [2] introduced a tractablestructural model, based on the demand and supply, where the electricity price is a linear trans-formation of an underlying factor, modelled as an Ornstein-Uhlenbeck (OU) process. Manyother models have been suggested for pricing in electricity markets, on which the aforemen-tioned models often had a strong influence (cf. [3], [14], [43], [5] for reviews on differentmodels).The method for modelling the electricity prices used in this paper is introduced by Fil-ipovic, Larsson and Ware (2018) [19]. They explained that with the use of a polynomial mapsuch prices can be accurately modeled, where a polynomial process is used for the underlyingfactors. By means of an increasing polynomial map, the desired dynamics for the price process,e.g. creating spikes, can be mimicked, while the underlying polynomial process is typically awell-known OU process. This structural model provides a framework that shows a generaland tractable relationship between the underlying factors of the polynomial process and theresulting electricity prices.A dynamic pricing algorithm is based on a problem formulation with conditional expec-tations of the discounted value of the contract under the risk-neutral measure, the so-calledrisk-neutral valuation formula. There are several ways to compute these conditional expec-tations. In this paper, we focus on a Fourier-based numerical integration technique, knownas the COS method, see [15], [16]. The main idea of the COS method is to approximate theconditional probability density function in the risk-neutral valuation formula by its Fouriercosine series expansion, where we make use of the closed-form relation between the coefficientsof the Fourier cosine expansion and the characteristic function. The characteristic function ofa polynomial process is generally not available, however, we use a variable transformation sothat the COS method can still be employed. Furthermore, with the COS method the Greekscan be easily computed for all points in time. Additionally, the Least Squares Monte Carlo(LSMC) method is used here to validate the results obtained with the COS method.In Section 2 we introduce the stochastic price model based on polynomial stochastic pro-cesses. Subsequently, in Section 2.2 we formalize the electricity storage contracts, takinginto account physical and operational constraints. Moreover, a dynamic pricing algorithm isgiven, allowing the contracts to be priced efficiently, and in Section 3 the valuation of financialderivatives on the electricity market is discussed. In Section 4 the COS method for pricingthe electricity storage contracts is presented, particularly, Subsection 4.1 details the pricingunder the polynomial process for the electricity prices. The numerical results of the electricitystorage contracts are discussed in Section 5, where, next to batteries, the concept Car-Parkas Power Plant and the cost optimization of charging electric vehicles is analyzed. Finally, inSection 6 we provide a conclusion reflecting on the results and the methods used.
2. Electricity prices and contracts
The electricity price model, based on a so-called polynomial stochastic process, is definedas follows. If X t follows a stochastic polynomial process, e.g. a Geometric Brownian Motion(GBM) or an Ornstein-Uhlenbeck (OU) process, then the spot price at time t can be modelledin the form [19]: S t = Φ( X t ) = H ( X t ) (cid:62) p , (2.1)where H ( X t ) denotes a vector of basis functions for the space of polynomials preserved bythe polynomial process (e.g. H ( X t ) = (1 , X t , X t , ..., X nt ) (cid:62) for a one-dimensional polynomial ofdegree n ∈ N ) and the elements in vector p are the coefficients that characterize the polynomialmap Φ. The construction of these increasing polynomial maps is explained in some detail inAppendix A. Moreover, seasonality can be included by making the coefficients in vector p time-dependent. 2tochastic polynomial processes are relevant in many financial applications, including fi-nancial market models for interest rates, commodities and electricity. Filipovic and Larsson[17] provide a mathematical basis for polynomial processes and describe the growing range ofapplications in finance. In this section, the electricity storage contracts will be defined. These are contracts whereelectricity can be sold/bought at fixed moments in time by trading on the electricity market.The holder of the electricity storage contract has to decide which action to take at each exercisemoment. The actions that can be taken are releasing electricity from the storage, storingelectricity or doing nothing.Compared to financial derivative contracts, as Bermudan or swing options, there are extrafeatures to be taken into account with an electricity storage contract:1. Physical limitations of the electricity storage, e.g. capacity and endurance.2. Efficiency of the electricity storage.3. Restrictions to the amount of electricity from the grid that can be released or stored.4. The possibility to have negative payoff by storing electricity.5. A storage contract often includes a penalty function that is activated if the holder of thecontract does not comply with the contract conditions.These additional features make the contract valuation a challenging problem. In the followingtwo subsections, an electricity storage contract with such features is defined and the dynamicpricing algorithm for the contract valuation is presented.
In this section, a storage contract is defined, in a similar way as in [4], where a gas storagecontract was discussed.Consider a contract with M exercise moments, where t is the initial time and { t , ..., t M } the set of exercise moments, where 0 = t < t < ... < t M = T and the time between exercisemoments is evenly distributed, ∆ t = ∆ t m := t m +1 − t m . The holder of the contract hasthe right to perform an action at any discrete exercise time . However, unlike most financialderivatives, the storage contract includes an extra time point, t M +1 , at which it is not possibleto exercise, the so-called settlement date of the contract.On the settlement date, the holder of the contract has to pay a penalty if the agreed amountof energy is not present in the storage. The penalty function at the settlement date is denotedby q s (cid:0) t M +1 , S t M +1 , e ( t M +1 ) (cid:1) , where the notation e ( t m ) represents the amount of energy leftin the storage at time t m , for all m ∈ { , ..., M + 1 } .As described, the holder of the contract can take three actions: do nothing, release elec-tricity or store electricity. These actions correspond to a change in energy level in the storage,denoted as ∆ e ( t m ) = e ( t m +1 ) − e ( t m ). If nothing is done at time t m , the energy level doesnot change, ∆ e ( t m ) = 0, releasing electricity is defined as a negative change, ∆ e ( t m ) <
0, andstoring electricity as a positive change, ∆ e ( t m ) >
0. By definition, at initial time t the holdercan not take any action and the energy change is zero, ∆ e ( t ) := 0.The payoff function depends on the action taken by the holder of the contract, and isdefined by: g ( t m , S t m , ∆ e ( t m )) = − ¯ c ( S t m )∆ e ( t m ) , ∆ e ( t m ) > , , ∆ e ( t m ) = 0 , − ¯ p ( S t m )∆ e ( t m ) , ∆ e ( t m ) < , (2.2)where ¯ c ( S t m ) and ¯ p ( S t m ) are respectively the cost of storing and the profit of releasing electricityas a function of the electricity spot price.Moreover, the contract includes the efficiency of the storage, e.g. an electric vehicle batteryis typically between 75% and 93% efficient depending on the battery type [1]. The efficiency3f the storage will be denoted by η , which means it converts η · c ( S t m ) = S t m η , ¯ p ( S t m ) = S t m . (2.3)These functions imply that we have to buy (1 /η ) > e min and e max . The energylevel of the electricity storage should satisfy, for all t m , m ∈ { , ..., M + 1 } : e min ≤ e ( t m ) ≤ e max . (2.4)Furthermore, there are operational restrictions on the minimum and maximum energy levelchanges at an exercise moment, respectively, i minop and i maxop . In most markets, there is alsoa required minimum energy to be released, denoted by i minmarket . So, the allowed energy levelchanges at an exercise moment are limited by:∆ e ( t m ) ∈ (cid:2) i minop , i minmarket (cid:3) ∪ (cid:2) , i maxop (cid:3) , ∀ m ∈ { , ..., M } , (2.5)where i minop ≤ i minmarket ≤ ≤ i maxop .In addition, the endurance of the storage facility should be taken into account, e.g. charg-ing/discharging a battery too quickly is known to reduce the battery lifetime [26]. That iswhy we have set an interval [ i minb , i maxb ] for energy changes in which the storage can last aslong as possible. A penalty function is introduced in case the holder of the contract wants tomake a change in the energy level that lies outside this interval. This penalty function for exer-cise moment t m depends only on the action ∆ e ( t m ) that is taken and is denoted as q b (∆ e ( t m )).Summarizing the storage limitations, the allowed actions, ∆ e ( t m ), are limited by the ca-pacity and the operational constraints. We define the set of allowed actions at time t m , for all m ∈ { , ..., M } , as follows: A ( t m , e ( t m )) = (cid:8) ∆ e (cid:12)(cid:12) e min ≤ e ( t m ) + ∆ e ≤ e max and ∆ e ∈ (cid:2) i minop , i minmarket (cid:3) ∪ (cid:2) , i maxop (cid:3)(cid:9) , (2.6)and the set of allowed actions without getting a penalty q b ( t m , ∆ e ( t m )) is defined by: D ( t m , e ( t m )) = (cid:8) ∆ e (cid:12)(cid:12) e min ≤ e ( t m ) + ∆ e ≤ e max and ∆ e ∈ (cid:2) i minb , i minmarket (cid:3) ∪ (cid:2) , i maxb (cid:3)(cid:9) . (2.7)Note that the set of allowed actions where a penalty is handed out is given by A \ D .The value of the contract at initial time t , denoted by v ( t , S t ), is given by the discountedfuture payoff and penalties, where the holder of the contract chooses the optimal allowed actionat each exercise moment. Thus, the following pricing problem is considered: v ( t , S t ) = max ∆ e ∗ E Q (cid:34) M (cid:88) m =1 e − rt m (cid:18) g (cid:0) t m , S t m , ∆ e ( t m ) (cid:1) + q b (cid:0) ∆ e ( t m ) (cid:1)(cid:19) + e − rt M +1 q s (cid:0) t M +1 , S t M +1 , e ( t M +1 ) (cid:1)(cid:35) , (2.8)where Q is the risk-neutral pricing measure, r the risk-neutral interest rate and the optimalactions are given by the set ∆ e ∗ = { ∆ e ∗ ( t ) , ..., ∆ e ∗ ( t M ) } .In Table 1, the characteristics of the electricity storage contract are described.4tart date t Number of exercise moments M Time to maturity T Time between exercise moments ∆ t = TM Settlement date t M +1 Start energy level e ( t )Min. capacity e min Max. capacity e max Min. energy level change i minop ≤ i maxop ≥ i minmarket Min. energy level change without penalty i minb ≤ i maxb ≥ q b (∆ e )Penalty at settlement date q s ( e )The efficiency of the storage η Table 1: The electricity storage contract characteristics
This section shows how the electricity storage contract, whose details are described inSection 2.2.1, can be priced by a dynamic pricing algorithm.First, the total electricity storage capacity is discretized into N e equally distributed energylevels, with δ := ( e max − e min ) /N e the change of energy between two consecutive energy levels.This results in the set E := { e min , e min + δ, e min + 2 δ, · · · , e max } of all possible energy levelsthat the storage can contain. It is assumed that the action ∆ e ( t m ) at each exercise moment t m is a multiple of δ .Furthermore, the pricing algorithm is computed backward in time and it is not knownbeforehand which energy levels will be processed. Therefore, the contract values need to becomputed for all possible energy levels e ∈ E . So, the notation of the contract value ateach exercise moment needs to be extended by the level of energy in storage at that time. Thecontract value at time t m , with e ( t m ) energy in storage and the electricity price S t m , is denotedby v ( t m , S t m , e ( t m )).On the settlement date, t M +1 , it is not possible for the holder to change the energy level inthe storage. However, if the agreed amount of energy is not present in the storage, the holderof the contract has to pay a penalty. Therefore, the value of the contract at time t M +1 is equalto the penalty function: v ( t M +1 , S t M +1 , e ) = q s ( t M +1 , S t M +1 , e ) , ∀ e ∈ E. (2.9)At all exercise moments t m , m ∈ { M, ..., } , before the settlement date, the holder of thecontract can choose an action ∆ e ( t m ) ∈ A . The holder will choose the action that ultimatelygives the highest value. To make this decision, continuation values are needed to determinethe expected value of the contract after choosing a specific action. So, the continuation valuedepends on the energy level in storage after the action is taken, e ( t m ) + ∆ e ( t m ) = e ( t m +1 ).Note that, at time t m , the continuation value can be the same for different levels of energy instorage, by choosing certain actions. As an example, the situation with energy level e ( t m ) = 2and action ∆ e ( t m ) = 1 will result in the same continuation value as the situation with energylevel e ( t m ) = 4 and action ∆ e ( t m ) = −
1, assuming that ∆ e ( t m ) ∈ A .The continuation value thus does not have to be determined for every energy level e ( t m )and all corresponding allowed actions ∆ e ( t m ) ∈ A , but only for the possible energy levels inthe set E . The continuation value, for all possible energy levels e ∈ E , can be computed withthe risk-neutral valuation formula, i.e., c ( t m , S t m , e ) = e − r ∆ t E Q [ v ( t m +1 , S t m +1 , e ) |F t m ] , ∀ e ∈ E, (2.10)where F t = σ ( S s ; s ≤ t ), r the constant interest rate and E Q [ · ] the expectation under therisk-neutral measure Q . 5ow that the continuation value has been defined, the value of the contract can be given.The holder of the contract at time t m with e ∈ E electricity in storage and electricity price S t m , chooses the action ∆ e ∈ A ( t m , e ) that gives the highest value, taking into account thepenalty function. This results in the following contract value at time t m : v ( t m , S t m , e ) = max ∆ e ∈A ( t m ,e ) (cid:110) g (cid:0) t m , S t m , ∆ e (cid:1) + c (cid:0) t m , S t m , e + ∆ e (cid:1) + q b (cid:0) ∆ e (cid:1)(cid:111) , ∀ e ∈ E. (2.11)The value of the contract at initial time t is equal to the continuation value, because noactions can be taken at this time: v ( t , S t , e ( t )) = c ( t , S t , e ( t )) = e − r ∆ t E Q [ v ( t , S t , e ( t )) |F t ] . (2.12)We now have all the building blocks to determine the contract value, backward in time.This is summarized in the following dynamic pricing algorithm: v ( t M +1 , S t M +1 , e ) = q s ( t M +1 , S t M +1 , e ) , ∀ e ∈ E,c ( t m , S t m , e ) = e − r ∆ t E Q [ v ( t m +1 , S t m +1 , e ) |F t m ] , ∀ e ∈ E and m ∈ { M, ...., } ,v ( t m , S t m , e ) = max ∆ e ∈A (cid:110) g (cid:0) t m , S t m , ∆ e (cid:1) + c (cid:0) t m , S t m , e + ∆ e (cid:1) + q b (cid:0) ∆ e (cid:1)(cid:111) , ∀ e ∈ E and m ∈ { M, ...., } ,v ( t , S t , e ) = c ( t , S t , e ) , ∀ e ∈ E. (2.13)
3. Valuation of the electricity storage contract
Option and contract pricing usually boils down to computing conditional expectations ofthe discounted values of the financial derivative under the risk-neutral measure, the so-calledrisk-neutral valuation formula. In financial mathematics, it is an important branch of researchto price financial derivatives as quickly and accurately as possible. For calibration of a model,the speed of the computation is essential, while for pricing specific derivative contracts theaccuracy and robustness appear crucial [35]. A comparison of various numerical approximationtechniques was made in [42], where Fourier-based integration techniques performed among thebest for many different, but rather basic, options.In the next section, we discuss a numerical Fourier-based valuation technique, known as theCOS method, to price the electricity storage contract, where the electricity prices are definedby the structural model based on the polynomial stochastic processes. In addition to the optionvalue, there is also relevant information in the option Greeks, i.e., the hedge parameters or risksensitivities.
In risk management, as well as in portfolio management, the probability space (Ω , F , P ) istypically considered, where Ω denotes the sample space of all scenarios, F the σ -algebra on Ωand P the probability measure that assigns a probability to every event ω ∈ Ω. The probabilitymeasure P is called the real-world measure.However, when pricing financial derivatives, the risk-neutral measure Q (also called mar-tingale measure) is used. Under the risk-neutral measure Q the discounted prices of assets aremartingales, which is not assured under the real-world measure P , i.e. ∀ t ≤ T : E Q [ e − r ( T − t ) S T |F t ] = S t , where F t ⊂ F is the natural filtration on (Ω , F ) and S t the price of an asset at time t .The fundamental theorems of asset pricing (FTAPs) give conditions for a market to befree of arbitrage and complete. The first FTAP states that a market is arbitrage free if andonly if there exists a risk-neutral probability measure equivalent to the real-world measure.The second FTAP states that an arbitrage-free market is complete, which means that everyrisky derivative can be hedged, if and only if there exists a unique risk-neutral measure that isequivalent to the real-world measure. Given a space (Ω , F ), two measures are equivalent on (Ω , F ) if they agree on which sets in F have probabilityzero. .1.1. Risk-neutral measure for the electricity market As described in the second FTAP, in complete markets the risk-neutral pricing measure isunique, ensuring only one arbitrage-free price of the option/storage contract. Furthermore, in acomplete market, risk can be removed to a large extent by the delta hedging trading strategy. Incontrast to a complete market, an incomplete market has many different risk-neutral measuresand the property to hedge the risk is not existent.The electricity market is a typical example of an incomplete market, due to its characteris-tics, and therefore the risk-neutral measure is not unique. In the literature, there are differentways to deal with this incompleteness. It is possible to define a risk-neutral probability mea-sure, Q , by means of the Esscher transform. Another commonly used approach is to assumethat the real-world measure is already risk-neutral, and then directly perform the pricing. Thislatter approach calibrates the model through implied parameters from a liquid market, howeverthe electricity market is illiquid.In this paper, we adopt another common approach (see e.g. [30], [4], [6]), it is assumedthere exists a risk premium to compensate for risk. This risk premium, defined by λσ ( t, X t ), issubtracted from the real drift of the electricity price process, where λ is the market price perunit risk linked to the state variable X t and σ ( t, X t ) the volatility parameter of the process.This risk premium, calibrated from observed market data, determines the choice of one specificrisk-neutral measure.
4. The COS method
The risk-neutral valuation formula can be written in integral form, where v denotes theoption value, x the state variable at time t and y at time T : v ( t, x ) = e − r ∆ t E Q [ v ( T, y ) |F t ] = e − r ∆ t (cid:90) ∞−∞ v ( T, y ) f ( y | T, t, x ) dy. (4.1)The COS method, [15][16], can be applied to (underlying) processes for which the charac-teristic function is available. For processes with affine dynamics, the characteristic functioncan be derived by solving the Ricatti differential equations [11]. Affinity is not invariant underpolynomial transformations [18] and therefore the characteristic function of the model basedon polynomial processes, S t = Φ( X t ), generally does not exist. However, by using a transfor-mation, the state variables can be conveniently chosen so that the characteristic function ofthe underlying process, X t , can be employed.The conditional probability density function in equation (4.1) is approximated by a trun-cated Fourier cosine expansion, which uses the characteristic function to recover the Fouriercoefficients, as follows [15]: f ( y | T, t, x ) ≈ b − a N − (cid:88) k =0 (cid:48) Re (cid:26) φ (cid:18) kπb − a (cid:12)(cid:12)(cid:12)(cid:12) ∆ t, x (cid:19) · e − ikπ ab − a (cid:27) cos (cid:18) kπ y − ab − a (cid:19) , (4.2)where ∆ t = T − t , φ ( u | ∆ t, x ) is the characteristic function of f ( y | T, t, x ), a and b the truncatedintegration interval, Re {·} the real part of the input argument and (cid:80) (cid:48) implies that the firstterm of the summation is multiplied by .By replacing the conditional density function f ( y | T, t, x ) in equation (4.1) by its Fouriercosine expansion approximation (4.2) and interchanging the integral and the summation, thefollowing formula is obtained, the so-called COS formula: v ( t, x ) ≈ e − r ∆ t N − (cid:88) k =0 (cid:48) Re (cid:26) φ (cid:18) kπb − a (cid:12)(cid:12)(cid:12)(cid:12) ∆ t, x (cid:19) e ikπ − ab − a (cid:27) V k , (4.3)where coefficients V k are defined by: V k = 2 b − a (cid:90) ba v ( T, y ) cos (cid:18) kπ y − ab − a (cid:19) dy, (4.4)7here the state variables x and y can be any function of respectively the asset prices S t and S T , e.g., in our case, x = Φ − ( S t ) = X t and y = Φ − ( S T ) = X T . A function of the assetprice that is invertible fits well with the use of the COS method, especially for the calibrationof the price process parameters. A closed-form solution of the coefficients V k is available forvarious payoff functions and several choices of the state variables x and y .The COS method will form the basis to value the electricity storage contracts. For a largeset of processes used in asset pricing models, including process X t with available characteristicfunction used in this paper, the characteristic function can be written in the following form: φ ( u | ∆ t, x ) = e iuβx φ ( u | ∆ t ) , (4.5)where φ ( u | ∆ t ) does not depend on x . For processes with independent and stationary incre-ments, e.g. exponential L´evy processes, it holds that β = 1. For the OU process it follows that β = e − κ ∆ t .The integration range [ a, b ] is defined as proposed in [35]:[ a, b ] := (cid:20) κ − ¯ L (cid:113) κ + √ κ , κ + ¯ L (cid:113) κ + √ κ (cid:21) , (4.6)where κ n is the n th -order cumulant of the X t process and ¯ L depending on the user-definedtolerance level, typically ¯ L ∈ [6 , In this section, the COS method for electricity storage contracts is discussed. The main ideais to approximate the continuation values, described in the dynamic pricing algorithm (2.13),with the COS formula (4.3) for each allowed energy level e ∈ E . However, the coefficients V k also depend on the energy level e ∈ E in storage. It follows that the COS formula forapproximating the continuation value with e ∈ E energy in storage is defined, for m ∈ { M +1 , ..., } , by: c ( t m − , x, e ) ≈ ˆ c ( t m − , x, e ) := e − r ∆ t N − (cid:88) k =0 (cid:48) Re (cid:26) φ (cid:18) kπb − a (cid:12)(cid:12)(cid:12)(cid:12) ∆ t, x (cid:19) e ikπ − ab − a (cid:27) V k ( t m , e ) , (4.7)where the coefficients V k ( t m , e ) are defined by: V k ( t m , e ) = 2 b − a (cid:90) ba v ( t m , y, e ) cos (cid:18) kπ y − ab − a (cid:19) dy, ∀ e ∈ E. (4.8)The value of the electricity contract with e ( t ) electricity can be approximated if the coef-ficients V k ( t , e ( t )) have been recovered, due to the fact that the initial value of the contractis equal to the continuation value at time t . V k ( t m , e ) for electricity storage contracts In this section, the coefficients V k ( t m , e ), ∀ m ∈ { M + 1 , ..., } and e ∈ E , defined by formula(4.8), are recovered by means of the dynamic pricing algorithm (2.13).First, the coefficients are recovered at the settlement date t M +1 , where the value of thecontract equals the penalty function, v ( t M +1 , y, e ) = q s ( t M +1 , y, e ). Therefore, the coefficients V k ( t M +1 , e ) are obtained for each energy level as follows: V k ( t M +1 , e ) = 2 b − a (cid:90) ba q s ( t M +1 , y, e ) cos (cid:18) kπ y − ab − a (cid:19) dy, ∀ e ∈ E. (4.9)At moments t m , m ∈ { M, ..., } , with energy level e ∈ E , the holder of the contract choosesan allowed action ∆ e ∈ A ( t m , e ) which ensures the maximum ultimate value. This decisionresults in the contract value v ( t m , y, e ) defined in (2.11). By substitution of this value in (4.8),the coefficients at time t m , m ∈ { M, ..., } , with e ∈ E energy in storage, are obtained by: V k ( t m , e ) = 2 b − a (cid:90) ba max ∆ e ∈A ( t m ,e ) (cid:110) g (cid:0) t m , y, ∆ e (cid:1) + c (cid:0) t m , y, e + ∆ e (cid:1) + q b (cid:0) ∆ e (cid:1)(cid:111) cos (cid:18) kπ y − ab − a (cid:19) dy, ∀ e ∈ E. (4.10)8he Fourier coefficients for Bermudan-style options are obtained by the integral of a max-imum of two functions [16]. Such an integral is divided into two parts by means of the early-exercise point. When determining the coefficients for the electricity contract, V k ( t m , e ) (4.10),a maximum of Dim (cid:0) A ( t m , e ) (cid:1) functions is taken.For this, the integration range [ a, b ] is split into intervals, [ a, x ] , [ x , x ] , [ x , x ] , ..., [ x n , b ],so that the maximum is always taken. This split is based on the so-called optimal actions,∆ e ∗ i ∈ A ( t m , e ), for which the contract value v ( t m , y, e ) results in the maximum on interval[ x i , x i +1 ], for i ∈ { , ..., n } .These intervals with corresponding optimal actions can be found by discretizing [ a, b ], i.e.[ a, a + δ, a + 2 δ, ..., b ], and compare the values v ( t m , y, e ) for all actions ∆ e ∈ A ( t m , e ) at eachelement in the discretization and choose the action which gives the maximum value.Once the intervals and the corresponding optimal actions for splitting the integral are found,the coefficients V k ( t m , e ) (4.10) for all e ∈ E at moment t m , m ∈ { M, ..., } , can be written as: V k ( t m , e ) = 2 b − a (cid:20) (cid:90) x a (cid:16) g (cid:0) t m , y, ∆ e ∗ (cid:1) + c (cid:0) t m , y, e + ∆ e ∗ (cid:1) + q b (cid:0) ∆ e ∗ (cid:1)(cid:17) cos (cid:16) kπ y − ab − a (cid:17) dy + (cid:90) x x (cid:16) g (cid:0) t m , y, ∆ e ∗ (cid:1) + c (cid:0) t m , y, e + ∆ e ∗ (cid:1) + q b (cid:0) ∆ e ∗ (cid:1)(cid:17) cos (cid:16) kπ y − ab − a (cid:17) dy + · · · + (cid:90) bx n (cid:16) g (cid:0) t m , y, ∆ e ∗ n (cid:1) + c (cid:0) t m , y, e + ∆ e ∗ n (cid:1) + q b (cid:0) ∆ e ∗ n (cid:1)(cid:17) cos (cid:16) kπ y − ab − a (cid:17) dy (cid:21) = n (cid:88) i =0 (cid:104) G k (cid:0) x i , x i +1 , ∆ e ∗ i (cid:1) + C k (cid:0) x i , x i +1 , e + ∆ e ∗ i (cid:1) + Q k (cid:0) x i , x i +1 , ∆ e ∗ i (cid:1)(cid:105) , (4.11)where x = a , x n +1 = b , ∆ e ∗ i ∈ A ( t m , e ) the optimal action on interval [ x i , x i +1 ] and theFourier cosine series coefficients of the payoff function, the continuation value and the penaltyfunction are respectively defined by, ∀ e ∈ E and ∆ e ∈ A ( t m , e ): G k (cid:0) x i , x i +1 , ∆ e (cid:1) = 2 b − a (cid:90) x i +1 x i g (cid:0) t m , y, ∆ e (cid:1) cos (cid:16) kπ y − ab − a (cid:17) dy, (4.12) C k (cid:0) x i , x i +1 , t m , e (cid:1) = 2 b − a (cid:90) x i +1 x i c (cid:0) t m , y, e (cid:1) cos (cid:16) kπ y − ab − a (cid:17) dy, (4.13) Q k (cid:0) x i , x i +1 , ∆ e (cid:1) = 2 b − a (cid:90) x i +1 x i q b (cid:0) ∆ e (cid:1) cos (cid:16) kπ y − ab − a (cid:17) dy. (4.14)Next, it is shown how the coefficients G k ( x i , x i +1 , ∆ e ), C k (cid:0) x i , x i +1 , e (cid:1) and Q k ( x i , x i +1 , ∆ e )are computed, where the electricity price follows the model based on the polynomial process, S t = Φ( X t ), as described in Section 2. Remark 4.1.
It is possible to improve the computational efficiency by adding a tolerance levelfor a minimum integration interval of [ x i , x i +1 ], ∀ i ∈ { , ..., n } , i.e. if x i +1 − x i < T OL then[ x i − , x i ] = [ x i − , x i +1 ] with optimal action ∆ e ∗ i − . By choosing a convenient tolerance level,the accuracy will be almost the same, while the computational efficiency will be significantlyenhanced. The coefficients G k ( x i , x i +1 , ∆ e )The coefficients G k ( x i , x i +1 , ∆ e ) are obtained by substitution of the payoff function of theelectricity contract (2.2) in equation (4.12), after which the integral can be calculated.Since the characteristic function of the polynomial model S t = Φ( X t ) is not available ingeneral, the characteristic function of the underlying process X t is used, which is assumed tobe available. Therefore, the state variables are defined as x = Φ − ( S t m − ) and y = Φ − ( S t m ).These state variables result in the following payoff function: g ( t m , y, ∆ e ) = − Φ( y ) η ∆ e, ∆ e > , , ∆ e = 0 , − Φ( y )∆ e, ∆ e < . (4.15)9y substitution of (4.15) in equation (4.12), the coefficients G k (cid:0) x i , x i +1 , ∆ e (cid:1) can be calculatedas follows: G k (cid:0) x i , x i +1 , ∆ e (cid:1) = b − a (cid:82) x i +1 x i − Φ( y ) η ∆ e cos (cid:16) kπ y − ab − a (cid:17) dy, ∆ e > , , ∆ e = 0 , b − a (cid:82) x i +1 x i − Φ( y )∆ e cos (cid:16) kπ y − ab − a (cid:17) dy, ∆ e < . (4.16)The coefficients G k (cid:0) x i , x i +1 , ∆ e (cid:1) can be computed in closed-form for any finite-order polyno-mial map Φ( · ). Example 4.2.
In this example, the coefficients G k (cid:0) x i , x i +1 , ∆ e (cid:1) are computed where theelectricity price follows a second-order polynomial model, by choosing K = 1, α = 0 and γ = γ in equation A.2, i.e., S t = Φ( X t ) = 1 − γ X t + γX t . (4.17)After substitution of the second-order polynomial model in (4.16), the integral can be com-puted. This results in the following coefficients G k (cid:0) x i , x i +1 , ∆ e (cid:1) : • For k = 0 G k (cid:0) x i , x i +1 , ∆ e (cid:1) = b − a ∆ eη (cid:104) − γ y − − γ y (cid:105) x i +1 x i , ∆ e > , , ∆ e = 0 , b − a ∆ e (cid:104) − γ y − − γ y (cid:105) x i +1 x i , ∆ e < , (4.18) • For k > G k (cid:0) x i , x i +1 , ∆ e (cid:1) = b − a ∆ eη π k ( a − b ) (cid:18) sin (cid:16) πk ( a − y ) a − b (cid:17)(cid:0) a ( γ − − ab ( γ −
1) + 2 b ( γ −
1) + π k y ( y − γy + 2 γ ) (cid:1) +2 πk ( a − b )( γ ( y − − y ) cos (cid:16) πk ( y − a ) a − b (cid:17)(cid:19) x i +1 x i , ∆ e > , , ∆ e = 0 , b − a ∆ e π k ( a − b ) (cid:18) sin (cid:16) πk ( a − y ) a − b (cid:17)(cid:0) a ( γ − − ab ( γ −
1) + 2 b ( γ −
1) + π k y ( y − γy + 2 γ ) (cid:1) +2 πk ( a − b )( γ ( y − − y ) cos (cid:16) πk ( y − a ) a − b (cid:17)(cid:19) x i +1 x i , ∆ e < . (4.19) The coefficients Q k ( x , x , ∆ e )The penalty function q b (∆ e ) depends only on the action ∆ e ∈ A taken and not on the electricityprice. Furthermore, the penalty function will only be nonzero if an action ∆ e ∈ A \ D is taken.Substitution of the penalty function in equation (4.14) results in the following Fourier cosineseries coefficients of the penalty function: • For k = 0 Q k (cid:0) x i , x i +1 , ∆ e (cid:1) = (cid:40) , ∆ e ∈ D , b − a q b (∆ e )( x i +1 − x i ) , ∆ e ∈ A \ D . (4.20) • For k > Q k (cid:0) x i , x i +1 , ∆ e (cid:1) = (cid:40) , ∆ e ∈ D , πk q b (∆ e ) (cid:16) sin (cid:16) πk ( a − x i +1 ) a − b (cid:17) − sin (cid:16) πk ( a − x i ) a − b (cid:17)(cid:17) , ∆ e ∈ A \ D . (4.21)10he Fourier coefficients of the penalty function at the settlement date (4.9) are computed sim-ilarly. The coefficients C k ( x , x , t m , e )To obtain the coefficients C k ( x , x , t m , e ), the approximation of the continuation value withthe COS formula, ˆ c ( t m , y, e ), defined in (4.7), is used. Inserting ˆ c ( t m , y, e ) in equation (4.13)and interchanging summation and integration gives the following approximation of coefficients C k ( x , x , t m , e ), where it is assumed that the characteristic function can be written as inEquation (4.5):ˆ C k ( x , x , t m , e ) = e − r ∆ t N − (cid:88) l =0 (cid:48) Re (cid:26) φ (cid:18) kπb − a (cid:12)(cid:12)(cid:12)(cid:12) ∆ t (cid:19) ˆ V l ( t m +1 , e ) M k,l ( x , x ) (cid:27) , (4.22)with the coefficients ˆ V l ( t m ), for m = { , ..., M } , given by:ˆ V l ( t m , e ) = n (cid:88) i =0 (cid:104) G l (cid:0) x i , x i +1 , ∆ e ∗ i (cid:1) + ˆ C l (cid:0) x i , x i +1 , e + ∆ e ∗ i (cid:1) + Q l (cid:0) x i , x i +1 , ∆ e ∗ i (cid:1)(cid:105) , (4.23)where, for m = M + 1, the coefficient ˆ V l ( t M +1 , e ) is computed as in (4.9), and the coefficients M k,l ( x , x ) are defined as: M k,l ( x , x ) = 2 b − a (cid:90) x x e ilπ βy − ab − a cos (cid:18) kπ y − ab − a (cid:19) dy. (4.24)From basic calculus, the coefficients M k,l ( x , x ) can be rewritten into two parts [45]: M k,l ( x , x ) = − iπ (cid:0) M sk,l ( x , x ) + M ck,l ( x , x ) (cid:1) , where it holds that M ck,l ( x , x ) = (cid:40) ( x − x ) πib − a , for k = j = 0 , lβ + k (cid:104) e (( jβ + k ) x − ( l + k ) a ) πib − a − e (( lβ + k ) x − ( l + k ) a ) πib − a (cid:105) , otherwise, (4.25) M sk,l ( x , x ) = (cid:40) ( x − x ) πib − a , for k = j = 0 , lβ − k (cid:104) e (( lβ − k ) x − ( l − k ) a ) πib − a − e (( lβ − k ) x − ( l − k ) a ) πib − a (cid:105) , otherwise. (4.26) Remark 4.3.
For characteristic functions where β = 1 in Equation (4.5), an efficient fastFourier transform (FFT) based algorithm for the computation of the Fourier coefficients ofthe continuation value ˆ C k (4.22) can be used [16]. The key to this efficient algorithm is thatthe matrices M s = {M sk,l ( x , x ) } N − k,l =0 and M c = {M ck,l ( x , x ) } N − k,l =0 have respectively aToeplitz and Hankel structure. The algorithm to compute these coefficients with this methodis described in Algorithm 2 in section 2.3 of [16]. The option Greeks can be calculated at almost no additional computational costs, by dif-ferentiating the COS formula. With v := v ( t , S t ), S := S t and X := X t , we have:∆ := ∂v∂S = ∂v∂X ∂X∂S , Γ := ∂ v∂S = ∂ v∂X (cid:18) ∂X∂S (cid:19) + ∂v∂X ∂ X∂S , ν := ∂v∂σ , (4.27)where the initial contract value of the electricity storage contract, v ( t , S t , e ) = c ( t , S t , e ),is defined by the COS formula (4.7). It follows that the cosine series expansions of the Greeks11, Γ and ν are given by:ˆ∆ ≈ ∂X∂S · e − r ∆ t N − (cid:88) k =0 (cid:48) Re (cid:26) φ (cid:18) kπb − a (cid:12)(cid:12)(cid:12)(cid:12) ∆ t, x (cid:19) e ikπ − ab − a ikπb − a β (cid:27) V k , ˆΓ ≈ (cid:18) ∂X∂S (cid:19) · e − r ∆ t N − (cid:88) k =0 (cid:48) Re (cid:40) φ (cid:18) kπb − a (cid:12)(cid:12)(cid:12)(cid:12) ∆ t, x (cid:19) e ikπ − ab − a (cid:18) ikπb − a β (cid:19) (cid:41) V k + ∂ X∂S · e − r ∆ t N − (cid:88) k =0 (cid:48) Re (cid:26) φ (cid:18) kπb − a (cid:12)(cid:12)(cid:12)(cid:12) ∆ t, x (cid:19) e ikπ − ab − a ikπb − a β (cid:27) V k , ˆ ν ≈ e − r ∆ t N − (cid:88) k =0 (cid:48) Re (cid:40) ∂φ (cid:0) kπb − a (cid:12)(cid:12) ∆ t, x (cid:1) ∂σ e ikπ − ab − a (cid:41) V k , (4.28)where the characteristic function is defined as in (4.5). This formula can also be used todetermine the Greeks at exercise moments in the future, after the optimal action has beentaken, making the option value equal to the continuation value, for fixed energy level andelectricity price.The Greeks (4.28) of the electricity storage contract, where the price process follows thepolynomial model S t = Φ( X t ), can be computed as we assume the inverse of the polynomialmap Φ − ( S t ) exists.
5. Results and discussion
In this section, several different electricity storage contracts are valued and the Greeks arecomputed by the COS method. In addition, by using the Least Squares Monte Carlo (LSMC)method [29, 10], the 95% confidence intervals of the contract values are determined and theaverage energy level in storage for each contract over time is given. The LSMC method canalso be applied to value the electricity storage contracts (see Appendix C).For pricing various electricity storage contracts the electricity price is modeled here accord-ing to a second-order polynomial model, as defined in equation (4.17), i.e., S t = Φ( X t ) = 1 − γ X t + γX t , (5.1)where the underlying process X t follows an OU process (see Appendix B for details): dX t = κ ( θ − X t ) dt + σdW t . To generate an electricity price with more extreme spikes, a jump process can be added to theOU process.The parameters for this model are chosen, so that the electricity price modeled is represen-tative of a real electricity market with different volatilities: γ = 0 . , κ = 0 . , θ = 10 . , σ ∈ { . , . , . , . } , X = 10 , S = Φ( X ) , (5.2)and the risk-free interest rate r = 0 . σ = 0 . σ = 0 . σ = 0 . σ = 1 . S t are shown for each volatility parameter.12 ow volatility Mid-low volatilityMid-high volatility High volatility Figure 1: Six simulations of the electricity price process simulated by the second-orderpolynomial model (5.1) with parameters defined in (5.2).
This section is focused on four different electricity storage contracts. There are manydifferent technologies for large scale electricity storage systems and each technology has its owntechnical characteristics. Each contract is based on some of these technical characteristics. Inaddition, expected improvements of the characteristics and new concepts to store electricityare considered, such as higher battery efficiency, the concept of the Car-Park as Power Plant[37] and cost optimization of charging electric vehicles.Some of the electricity storage contract characteristics will be the same for each contract,these are shown in Table 2. The other contract characteristics will differ per contract.Start date t T M t / t M +1 T + ∆ t Required min. release in market i minmarket − . Table 2: The general electricity storage parameters, used for each discussed contract.
Contract 1: Standard electricity storage
The most widely used electricity storage system is the rechargeable battery. The currentrechargeable battery storage facilities often have capacities between 0.25-50 MWh and anoutput of 0.1-20 MW, depending on the battery [36]. Furthermore, rechargeable batteries canhave an efficiency up to ∼
95% [22]. In addition, charging and discharging a rechargeablebattery too rapidly reduces the battery lifetime, so there is a penalty q b (∆ e ) < e ( t ) 7 MWhMin. capacity e min e max
15 MWhMin. energy level change i minop − i maxop i minb − i maxb η
95 %Penalty of charging/discharging too rapidly q b (∆ e ) for ∆ e ∈ A \ D − e Penalty at settlement date q s ( e ) for e < e ( t ) − e Table 3: The electricity storage characteristics for contract 1.
Contract 2: Highly efficient electricity storage
In this contract, a highly efficient electricity storage is considered, with an efficiency of 100%.Batteries of ∼ η = 100%. Contract 3: Car-Park as Power Plant (CPPP)
The Car-Park as Power Plant is a concept that can be seen as a solution for the highly variableoutput of renewable energy sources. Cars are parked for an average of 96% of the time [24],therefore there is potential for using the batteries of electric vehicles for electricity storage. Asingle car cannot participate on the electricity market because there are minimal grid injectionrules of 100 KWh, while a full electric car has an average capacity of 80 KWh. In addition,a single car may not be continuously available, while the presence of a large number of carsis highly predictable [25]. That is why we consider multiple electric vehicles at the same timethat are seen as one storage facility, e.g. a car park can be used [37].For this contract, a car park of 150 electric vehicles and an average of 80 KWh capacityand 90% efficiency per electric vehicle is considered. At each exercise moment, the vehicle’sbattery can change 25% of its energy level without getting a penalty. This results in a totalcapacity of 12 MWh and an energy change without any penalty of 3 MWh.At the end of the contract, in our setting, the owner of the car must be able to drive, whichis why there is a high penalty if there is not enough energy in the car battery.Start energy level e ( t ) 6 MWhMin. capacity e min e max
12 MWhMin. energy level change i minop − i maxop i minb − i maxb η
90 %Penalty of charging/discharging too rapidly q b (∆ e ) for ∆ e ∈ A \ D − e Penalty at settlement date q s ( e ) for e < e ( t ) − e Table 4: The electricity storage characteristics for contract 3.
Contract 4: Cost optimization of charging electric vehicles
Besides viewing the battery of electric vehicles as a storage facility, the contract can also beused to examine how the battery can be charged in a cheap way. In addition to the savings,the cost optimization eases the pressure on the electricity grid during peak power demand,because the electricity price is high then and lower during periods of low demand [33].This price-based charging approach can be applied to places where multiple electric vehiclesare continuously parked and need to be charged at the end of a period. In this contract a similarcar park as in contract 3 is assumed. 14t the end of the contract, the car batteries should be fully charged, therefore there is anincreasing penalty if the battery has less than the maximum capacity on the settlement dateand if the energy is lower than a certain level, defined by e fix , there is a fixed higher penaltyto pay. This penalty function, where we set e fix = 6, is defined by: q s ( e ) = (cid:40) − · e max − ee max − e fix , e ≥ e fix , − , e < e fix . (5.3)The characteristics of this contract are defined in Table 5. In the numerical results, we alsocompare this contract with a contract where the cars are unable to “release” energy, i.e. i minop = i minb = 0, while the other characteristics remain the same.Start energy level e ( t ) 2 MWhMin. capacity e min e max
12 MWhMin. energy level change i minop − i maxop i minb − i maxb η
90 %Penalty of charging/discharging too rapidly q b (∆ e ) for ∆ e ∈ A \ D − e Penalty at settlement date q s ( e ) see (5.3) e Table 5: The electricity storage characteristics for contract 4.
This section shows the numerical results of the valuation of the electricity storage contracts,which are defined in Section 5.1, obtained with the COS method. Moreover, the 95% confidenceintervals of the contract values are given, derived with the LSMC method. The confidenceinterval is constructed as follows:Confidence Interval = (cid:20) V − z α/ (cid:18) ¯ σ √ (cid:19) , V + z α/ (cid:18) ¯ σ √ (cid:19)(cid:21) , (5.4)where V is the sample mean of ten experiments, ¯ σ the standard deviation and z α/ the criticalZ-value ( z α/ = 1 .
96 for a 95% confidence interval). The resulting confidence intervals arecomputed by ten runs with the LSMC method of 25 000 trajectories.Furthermore, the average energy levels and the 95% confidence intervals of the energylevels in the storage are shown. The average energy level gives an indication to the holder ofa contract which energy levels must be maintained to get maximum profit. Additionally, themaximum and minimum energy levels of all trajectories used are given.For the COS method, we use the terms N ∈ { , , } in the Fourier series expansionand ¯ L = 10 to define the integration interval. In addition, the maximum and minimum energylevels at each exercise moment are taken over all trajectories of the ten runs of the LSMCmethod.Python 3.7.1 is used for all the numerical experiments and the CPU is an Intel(R) Core(TM)i7-8750H CPU (2.20GHz, 2208 Mhz, 6 Cores, 12 Logical Processors). Contract 1: Standard electricity storage
Table 6 shows the numerical results of contract 1, with characteristics defined in Table 3. Inaddition, in Figure 2 the average energy levels in storage are given, together with the 95%confidence intervals. 15ricing method LSMC 95% c.i. COS σ N=25000 N=100 N=150 N=2000.3 [0.0000 , 0.0000] 0.0174 -0.0003 0.00000.6 [-0.0005 , 0.0014] -0.0002 0.0000 0.00000.9 [-0.0051 , 0.0222] 0.0087 0.0091 0.00911.2 [0.1399 , 0.1943] 0.1388 0.1434 0.1433
Table 6: The values of contract 1 obtained with the COS method and the LSMC method.Figure 2: The average energy levels in storage for contract 1, respectively from left toright, low volatility, mid-low volatility, mid-high volatility and high volatility obtainedwith the LSMC method.
Table 6 shows that the values obtained with the COS method converge and are alwaysin the 95% LSMC confidence interval. Furthermore, the values of contract 1 are relativelylow, especially for less volatile price processes. The reason for this is that the efficiency of theelectricity storage is 95%, so the discounted electricity price must increase by (1 / . − · σ = 0 .
3) this did not occur, as shown in Figure 2.Although there is little energy change taking place, a zoom of Figure 2 shows that thestrategy is to slowly start storing electricity around t = 0 .
1, if this is expected to yield a profit.In the end, the energy is usually sold until the starting energy level is reached to avoid beingfined 350 euros at the settlement date. In Figure 2, with mid-high and high volatility levels,all energy is sometimes sold from the battery even if a penalty has to be paid, because thisstrategy appears to yield more profit then.
Contract 2: Highly efficient electricity storage
Table 7 and Figures 3 and 4 state the numerical results for contract 2.16ricing method LSMC 95% c.i. COS σ N=25000 N=100 N=150 N=2000.3 [1.8550 , 1.9254] 1.9301 1.8624 1.86300.6 [3.4642 , 3.6050] 3.4770 3.4640 3.46410.9 [5.2075 , 5.4154] 5.2304 5.2291 5.22911.2 [7.1293 , 7.3802] 7.1323 7.1465 7.1464
Table 7: The values of contract 2 obtained with the COS method and LSMC method.Figure 3: The average energy level in storage for contract 2, respectively from left toright, low volatility, mid-low volatility, mid-high volatility and high volatility obtainedwith the LSMC method. igure 4: The average number of uses for the various actions ∆ e ∈ A obtained with theLSMC method. As shown in Table 7, the values obtained with the COS method converge and are withinthe 95% LSMC confidence intervals. The values of contract 2 are significantly higher thanthose of contract 1, because the 100% battery efficiency allows the holder of the contract tomake greater use of small fluctuations in the price process.The strategy is different from contract 1 and depends heavily on the volatility parameter.We see a gradual shift in Figure 3 from first releasing electricity and then loading at the lowvolatility parameter to loading electricity first and then releasing at the high volatility param-eter. Figure 4 shows that a very similar number of actions is taken for the different volatilityparameters. As the settlement date approaches, it is ensured that the energy level returns tothe starting energy level to avoid the penalty. The maximum and minimum energy levels inFigure 3 show that at mid-high and high volatility electricity price processes all electricity issometimes sold and the penalty is accepted.
Contract 3: Car-Park as Power Plant (CPPP)
In Table 8 and Figure 5 the numerical results are shown for contract 3, intended for the conceptCPPP. Pricing method LSMC 95% c.i. COS σ N=25000 N=100 N=150 N=2000.3 [0.0000 , 0.0000] 0.0114 -0.0002 0.00000.6 [-0.0001 , 0.0000] 0.0000 0.0000 0.00000.9 [-0.0008 , 0.0012] -0.0001 0.0000 0.00001.2 [-0.0044 , 0.0020] 0.0000 0.0004 0.0004
Table 8: The values of contract 3 obtained with the COS method and the LSMC method. igure 5: The average energy level in storage for contract 3 for high volatility obtainedwith the LSMC method. Table 8 shows that the CPPP has almost no profit by trading on the electricity market,where the electricity prices are generated by the dynamics described in Equations (5.1)-(5.2).The reason is that the efficiency is not high enough to make profitable use of the fluctuations ofthe electricity price (this can change due to technical improvements). However, an electricitystorage has additional value, e.g. guaranteeing a more stable and reliable energy system.In Figure 5, the average energy level in the storage is shown, where the high volatilityprice process is used. The strategy is similar to that of contract 1, where the storage has anefficiency of 95%. A difference with the results for contract 1 is that, due to the high penaltyat the settlement date, there are no trajectories in which all electricity has been sold in the end.
Contract 4: Cost optimization of charging electric vehicles
Table 9 and Figures 6 and 7 give the numerical results of contract 4, where the characteristicsare described in Table 5.Pricing method LSMC 95% c.i. COS σ N=25000 N=100 N=150 N=2000.3 [-331.3365 , -331.2007] -331.3426 -331.3153 -331.31600.6 [-330.7876 , -330.5472] -330.7729 -330.7741 -330.77420.9 [-330.3961 , -330.0825] -329.3769 -330.3782 -330.37821.2 [-330.1435 , -329.7515] -330.1309 -330.1443 -330.1442
Table 9: The values of contract 4 obtained with the COS method and the LSMC method. igure 6: The average energy level in storage for contract 4, respectively from left toright, low volatility, mid-low volatility, mid-high volatility and high volatility obtainedwith the LSMC method.Figure 7: The average number of uses for the various actions ∆ e ∈ A obtained with theLSMC method. Table 9 shows that the values obtained with the COS method are always within the 95%LSMC confidence intervals and are converging. The more volatile the electricity price, thecheaper the holder of the contract can charge the batteries by using the fluctuations in theprice process. If we charge all cars immediately the costs will be 333 .
33 euros, taking intoaccount the physical and operational constraints, 2-3 euros can be saved by the price-basedcharging approach if the electricity price follows the dynamics described in Equations (5.1)and (5.2). 20he average strategy of contract 4 is to gradually buy electricity and fully charge thebattery before the settlement date so that no penalty has to be paid. Moreover, with a low tomid-low volatility parameter only electricity is bought (not sold), due to the low efficiency. Incontrast, with mid-high and high volatility price processes electricity is also sold, but not to alarge extent as shown in Figure 7.
An advantage of the COS method is that the Greeks can be calculated at almost no addi-tional computational costs on top of the computation of the option value, unlike Monte Carlomethods. Tables 10 and 11 show the Greeks of the electricity contracts at initial time t , withrespectively an underlying mid-low and high volatile electricity price. Moreover, in Figure 8the Greeks of contract 2 are shown, as a heat map, at half the contract time, after the actionshave been taken, at different fixed electricity prices with the high volatility price parameterand energy levels. Contract 1 Contract 2 Contract 3 Contract 4ˆ∆ 0.0000 0.1663 0.0000 -9.1176ˆΓ 0.0001 0.8336 0.0000 0.4957ˆ ν Table 10: The Greeks of the electricity contracts at initial time, where mid-low volatility( σ = 0 .
6) electricity prices are used, obtained with the COS method.
Contract 1 Contract 2 Contract 3 Contract 4ˆ∆ -0.0443 -0.2294 -0.0003 -9.3865ˆΓ 0.0516 0.4055 0.0003 0.3245ˆ ν Table 11: The Greeks of the electricity contracts at initial time, where high volatility( σ = 1 .
2) electricity prices are used, obtained with the COS method. igure 8: The Greeks of contract 2 at half the contract time t M/ , where high volatility( σ = 1 .
2) electricity prices are used, obtained with the COS method.
Delta, ˆ∆, measures the change in the contract value resulting from a change in the electricityprice. With an electricity storage contract the energy can be bought and sold, and it may thusbenefit from increasing and decreasing prices, therefore Delta can be both positive and negative.With a positive Delta, such as with contract 2 with a mid-low volatility parameter, the contractbenefits from a price increase, and vice versa with a negative Delta.If the starting energy level is lower than the level for which the holder receives a penalty,e.g. with contract 4, the Delta drops rapidly and is negative. The opposite is true if thestarting level is higher than the electricity level for which the penalty is obtained.Gamma ˆΓ is a measure of the change of Delta resulting from a change in the electricityprice. For values of Gamma far from zero, the Delta changes drastically when the price changes,making the contract behave differently after a next price change.The contracts where a highly efficient storage is considered have a higher Gamma value. Inaddition, the Gamma also becomes slightly higher if the holder may choose actions with highenergy level changes. Gamma is highest for contracts with highly efficient electricity storageand for highly volatility price processes.Vega ˆ ν measures the impact of a change in the volatility parameter σ on the contract value.The results in Tables 10 and 11 show that the Vega is highest for contracts that benefit froma highly volatile price process, i.e. high efficiency and allowing bigger energy changes.Furthermore, the high values of Vega and Gamma at S t M/ = 49 . .4. Insights and implications of the numerical results The numerical results of the contract value obtained with the COS method are all within theLSMC 95% confidence intervals. From this it can be implied that the COS method accuratelyapproximates the contract values.A numerical check of the method’s accuracy shows that the values of the considered con-tracts determined by the COS method converge exponentially, e.g. shown in Figure 9 forcontract 2, which confirms the theoretical findings for early-exercise options in [16].
Figure 9: The value of contract 2 obtained with the COS method for different numberof terms N in the Fourier series expansion. Although the contract value is accurately approximated by means of the COS method, theCPU time is relatively high compared to pricing other financial derivatives, such as Europeanand Bermudan-style options. The reason for this is that when calculating the electricity storagecontract, the integral is split into many parts to determine the maximum of all actions and thecontinuation value must be calculated for each split, which takes most time. However, eachintegral can be computed in parallel.The electricity price dynamics (5.1) have an impact on the electricity contract value. Thelarger the volatility parameter σ , the more the holder of the electricity storage contract cantake advantage of large price fluctuations, resulting in a higher contract value. It is alsoclearly visible in the Greek ˆ ν that the parameter σ of the underlying price process has a majorinfluence on the contract value, see Tables 10-11 and Figure 8. In addition, the rate of meanreversion parameter κ was examined, which shows that the contract values typically decreaseas κ increases, because for a higher κ the price process converges faster to the mean value.The contract value also depends on the electricity storage contract characteristics, defined inTables 2-4. The numerical results clearly show that the efficiency has a significant influence, thisis because when the electricity is bought the discounted price must increase by (1 /η − ·
6. Conclusion
In this paper, we formalized mathematically contracts for storing electricity by tradingon the electricity market, allowing the contracts to be valued efficiently. The contract takes23nto account the physical limitations of an electricity storage, the operational constraints ofthe electricity grid and the subsequent actions that a holder of the contract can make on theelectricity market.To value the electricity storage contract, we employed an efficient and robust Fourier-basedpricing method, where the electricity prices follow a structural model based on a polynomialprocess. In particular, we focused on the well-known Ornstein-Uhlenbeck (OU) process as theunderlying stochastic process. The numerical results show that the COS method performs welland can price the discussed electricity storage contract accurately. Moreover, with the COSmethod the Greeks can be easily computed, at any time during the contract’s lifetime.Different types of electricity storage contracts are valued, where each contract had itsown characteristics (i.e. storage efficiency, charge/discharge rate, endurance of the storage).Moreover, we looked at the concept Car Park as Power Plant (CPPP), where the numericalresults indicate that the efficiency of electric vehicle batteries strongly impacts the profitability.Furthermore, the contract can be used to reduce the costs of charging electricity storage orelectric vehicles. With these contracts, it appears that the efficiency of electricity storage andthe volatility of the electricity price are of great influence on whether it is useful to applyelectricity storage to make a profit by trading on the electricity market.
References [1] A. Ahmadian, M. Sedghi, A. Elkamel, M. Fowler, and M. A. Golkar. Plug-in electric vehicle batteriesdegradation modeling for smart grid studies: Review, assessment and conceptual framework.
Renewableand Sustainable Energy Reviews , 81:2609–2624, 2018.[2] M. Barlow. A diffusion model for electricity prices.
Mathematical finance , 12(4):287–298, 2002.[3] F. Benth, J. Benth, and S. Koekebakker.
Stochastic modelling of electricity and related markets , volume 11.World Scientific, 2008.[4] A. Boogert and C. De Jong. Gas storage valuation using a monte carlo method.
The journal of derivatives ,15(3):81–98, 2008.[5] R. Carmona and M. Coulon. A survey of commodity markets and structural models for electricity prices.In
Quantitative Energy Finance , pages 41–83. Springer, 2014.[6] A. Cartea and M. Figueroa. Pricing in electricity markets: a mean reverting jump diffusion model withseasonality.
Applied Mathematical Finance , 12(4):313–335, 2005.[7] H. Chen, T. Cong, W. Yang, C. Tan, Y. Li, and Y. Ding. Progress in electrical energy storage system: Acritical review.
Progress in natural science , 19(3):291–312, 2009.[8] D. Connolly, H. Lund, P. Finn, B. Mathiesen, and M. Leahy. Practical operation strategies for pumpedhydroelectric energy storage (phes) utilising electricity price arbitrage.
Energy Policy , 39(7):4189–4196,2011.[9] J. Doob. The brownian movement and stochastic equations.
Annals of Mathematics , pages 351–369, 1942.[10] U. D¨orr. Valuation of swing options and examination of exercise strategies by monte carlo techniques.
Mathematical Finance , 10:27, 2003.[11] D. Duffie, J. Pan, and K. Singleton. Transform analysis and asset pricing for affine jump-diffusions.
Econometrica , 68(6):1343–1376, 2000.[12] B. Dunn, H. Kamath, and J.-M. Tarascon. Electrical energy storage for the grid: a battery of choices.
Science , 334(6058):928–935, 2011.[13] A. Evans, V. Strezov, and T. J. Evans. Assessment of utility energy storage options for increased renewableenergy penetration.
Renewable and Sustainable Energy Reviews , 16(6):4141–4147, 2012.[14] A. Eydeland and K. Wolyniec.
Energy and power risk management: New developments in modeling,pricing, and hedging , volume 206. John Wiley & Sons, 2003.[15] F. Fang and C. Oosterlee. A novel pricing method for european options based on fourier-cosine seriesexpansions.
SIAM Journal on Scientific Computing , 31(2):826–848, 2009.[16] F. Fang and C. Oosterlee. Pricing early-exercise and discrete barrier options by fourier-cosine seriesexpansions.
Numerische Mathematik , 114(1):27, 2009.[17] D. Filipovi´c and M. Larsson. Polynomial diffusions and applications in finance.
Finance and Stochastics ,20(4):931–972, 2016.[18] D. Filipovi´c and M. Larsson. Polynomial jump-diffusion models.
Swiss Finance Institute Research Paper ,(17-60), 2019.[19] D. Filipovi´c, M. Larsson, and T. Ware. Polynomial processes for power prices.
Swiss Finance InstituteResearch Paper , (18-34), 2018.[20] G. Girish and S. Vijayalakshmi. Determinants of electricity price in competitive power market.
Interna-tional Journal of Business and Management , 8(21):70, 2013.[21] F. Graves, T. Jenkin, and D. Murphy. Opportunities for electricity storage in deregulating markets.
TheElectricity Journal , 12(8):46–56, 1999.[22] V. J¨ulch. Comparison of electricity storage options using levelized cost of storage (lcos) method.
AppliedEnergy , 183:1594–1606, 2016.[23] S. Karlin and L. Shapley. Geometry of reduced moment spaces.
Proceedings of the National Academy ofSciences of the United States of America , 35(12):673, 1949.
24] W. Kempton and S. Letendre. Electric vehicles as a new power source for electric utilities.
TransportationResearch Part D: Transport and Environment , 2(3):157–175, 1997.[25] W. Kempton, J. Tomic, S. Letendre, A. Brooks, and T. Lipman. Vehicle-to-grid power: battery, hybrid,and fuel cell vehicles as resources for distributed electric power in california. 2001.[26] M. Koller, T. Borsche, A. Ulbig, and G. Andersson. Defining a degradation cost function for optimalcontrol of a battery energy storage system. In , pages 1–6. IEEE, 2013.[27] G. Kyriakopoulos and G. Arabatzis. Electrical energy storage systems in electricity generation: Energypolicies, innovative technologies, and regulatory regimes.
Renewable and Sustainable Energy Reviews ,56:1044–1067, 2016.[28] S. Lechtenb¨ohmer, L. Nilsson, M. ˚Ahman, and C. Schneider. Decarbonising the energy intensive basicmaterials industry through electrification–implications for future eu electricity demand.
Dubrovnik , 2015.[29] F. Longstaff and E. Schwartz. Valuing american options by simulation: a simple least-squares approach.
The review of financial studies , 14(1):113–147, 2001.[30] J. Lucia and E. Schwartz. Electricity prices and power derivatives: Evidence from the nordic powerexchange.
Review of derivatives research , 5(1):5–50, 2002.[31] H. Lund and W. Kempton. Integration of renewable energy into the transport and electricity sectorsthrough v2g.
Energy policy , 36(9):3578–3587, 2008.[32] H. Lund and G. Salgi. The role of compressed air energy storage (caes) in future sustainable energysystems.
Energy conversion and management , 50(5):1172–1179, 2009.[33] B. Lunz, Z. Yan, J. Gerschler, and D. Sauer. Influence of plug-in hybrid electric vehicle charging strategieson charging and battery degradation costs.
Energy Policy , 46:511–519, 2012.[34] X. Luo, J. Wang, M. Dooner, and J. Clarke. Overview of current development in electrical energy storagetechnologies and the application potential in power system operation.
Applied energy , 137:511–536, 2015.[35] C. Oosterlee and L. Grzelak.
Mathematical Modeling and Computation in Finance . World Scientific, firstedition, November 2019. ISBN 978-1-78634-794-7.[36] O. Palizban and K. Kauhaniemi. Energy storage systems in modern grids—matrix of technologies andapplications.
Journal of Energy Storage , 6:248–259, 2016.[37] E. Park Lee.
A socio-technical exploration of the Car as Power Plant . PhD thesis, Delft University ofTechnology, 2019.[38] A. Poullikkas. A comparative overview of large-scale battery systems for electricity storage.
Renewableand Sustainable Energy Reviews , 27:778–788, 2013.[39] E. Schwartz and J. Smith. Short-term variations and long-term dynamics in commodity prices.
Manage-ment Science , 46(7):893–911, 2000.[40] I. Staffell and M. Rustomji. Maximising the value of electricity storage.
Journal of Energy Storage ,8:212–225, 2016.[41] Z. Sun. Polynomial processes for energy markets. 2019. PowerPoint presentation.[42] L. von Sydow, L. Josef H¨o¨ok, E. Larsson, E. Lindstr¨om, S. Milovanovi´c, J. Persson, V. Shcherbakov, Y. Sh-polyanskiy, S. Sir´en, J. Toivanen, et al. Benchop–the benchmarking project in option pricing.
InternationalJournal of Computer Mathematics , 92(12):2361–2379, 2015.[43] R. Weron. Electricity price forecasting: A review of the state-of-the-art with a look into the future.
International journal of forecasting , 30(4):1030–1081, 2014.[44] B. Zhang, L. Grzelak, and C. Oosterlee. Efficient pricing of commodity options with early-exercise underthe ornstein–uhlenbeck process.
Applied Numerical Mathematics , 62(2):91–111, 2012.[45] B. Zhang and C. Oosterlee. An efficient pricing algorithm for swing options based on fourier cosineexpansions.
Journal of Computational Finance , 16(4):1–32, 2013. ppendix A. Increasing polynomial map To model spot prices with the polynomial model a polynomial map is needed that providesa 1-1 map from a polynomial process X t to the spot price S t . In order to capture the char-acteristics of the electricity commodity market, increasing polynomial maps Φ( · ) from [0 , ∞ )onto [0 , ∞ ) will be used.Such increasing polynomial maps can be constructed from non-negative polynomials on[0 , ∞ ) [23]. To create these increasing polynomial maps the multiplication of positive quadraticswill be used. These quadratics are defined as [41]:˜ q α,γ ( x ) = α x + (1 − α − γ ) x + γ, (A.1)where the quadratic polynomials ˜ q α,γ are normalized such that: (cid:90) ∞ e − x ˜ q α,γ ( x ) dx = 1 . To produce all the possible quadratic polynomials that are on [0 , ∞ ), the parameter set ( α, γ )is chosen as follows: α = ˆ r cos( ξ ) ,γ = ˆ r sin( ξ ) , where ξ ∈ [0 , π/
2] and ˆ r ∈ (cid:104) , cos( ξ ) + sin( ξ ) + (cid:112) sin(2 ξ ) (cid:105) .With the quadratic polynomials as described above, the increasing polynomial map, Φ( · ) :[0 , ∞ ) (cid:55)−→ [0 , ∞ ), can be constructed from pairs of parameters ( α , γ ) , ..., ( α K , γ K ) where K ∈ N + : Φ( x ) = (cid:90) x K (cid:89) k =1 ˜ q α k ,γ k ( u ) du. (A.2)Note that if α k (cid:54) = 0 ∀ k ∈ [1 , K ] the degree of Φ is 2 K + 1. Furthermore, under this generalconstruction of the polynomial map (A.2) a polynomial process with even order can be obtainedby setting α k = 0 for some k ∈ [1 , K ]. Appendix B. Ornstein-Uhlenbeck process
The Ornstein-Uhlenbeck (OU) process has the following SDE: dX t = κ ( θ − X t ) dt + σdW t , (B.1)where t ≥ κ the rate of mean reversion, θ the long-run mean, σ the volatility of the processand W t the Brownian motion under real-world measure P .Integration from t to t on both sides of the SDE (B.1) yields the following solution of theOU process: X t = X t e − κ ( t − t ) + θ (cid:16) − e − κ ( t − t ) (cid:17) + σe − κ ( t − t ) (cid:90) tt e κs dW s . (B.2)By using a scaled time-transformed Brownian motion an analytic solution of the integral inequation (B.2) can be computed, see [9] for further details.The OU process is normally distributed, i.e. X t ∼ N (cid:0) E [ X t ] , V ar [ X t ] (cid:1) , where the expectedvalue and variance are given by: E [ X t |F ] = X t e − κ ( t − t ) + θ (1 − e − κ ( t − t ) ) ,V ar [ X t |F ] = σ κ (1 − e − κ ( t − t ) ) , (B.3)with initial position X t .The well-known characteristic function of the OU process is given by [44]: φ ( u | x, ∆ t ) = e iuxe − κ ∆ t + A ( u, ∆ t ) , (B.4)where A ( u, ∆ t ) = 14 κ (cid:0) e − κ ∆ t − e − κ ∆ t (cid:1) (cid:0) u σ + ue κ ∆ t (cid:0) uσ − iκθ (cid:1)(cid:1) . ppendix C. LSMC algorithm In this appendix, we briefly explain the Least Squares Monte Carlo (LSMC) algorithm toprice the electricity storage contract. This LSMC method presents a simple, intuitive, yetpowerful approximation to value early-exercise options, e.g. Bermudan options, swing optionsand electricity storage contracts. The main insight behind LSMC is that the conditionalexpectation to compute the continuation value can be approximated by the use of least squaresregression as a function of the simulated price process (see [29] and [10]). An advantage ofLSMC is that it is not model dependent, a multi-factor pricing model can be valued just aseasily as a simpler model.The following notation is used for the algorithm: • S im : The asset price of trajectory i at time t m . • CV i,em : The continuation value with e ∈ E energy level in storage after action ∆ e ∈ A istaken at time t m . • CF i,em : The cash flow at time t m with e ∈ E energy level in storage. • ACF i,em : The accumulated value of the cash flows for the optimal actions realised formoments t k , k ∈ { m, ..., M + 1 } , with energy level e ∈ E . • DACF i,em := e − r ∆ t ACF i,em +1 , the discounted accumulated value of cash flows with e ∈ E energy in storage. • Q ∆ em := q b (∆ e ), the penalty function, which is imposed if an action ∆ e ∈ A ( t m , e ) \D ( t m , e ) is taken. • P O i, ∆ em := g ( t m , S im , ∆ e ), the payoff function if action ∆ e ∈ A is taken, as defined informula (2.2).Moreover, two functions of the Python package NumPy are used for the regression, polyfitand polyval. The polyfit function applies a polynomial regression which minimizes the error ina least squares sense. After the regression coefficients are obtained with the polyfit function,we can use polyval to evaluate the polynomial function on other data points to get a prediction. Algorithm 1
The LSMC algorithm for the electricity contract. for e = 0 , ..., N cap do CF i,eM +1 = q s ( t M +1 , S i ( t M +1 ) , e ) , ∀ i = 1 , ...N ACF i,eM +1 = CF i,eM +1 , ∀ i = 1 , ...N end for for m = M, ..., do X i = S i ( t m ) , ∀ i = 1 , ..., N for e = 0 , ..., N cap do DACF i,e = e − r ∆ t ACF i,em +1 , ∀ i = 1 , ..., N p e = polyfit( X, DACF e , CV i,em = polyval( p e , X im ) , ∀ i = 1 , ..., N end for for e = 0 , ..., N cap do for i = 1 , ..., N do ∆ e i, ∗ = arg max ∆ e ∈A ( t m ,e ) (cid:8) P O i, ∆ em + CV i,e +∆ em + Q ∆ e (cid:9) CF i,em = P O i, ∆ e i, ∗ m + Q ∆ e i, ∗ ACF i,em = CF i,em + e − r ∆ t ACF i,e +∆ e i, ∗ m +1 end for end for end for for e = 0 , ..., N cap do v ( t , S ( t ) , e ) = N (cid:80) Ni =1 DACF i,e = N (cid:80) Ni =1 e − r ∆ t ACF i,e end forend for