Valuation of electricity storage contracts using the COS method

Storage of electricity has become increasingly important, due to the gradual replacement of fossil fuels by more variable and uncertain renewable energy sources. In this paper, we provide details on how to mathematically formalize a corresponding electricity storage contract, taking into account the physical limitations of a storage facility and the operational constraints of the electricity grid. We give details of a valuation technique to price these contracts, where the electricity prices follow a structural model based on a stochastic polynomial process. In particular, we show that the Fourier-based COS method can be used to price the contracts accurately and efficiently.


Introduction
One of the main options for the reduction of greenhouse gasses is the use of renewable energy [28], like wind and solar energy. The current, highly variable, output of these energy sources however creates a great challenge in maintaining a balance between demand and supply and ensuring a reliable and stable electricity network [34]. Among all solutions for the highly variable 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 each technology has its own technical characteristics (see [34], [27], [38]). In addition, new techniques and concepts are being developed that can be used for electricity storage (e.g. Car as Power Plant [31], [37]).
The rapid technological improvements of electricity storage are also increasingly interesting from a financial point of view, i.e., storing electricity when there is a lot of supply (and therefore a low price) and selling when the demand is high (and therefore a high price). The businesseconomic consequences, profitability analyses, technological developments and applications of electricity storages have been thoroughly researched, [40], [32], [21], [8], [13], [12], [7]. In this paper, quantitative research is conducted into the valuation of contracts for storing electrical energy by trading on the electricity market. We formalize the contracts mathematically, taking into account the physical limitations of the electricity storage and the operational constraints of the electricity grid, while allowing the contracts to be valued efficiently by a dynamic pricing algorithm.
Prior to the liberalization of the electricity markets, price fluctuations were minimal and often 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 Pool AS, NYMEX). Electricity prices behave differently from other commodities because electricity cannot yet be stored on a large scale (technological progress could change this) causing the prices to respond more directly to the relationship between supply and demand [20]. This feature 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 Schwartz and Smith (2000) [39] and Lucia and Schwartz (2002) [30], where a two-factor model with shortterm mean reverting effects was presented. Further, Barlow (2002) [2] introduced a tractable structural model, based on the demand and supply, where the electricity price is a linear transformation of an underlying factor, modelled as an Ornstein-Uhlenbeck (OU) process. Many other models have been suggested for pricing in electricity markets, on which the aforementioned models often had a strong influence (cf. [3], [14], [43], [5] for reviews on different models).
The method for modelling the electricity prices used in this paper is introduced by Filipovic, Larsson and Ware (2018) [19]. They explained that with the use of a polynomial map such prices can be accurately modeled, where a polynomial process is used for the underlying factors. 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 a well-known OU process. This structural model provides a framework that shows a general and tractable relationship between the underlying factors of the polynomial process and the resulting electricity prices.
A dynamic pricing algorithm is based on a problem formulation with conditional expectations of the discounted value of the contract under the risk-neutral measure, the so-called risk-neutral valuation formula. There are several ways to compute these conditional expectations. In this paper, we focus on a Fourier-based numerical integration technique, known as the COS method, see [15], [16]. The main idea of the COS method is to approximate the conditional probability density function in the risk-neutral valuation formula by its Fourier cosine series expansion, where we make use of the closed-form relation between the coefficients of the Fourier cosine expansion and the characteristic function. The characteristic function of a polynomial process is generally not available, however, we use a variable transformation so that the COS method can still be employed. Furthermore, with the COS method the Greeks can 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 processes. Subsequently, in Section 2.2 we formalize the electricity storage contracts, taking into account physical and operational constraints. Moreover, a dynamic pricing algorithm is given, allowing the contracts to be priced efficiently, and in Section 3 the valuation of financial derivatives on the electricity market is discussed. In Section 4 the COS method for pricing the electricity storage contracts is presented, particularly, Subsection 4.1 details the pricing under the polynomial process for the electricity prices. The numerical results of the electricity storage contracts are discussed in Section 5, where, next to batteries, the concept Car-Park as Power Plant and the cost optimization of charging electric vehicles is analyzed. Finally, in Section 6 we provide a conclusion reflecting on the results and the methods used.

Electricity price dynamics
The electricity price model, based on a so-called polynomial stochastic process, is defined as 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 modelled in the form [19]: where H(X t ) denotes a vector of basis functions for the space of polynomials preserved by the polynomial process (e.g. H(X t ) = (1, X t , X 2 t , ..., X n t ) for a one-dimensional polynomial of degree n ∈ N) and the elements in vector p are the coefficients that characterize the polynomial map Φ. The construction of these increasing polynomial maps is explained in some detail in Appendix A. Moreover, seasonality can be included by making the coefficients in vector p time-dependent.
Stochastic polynomial processes are relevant in many financial applications, including financial market models for interest rates, commodities and electricity. Filipovic and Larsson [17] provide a mathematical basis for polynomial processes and describe the growing range of applications in finance.

Electricity storage contracts
In this section, the electricity storage contracts will be defined. These are contracts where electricity 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 exercise moment. The actions that can be taken are releasing electricity from the storage, storing electricity or doing nothing.
Compared to financial derivative contracts, as Bermudan or swing options, there are extra features 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 the contract does not comply with the contract conditions.
These additional features make the contract valuation a challenging problem. In the following two subsections, an electricity storage contract with such features is defined and the dynamic pricing algorithm for the contract valuation is presented.

Details electricity storage contract
In this section, a storage contract is defined, in a similar way as in [4], where a gas storage contract was discussed.
Consider a contract with M exercise moments, where t 0 is the initial time and {t 1 , ..., t M } the set of exercise moments, where 0 = t 0 < t 1 < ... < t M = T and the time between exercise moments is evenly distributed, ∆t = ∆t m := t m+1 − t m . The holder of the contract has the right to perform an action at any discrete exercise time. However, unlike most financial derivatives, the storage contract includes an extra time point, t M +1 , at which it is not possible to 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 amount of energy is not present in the storage. The penalty function at the settlement date is denoted by q s t M +1 , S t M +1 , e(t M +1 ) , where the notation e(t m ) represents the amount of energy left in the storage at time t m , for all m ∈ {0, ..., M + 1}.
As described, the holder of the contract can take three actions: do nothing, release electricity 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 does not change, ∆e(t m ) = 0, releasing electricity is defined as a negative change, ∆e(t m ) < 0, and storing electricity as a positive change, ∆e(t m ) > 0. By definition, at initial time t 0 the holder can not take any action and the energy change is zero, ∆e(t 0 ) := 0.
The payoff function depends on the action taken by the holder of the contract, and is defined by: wherec(S tm ) andp(S tm ) are respectively the cost of storing and the profit of releasing electricity as a function of the electricity spot price.
Moreover, the contract includes the efficiency of the storage, e.g. an electric vehicle battery is typically between 75% and 93% efficient depending on the battery type [1]. The efficiency of the storage will be denoted by η, which means it converts η · 100% of the purchased amount into electricity which can be sold. Therefore, the following cost and profit functions are used: These functions imply that we have to buy (1/η) > 1 units of electrical energy in order to sell 1 unit of electrical energy. There are physical limitations to the capacity of the storage, for which we introduce a maximum and minimum energy level for the storage, respectively, e min and e max . The energy level of the electricity storage should satisfy, for all t m , m ∈ {0, ..., M + 1}: (2.4) Furthermore, there are operational restrictions on the minimum and maximum energy level changes at an exercise moment, respectively, i min op and i max op . In most markets, there is also a required minimum energy to be released, denoted by i min market . So, the allowed energy level changes at an exercise moment are limited by: In addition, the endurance of the storage facility should be taken into account, e.g. charging/discharging a battery too quickly is known to reduce the battery lifetime [26]. That is why we have set an interval [i min b , i max b ] for energy changes in which the storage can last as long as possible. A penalty function is introduced in case the holder of the contract wants to make a change in the energy level that lies outside this interval. This penalty function for exercise 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 capacity and the operational constraints. We define the set of allowed actions at time t m , for all m ∈ {1, ..., M }, as follows: 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 0 , denoted by v(t 0 , S t0 ), is given by the discounted future payoff and penalties, where the holder of the contract chooses the optimal allowed action at each exercise moment. Thus, the following pricing problem is considered: where Q is the risk-neutral pricing measure, r the risk-neutral interest rate and the optimal actions are given by the set ∆e * = {∆e * (t 1 ), ..., ∆e * (t M )}.
In Table 1, the characteristics of the electricity storage contract are described. ≥ 0 Penalty of charging/discharging too rapidly q b (∆e) Penalty at settlement date q s (e) The efficiency of the storage η First, the total electricity storage capacity is discretized into N e equally distributed energy levels, 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 levels that 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 known beforehand which energy levels will be processed. Therefore, the contract values need to be computed for all possible energy levels e ∈ E. So, the notation of the contract value at each exercise moment needs to be extended by the level of energy in storage at that time. The contract value at time t m , with e(t m ) energy in storage and the electricity price S tm , is denoted by v(t m , S tm , e(t m )).
On the settlement date, t M +1 , it is not possible for the holder to change the energy level in the storage. However, if the agreed amount of energy is not present in the storage, the holder of the contract has to pay a penalty. Therefore, the value of the contract at time t M +1 is equal to the penalty function: At all exercise moments t m , m ∈ {M, ..., 1}, before the settlement date, the holder of the contract can choose an action ∆e(t m ) ∈ A. The holder will choose the action that ultimately gives the highest value. To make this decision, continuation values are needed to determine the expected value of the contract after choosing a specific action. So, the continuation value depends 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 in storage, by choosing certain actions. As an example, the situation with energy level e(t m ) = 2 and action ∆e(t m ) = 1 will result in the same continuation value as the situation with energy level 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 in the set E. The continuation value, for all possible energy levels e ∈ E, can be computed with the risk-neutral valuation formula, i.e., where F t = σ(S s ; s ≤ t), r the constant interest rate and E Q [ · ] the expectation under the risk-neutral measure Q.
Now 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 tm , chooses the action ∆e ∈ A(t m , e) that gives the highest value, taking into account the penalty function. This results in the following contract value at time t m : v(t m , S tm , e) = max ∆e∈A(tm,e) g t m , S tm , ∆e + c t m , S tm , e + ∆e + q b ∆e , ∀e ∈ E. (2.11) The value of the contract at initial time t 0 is equal to the continuation value, because no actions can be taken at this time: (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: (2.13)

Valuation of the electricity storage contract
Option and contract pricing usually boils down to computing conditional expectations of the discounted values of the financial derivative under the risk-neutral measure, the so-called risk-neutral valuation formula. In financial mathematics, it is an important branch of research to 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 the accuracy and robustness appear crucial [35]. A comparison of various numerical approximation techniques was made in [42], where Fourier-based integration techniques performed among the best for many different, but rather basic, options.
In the next section, we discuss a numerical Fourier-based valuation technique, known as the COS method, to price the electricity storage contract, where the electricity prices are defined by the structural model based on the polynomial stochastic processes. In addition to the option value, there is also relevant information in the option Greeks, i.e., the hedge parameters or risk sensitivities.

Risk-neutral measure
In risk management, as well as in portfolio management, the probability space (Ω, F, P) is typically 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 probability measure P is called the real-world measure.
However, when pricing financial derivatives, the risk-neutral measure Q (also called martingale measure) is used. Under the risk-neutral measure Q the discounted prices of assets are martingales, which is not assured under the real-world measure P, i.e. ∀t ≤ 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 be free of arbitrage and complete. The first FTAP states that a market is arbitrage free if and only if there exists a risk-neutral probability measure equivalent 1 to the real-world measure. The second FTAP states that an arbitrage-free market is complete, which means that every risky derivative can be hedged, if and only if there exists a unique risk-neutral measure that is equivalent to the real-world measure.

Risk-neutral measure for the electricity market
As described in the second FTAP, in complete markets the risk-neutral pricing measure is unique, ensuring only one arbitrage-free price of the option/storage contract. Furthermore, in a complete market, risk can be removed to a large extent by the delta hedging trading strategy. In contrast to a complete market, an incomplete market has many different risk-neutral measures and the property to hedge the risk is not existent.
The electricity market is a typical example of an incomplete market, due to its characteristics, and therefore the risk-neutral measure is not unique. In the literature, there are different ways to deal with this incompleteness. It is possible to define a risk-neutral probability measure, Q, by means of the Esscher transform. Another commonly used approach is to assume that the real-world measure is already risk-neutral, and then directly perform the pricing. This latter approach calibrates the model through implied parameters from a liquid market, however the electricity market is illiquid.
In this paper, we adopt another common approach (see e.g. [30], [4], [6]), it is assumed there exists a risk premium to compensate for risk. This risk premium, defined by λσ(t, X t ), is subtracted from the real drift of the electricity price process, where λ is the market price per unit 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 specific risk-neutral measure.

The COS method
The risk-neutral valuation formula can be written in integral form, where v denotes the option value, x the state variable at time t and y at time T : The COS method, [15][16], can be applied to (underlying) processes for which the characteristic function is available. For processes with affine dynamics, the characteristic function can be derived by solving the Ricatti differential equations [11]. Affinity is not invariant under polynomial transformations [18] and therefore the characteristic function of the model based on polynomial processes, S t = Φ(X t ), generally does not exist. However, by using a transformation, the state variables can be conveniently chosen so that the characteristic function of the underlying process, X t , can be employed.
The conditional probability density function in equation (4.1) is approximated by a truncated Fourier cosine expansion, which uses the characteristic function to recover the Fourier coefficients, as follows [15]: where ∆t = T − t, φ(u|∆t, x) is the characteristic function of f (y|T, t, x), a and b the truncated integration interval, Re{·} the real part of the input argument and implies that the first term of the summation is multiplied by 1 2 . By replacing the conditional density function f (y|T, t, x) in equation (4.1) by its Fourier cosine expansion approximation (4.2) and interchanging the integral and the summation, the following formula is obtained, the so-called COS formula: where coefficients V k are defined by: where 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 = Φ −1 (S t ) = X t and y = Φ −1 (S T ) = X T . A function of the asset price that is invertible fits well with the use of the COS method, especially for the calibration of the price process parameters. A closed-form solution of the coefficients V k is available for various 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 large set of processes used in asset pricing models, including process X t with available characteristic function used in this paper, the characteristic function can be written in the following form: where φ(u|∆t) does not depend on x. For processes with independent and stationary increments, e.g. exponential Lévy 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]: where κ n is the n th -order cumulant of the X t process andL depending on the user-defined tolerance level, typicallyL ∈ [6, 12].

The COS method for electricity storage contracts
In this section, the COS method for electricity storage contracts is discussed. The main idea is 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 for approximating the continuation value with e ∈ E energy in storage is defined, for m ∈ {M + 1, ..., 1}, by: where the coefficients V k (t m , e) are defined by: The value of the electricity contract with e(t 0 ) electricity can be approximated if the coefficients V k (t 1 , e(t 0 )) have been recovered, due to the fact that the initial value of the contract is equal to the continuation value at time t 0 .

The coefficients V k (t m , e) for electricity storage contracts
In this section, the coefficients V k (t m , e), ∀m ∈ {M + 1, ..., 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 the contract 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: (4.10) The Fourier coefficients for Bermudan-style options are obtained by the integral of a maximum of two functions [16]. Such an integral is divided into two parts by means of the earlyexercise point. When determining the coefficients for the electricity contract, V k (t m , e) (4.10), a maximum of Dim A(t m , e) functions is taken.
For this, the integration range [a, b] is split into intervals, [a, , 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 ∈ {0, ..., 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 each element 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, ..., 1}, can be written as: where x 0 = a, x n+1 = b, ∆e * i ∈ A(t m , e) the optimal action on interval [x i , x i+1 ] and the Fourier cosine series coefficients of the payoff function, the continuation value and the penalty function are respectively defined by, ∀e ∈ E and ∆e ∈ A(t m , e): xi g t m , y, ∆e cos kπ y − a b − a dy, (4.12) (4.14) Next, it is shown how the coefficients G k (x i , x i+1 , ∆e), C k x i , x i+1 , e 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.
] with optimal action ∆e * i−1 . By choosing a convenient tolerance level, the accuracy will be almost the same, while the computational efficiency will be significantly enhanced.
The coefficients G k (x i , x i+1 , ∆e) are obtained by substitution of the payoff function of the electricity 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 in general, the characteristic function of the underlying process X t is used, which is assumed to be available. Therefore, the state variables are defined as x = Φ −1 (S tm−1 ) and y = Φ −1 (S tm ). These state variables result in the following payoff function: η ∆e, ∆e > 0, 0, ∆e = 0, −Φ(y)∆e, ∆e < 0. By substitution of (4.15) in equation (4.12), the coefficients G k x i , x i+1 , ∆e can be calculated as follows: (4.16) The coefficients G k x i , x i+1 , ∆e can be computed in closed-form for any finite-order polynomial map Φ(·).
Example 4.2. In this example, the coefficients G k x i , x i+1 , ∆e are computed where the electricity price follows a second-order polynomial model, by choosing K = 1, α 1 = 0 and γ 1 = γ in equation A.2, i.e., After substitution of the second-order polynomial model in (4.16), the integral can be computed. This results in the following coefficients G k x i , x i+1 , ∆e : The coefficients Q k (x 1 , x 2 , ∆e) The penalty function q b (∆e) depends only on the action ∆e ∈ A taken and not on the electricity price. 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 cosine series coefficients of the penalty function: (4.20) • For k > 0  The Fourier coefficients of the penalty function at the settlement date (4.9) are computed similarly.
The coefficients C k (x 1 , x 2 , t m , e) To obtain the coefficients C k (x 1 , x 2 , t m , e), the approximation of the continuation value with the COS formula,ĉ(t m , y, e), defined in (4.7), is used. Insertingĉ(t m , y, e) in equation (4.13) and interchanging summation and integration gives the following approximation of coefficients C k (x 1 , x 2 , t m , e), where it is assumed that the characteristic function can be written as in Equation (4.5): with the coefficientsV l (t m ), for m = {1, ..., M }, given by: where, for m = M + 1, the coefficientV l (t M +1 , e) is computed as in (4.9), and the coefficients M k,l (x 1 , x 2 ) are defined as: From basic calculus, the coefficients M k,l (x 1 , x 2 ) can be rewritten into two parts [45]: where it holds that  [16]. The key to this efficient algorithm is that the matrices M s = {M s k,l (x 1 , x 2 )} N −1 k,l=0 and M c = {M c k,l (x 1 , x 2 )} N −1 k,l=0 have respectively a Toeplitz and Hankel structure. The algorithm to compute these coefficients with this method is described in Algorithm 2 in section 2.3 of [16].

The Greeks with the COS method
The option Greeks can be calculated at almost no additional computational costs, by differentiating the COS formula. With v := v(t 0 , S t0 ), S := S t0 and X := X t0 , we have: where the initial contract value of the electricity storage contract, v(t 0 , S t0 , e) = c(t 0 , S t0 , e), is defined by the COS formula (4.7). It follows that the cosine series expansions of the Greeks ∆, Γ and ν are given by: where the characteristic function is defined as in (4.5). This formula can also be used to determine the Greeks at exercise moments in the future, after the optimal action has been taken, making the option value equal to the continuation value, for fixed energy level and electricity price. The Greeks (4.28) of the electricity storage contract, where the price process follows the polynomial model S t = Φ(X t ), can be computed as we assume the inverse of the polynomial map Φ −1 (S t ) exists.

Results and discussion
In this section, several different electricity storage contracts are valued and the Greeks are computed 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 the average energy level in storage for each contract over time is given. The LSMC method can also be applied to value the electricity storage contracts (see Appendix C).
For pricing various electricity storage contracts the electricity price is modeled here according to a second-order polynomial model, as defined in equation (4.17), i.e., where the underlying process X t follows an OU process (see Appendix B for details): To generate an electricity price with more extreme spikes, a jump process can be added to the OU process.
The parameters for this model are chosen, so that the electricity price modeled is representative of a real electricity market with different volatilities: (5.2) and the risk-free interest rate r = 0.01.
The valuation of the electricity storage contract is done for σ = 0.3 (low volatility), σ = 0.6 (mid-low volatility), σ = 0.9 (mid-high volatility) and σ = 1.2 (high volatility). Note that for a stock market these are fairly high volatility parameters, however the electricity market is much more volatile than the stock market. In Figure 1, six simulations of the electricity price process trajectory S t are shown for each volatility parameter.

Mid-low volatility
Mid-high volatility High volatility

The electricity storage contracts
This section is focused on four different electricity storage contracts. There are many different technologies for large scale electricity storage systems and each technology has its own technical characteristics. Each contract is based on some of these technical characteristics. In addition, expected improvements of the characteristics and new concepts to store electricity are 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.

Contract 1: Standard electricity storage
The most widely used electricity storage system is the rechargeable battery. The current rechargeable battery storage facilities often have capacities between 0.25-50 MWh and an output of 0.1-20 MW, depending on the battery [36]. Furthermore, rechargeable batteries can have an efficiency up to ∼ 95% [22]. In addition, charging and discharging a rechargeable battery too rapidly reduces the battery lifetime, so there is a penalty q b (∆e) < 0 for this. Moreover, the settlement penalty of 350 euro is activated when the energy level in the battery is lower than the level when the contract started.
The characteristics for contract 1 are given in Table 3.

Contract 2: Highly efficient electricity storage
In this contract, a highly efficient electricity storage is considered, with an efficiency of 100%. Batteries of ∼ 100% efficiency already exist, but due to the high manufacturing costs, they are not yet used for electricity storage [7]. The contract characteristics for contract 2 are the same as for contract 1 (Table 3), only the efficiency of the electricity storage is increased to η = 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 variable output 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. A single car cannot participate on the electricity market because there are minimal grid injection rules 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 cars is highly predictable [25]. That is why we consider multiple electric vehicles at the same time that 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 capacity and 90% efficiency per electric vehicle is considered. At each exercise moment, the vehicle's battery can change 25% of its energy level without getting a penalty. This results in a total capacity 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, which is why there is a high penalty if there is not enough energy in the car battery.

Contract 4: Cost optimization of charging electric vehicles
Besides viewing the battery of electric vehicles as a storage facility, the contract can also be used 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 vehicles are continuously parked and need to be charged at the end of a period. In this contract a similar car park as in contract 3 is assumed.
At the end of the contract, the car batteries should be fully charged, therefore there is an increasing penalty if the battery has less than the maximum capacity on the settlement date and if the energy is lower than a certain level, defined by e f ix , there is a fixed higher penalty to pay. This penalty function, where we set e f ix = 6, is defined by: e < e f ix .
The characteristics of this contract are defined in Table 5. In the numerical results, we also compare this contract with a contract where the cars are unable to "release" energy, i.e.

The numerical contract values
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% confidence intervals of the contract values are given, derived with the LSMC method. The confidence interval is constructed as follows: where V is the sample mean of ten experiments,σ the standard deviation and z α/2 the critical Z-value (z α/2 = 1.96 for a 95% confidence interval). The resulting confidence intervals are computed by ten runs with the LSMC method of 25 000 trajectories. Furthermore, the average energy levels and the 95% confidence intervals of the energy levels in the storage are shown. The average energy level gives an indication to the holder of a contract which energy levels must be maintained to get maximum profit. Additionally, the maximum and minimum energy levels of all trajectories used are given.
For the COS method, we use the terms N ∈ {100, 150, 200} in the Fourier series expansion andL = 10 to define the integration interval. In addition, the maximum and minimum energy levels at each exercise moment are taken over all trajectories of the ten runs of the LSMC method.
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). Table 6 shows the numerical results of contract 1, with characteristics defined in Table 3. In addition, in Figure 2 the average energy levels in storage are given, together with the 95% confidence intervals.    Table 6 shows that the values obtained with the COS method converge and are always in the 95% LSMC confidence interval. Furthermore, the values of contract 1 are relatively low, especially for less volatile price processes. The reason for this is that the efficiency of the electricity storage is 95%, so the discounted electricity price must increase by (1/0.95−1)·100%, to make a profit. This price change happens more often when the volatility of the price is high, with the low volatility price process (σ = 0.3) this did not occur, as shown in Figure 2.

Contract 1: Standard electricity storage
Although there is little energy change taking place, a zoom of Figure 2 shows that the strategy 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 being fined 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 this strategy appears to yield more profit then.    As shown in Table 7, the values obtained with the COS method converge and are within the 95% LSMC confidence intervals. The values of contract 2 are significantly higher than those of contract 1, because the 100% battery efficiency allows the holder of the contract to make greater use of small fluctuations in the price process.

Contract 2: Highly efficient electricity storage
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 low volatility parameter to loading electricity first and then releasing at the high volatility parameter. Figure 4 shows that a very similar number of actions is taken for the different volatility parameters. As the settlement date approaches, it is ensured that the energy level returns to the starting energy level to avoid the penalty. The maximum and minimum energy levels in Figure 3 show that at mid-high and high volatility electricity price processes all electricity is sometimes sold and the penalty is accepted.

Contract 3: Car-Park as Power Plant (CPPP)
In Table 8 and Figure 5 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 of the electricity price (this can change due to technical improvements). However, an electricity storage 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 volatility price process is used. The strategy is similar to that of contract 1, where the storage has an efficiency of 95%. A difference with the results for contract 1 is that, due to the high penalty at the settlement date, there are no trajectories in which all electricity has been sold in the end. Table 9 and Figures 6 and 7 give the numerical results of contract 4, where the characteristics are described in Table 5 Table 9: The values of contract 4 obtained with the COS method and the LSMC 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, the cheaper the holder of the contract can charge the batteries by using the fluctuations in the price process. If we charge all cars immediately the costs will be 333.33 euros, taking into account the physical and operational constraints, 2-3 euros can be saved by the price-based charging approach if the electricity price follows the dynamics described in Equations (5.1) and (5.2).

Contract 4: Cost optimization of charging electric vehicles
The average strategy of contract 4 is to gradually buy electricity and fully charge the battery before the settlement date so that no penalty has to be paid. Moreover, with a low to mid-low volatility parameter only electricity is bought (not sold), due to the low efficiency. In contrast, with mid-high and high volatility price processes electricity is also sold, but not to a large extent as shown in Figure 7.

The Greeks of the electricity storage contracts
An advantage of the COS method is that the Greeks can be calculated at almost no additional computational costs on top of the computation of the option value, unlike Monte Carlo methods. Tables 10 and 11 show the Greeks of the electricity contracts at initial time t 0 , with respectively an underlying mid-low and high volatile electricity price. Moreover, in Figure 8 the Greeks of contract 2 are shown, as a heat map, at half the contract time, after the actions have been taken, at different fixed electricity prices with the high volatility price parameter and energy levels.    Delta,∆, measures the change in the contract value resulting from a change in the electricity price. With an electricity storage contract the energy can be bought and sold, and it may thus benefit 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 contract benefits 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 the starting 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 electricity price. 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. In addition, the Gamma also becomes slightly higher if the holder may choose actions with high energy level changes. Gamma is highest for contracts with highly efficient electricity storage and 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 from a 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/2 = 49.8 in Figure 8 are because around this electricity price a choice is made to either sell all and accept the penalty or purchase a high enough energy level to avoid the penalty.

Insights and implications of the numerical results
The numerical results of the contract value obtained with the COS method are all within the LSMC 95% confidence intervals. From this it can be implied that the COS method accurately approximates the contract values.
A numerical check of the method's accuracy shows that the values of the considered contracts determined by the COS method converge exponentially, e.g. shown in Figure 9 for contract 2, which confirms the theoretical findings for early-exercise options in [16]. Although the contract value is accurately approximated by means of the COS method, the CPU time is relatively high compared to pricing other financial derivatives, such as European and Bermudan-style options. The reason for this is that when calculating the electricity storage contract, the integral is split into many parts to determine the maximum of all actions and the continuation value must be calculated for each split, which takes most time. However, each integral can be computed in parallel.
The electricity price dynamics (5.1) have an impact on the electricity contract value. The larger the volatility parameter σ, the more the holder of the electricity storage contract can take advantage of large price fluctuations, resulting in a higher contract value. It is also clearly visible in the Greekν that the parameter σ of the underlying price process has a major influence on the contract value, see Tables 10-11 and Figure 8. In addition, the rate of mean reversion parameter κ was examined, which shows that the contract values typically decrease as κ 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 in Tables 2-4. The numerical results clearly show that the efficiency has a significant influence, this is because when the electricity is bought the discounted price must increase by (1/η − 1) · 100% in order not to make a loss. As a result, the price of the storage contract, which only takes into account trading on the electricity market by buying and selling electricity, is lower for less efficient storage. In addition, the magnitude of energy change per exercise moment is important for the value, the greater the permitted energy change, the greater the value.
Besides the dependence on the contract values, the characteristics have a major influence on the strategy of charging/selling electricity at each exercise moment. When the characteristics of the contract change, the strategy can change completely, as can be seen by comparing the strategies of contract 1 and contract 2. The price dynamics also have an effect on this. Especially the volatility parameter, not only does the strategy completely change by taking a different volatility parameter, see Figure 3, the confidence intervals also become smaller with a smaller volatility parameter.

Conclusion
In this paper, we formalized mathematically contracts for storing electricity by trading on the electricity market, allowing the contracts to be valued efficiently. The contract takes into account the physical limitations of an electricity storage, the operational constraints of the electricity grid and the subsequent actions that a holder of the contract can make on the electricity market.
To value the electricity storage contract, we employed an efficient and robust Fourier-based pricing method, where the electricity prices follow a structural model based on a polynomial process. In particular, we focused on the well-known Ornstein-Uhlenbeck (OU) process as the underlying stochastic process. The numerical results show that the COS method performs well and can price the discussed electricity storage contract accurately. Moreover, with the COS method 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 its own 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 numerical results 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 or electric vehicles. With these contracts, it appears that the efficiency of electricity storage and the volatility of the electricity price are of great influence on whether it is useful to apply electricity storage to make a profit by trading on the electricity market.