Capturing the power options smile by an additive two-factor model for overlapping futures prices
aa r X i v : . [ q -f i n . M F ] O c t Capturing the power options smile by an additivetwo-factor model for overlapping futures prices
Marco Piccirilli , Maren Diane Schmeck , and Tiziano Vargiolu Dublin City University, School of Mathematical Sciences, Glasnevin, Dublin9, Ireland Center for Mathematical Economics, University of Bielefeld,Universit¨atsstaße 25, D-33615 Bielefeld, Germany University of Padova, Department of Mathematics, via Trieste 63, TorreArchimede, I-35121 Padova, ItalyOctober 3, 2019
Abstract
In this paper we introduce an additive two-factor model for electricity fu-tures prices based on Normal Inverse Gaussian L´evy processes, that fulfillsa no-overlapping-arbitrage (NOA) condition. We compute European optionprices by Fourier transform methods, introduce a specific calibration procedurethat takes into account no-arbitrage constraints and fit the model to power op-tion settlement prices of the European Energy Exchange (EEX). We show thatour model is able to reproduce the different levels and shapes of the impliedvolatility (IV) profiles displayed by options with a variety of delivery periods.
Keywords:
Volatility Smile, Overlapping Delivery Periods, Arbitrage, AdditiveModels, Power Options, FFT
JEL Classification:
C13, C14, C32, Q40, G13.
One of the big challenges of this age is to develop efficient storage possibilities forelectricity. As the available storage possibilities are very limited in efficiency andcapacity, a major breakthrough of technology is necessary to make the storabilityof electricity comparable to those of commodities as corn or oil. Thus, at this pointof time one can say (simplifying) that electricity is not storable, at least not in thesense of classical commodities. This has important consequences. One of them isthat electricity futures contracts deliver the underlying not at a point of time, butrather over a specified delivery period of different length: typically one can findmonthly, quarterly and yearly delivery periods. Obviously, one quarter consist ofthree months and one year of four quarters, such that there might be futures with1verlapping delivery periods offered in the market. Compared to other markets,this gives rise to an additional dimension of arbitrage opportunities, which appearswhen trading in futures with overlapping delivery periods. Recently, [10] providenew theoretical results concerning the no-arbitrage theory of overlapping additivefutures models. In this framework, the present paper addresses the arbitrage freepricing and calibration of options on overlapping electricity futures.Trading in electricity derivatives is living a significant expansion in Europe. InJune 2018, the European Energy Exchange increased volumes on its power deriva-tives markets by 28% from 181.2 TWh (June 2017) to 231.1 TWh [24]. In particular,the trading volume in power options experienced a boost of 45%. Additionally, thesecontracts are traded over the counter. Thus, it is not a surprise that also the lit-erature considering electricity options as well as options on other commodities isbooming. One stream of literature considers e.g. seasonality, stochastic volatilityand Samuelson-like effects in pricing of power options (e.g. [11, 12, 26, 31, 32, 33, 37])or commodity options (e.g. [3, 4, 34, 35]). However, the issue of overlapping deliveryperiods is not addressed in the above mentioned articles.When building no-arbitrage models with overlapping deliver periods (from nowon called “no-overlapping-arbitrage” or shorter NOA) the focus must be twofold.First, one has to work in a risk neutral setting where traded contracts are mar-tingales. When such a risk neutral measure exists, the first fundamental theoremof asset pricing guarantees that there are no arbitrage possibilities. In fact, whenconsidering overlapping periods, one must additionally take into account possiblearbitrages arising from trading futures with different delivery periods (cf. [30]). Thefirst paper that takes into account this fact is [29], where the authors apply an ap-proach inspired by LIBOR market models [19]. To the best of our knowledge, thisis also the first attempt to fit a consistent option pricing model for power markets.However, the generic LIBOR approach used in [19] allows to model directly onlyshortest delivery, i.e. monthly, contracts, while resorting to approximation to thedistribution of longest delivery futures. The two-factor model of [29] has been gen-eralized by [17, Section 5] to L´evy models, allowing to describe the implied volatility(IV) surface of option prices. Furthermore, a more general approximation procedureis proposed there, but the results, though accurate for single deliveries, are partiallysatisfactory for longest maturity contracts, probably due to the required approxi-mations. Here, we see a strength of our additive approach, making the necessaryapproximations of the geometric approach useless.Among recent contributions in power options modeling, we mention [37], whouses a structural model that explains the formation of IV skews, and [32], who de-velop an extensive sensitivity analysis of IV patterns. However, differently from oursetting, both papers consider options written on electricity spot prices. Anothercharacteristic that we want to reproduce in our model is prices’ seasonality (see also[18]). Here, we address seasonality in the sense of dependence on the delivery of theunderlying futures. This has been addressed by [26], though in a Gaussian setting,2hich does not allow to capture the different IV levels displayed by options withdifferent strikes.In order to compute option prices, we introduce an arbitrage-free two-factoradditive model for futures prices. In contrast to spot-based models, where thefutures prices are computed by taking the discounted expectation of spot pricesunder a pricing measure, we start by specifying the stochastic dynamics of futuresprices under a risk-neutral measure Q . This is known to literature as the Heath-Jarrow-Morton approach applied to commodity markets (see e.g. [25] for a literaturereview in the case of electricity). Since the parameters of our model will be estimatedfrom the option market, this approach has several advantages as already pointedout e.g. by [29]. For instance, we do not need to choose a pricing measure, or,equivalently, estimate the market price of risk that determines how the stochasticmodel changes from the physical measure to the risk-neutral measure under whichthe options are priced in the market. In fact, we can call our risk neutral measure Q , the risk-neutral measure, since this is implicitly and univocally determined by theoption prices observed in the market. We refer to [12] for an empirical discussionon pricing measures for electricity derivatives.Since our futures price model is based on L´evy processes, it allows to capturethe implied volatility profile described by options with different strikes. The dy-namics arises as the natural generalization of its Gaussian counterpart introducedin [30]. It assumes that futures prices are stirred by two stochastic factors builton Normal Inverse Gaussian (NIG) L´evy processes modulated by deterministic co-efficients depending both on time and delivery period. The NIG distribution is aflexible family of distributions that is very popular in financial modeling (see e.g.[6] for general applications of NIG processes to finance, [7] for an electricity spotprice NIG model and [2] for modeling electricity forwards). The first factor has adelivery-averaged exponential behavior, meant to reproduce the Samuelson effect[27]. The second factor is independent of current time, but varies for contracts withdifferent delivery in a seasonal and no-arbitrage way (cf. [30]) and accounts for afiner reproduction of the term structure of futures’ volatilities. We use an additive model, meaning that we do not consider the log-prices, but we instead model di-rectly the prices. This class of models, as opposed to geometric models, has recentlygained an increasing attention in literature due to many modeling advantages (see,for example, [8, 9, 10, 23, 28]). In the context of option pricing, [9] exploit theadditive structure of their spot dynamics for pricing Asian and spread options byfast Fourier techniques. In our case, the additivity property will allow us to fit ourmodel consistently to all the delivery periods traded in the market. Instead, thecalibration results of [17, Section 5], where the author studied a geometric versionof our model (still based on two NIG factors) were only partially satisfactory. Asalready mentioned at the beginning, this is likely due to the necessity to introducean approximation procedure for the distribution of contracts covering periods ofmore than one month. In our case instead, by considering an additive model, we donot need to introduce such approximations, as we have an exact expression for thecontract’s distribution. 3or pricing options, we follow the classical method of [20], which consists, roughlyspeaking, in computing the Fourier transform of (a suitable modification of) the callpayoff as a function of the strike price so to recover the option value by inversetransform. We recall from [20] two different approaches, that we adapt to the caseof additive models. The first approach is based on the use of an exponential dampingfactor in order to make the option value integrable on the whole real line. Instead,the second approach consists in substracting the time value of the option from itspayoff. In this way we derive semi-analytical expressions (analytical up to numericalintegration) for the option prices, that depend on the characteristic function of theunderlying (see also [13]).We discuss a calibration methodology that we will apply in our empirical study.The calibration happens statically, in the sense that we fix a trading date and ob-serve the market option prices for different strikes and deliveries. We then find theparameters that minimize an objective function representing the distance from ob-served prices to model prices. In the same way, one can alternatively use IVs insteadof option prices. The futures model under consideration is defined in such a waythat there is no possibility of arbitrages from trading in overlapping delivery periods.Because no-arbitrage implies certain relations on the coefficients, this translates intoparameter constraints (cf. the calibration procedure presented in [30]).Since there is not enough liquidity in the market in order to extract informationon the IV surface from traded market quotes, we consider the options settlementprices, that are available for a sufficiently large range of strike prices. Though settle-ment prices do not necessarily represent trades that take place in the market, theycontain information on market expectations. We perform the calibration proceduredescribed above, first, for a one-factor model derived by the two-factor one by set-ting one coefficient to 0, and then for the general two-factor model. We comparethe IVs of both models to the empirical IVs and the one (constant across strikes)generated by Black’s model [11, 16]. As a by-product from the estimation of theone-factor model, we derive that, under the risk-neutral measure, futures prices areleptokurtic and have significantly positive skewness. This is reflected also into theshape of empirical IVs, which display a forward skew (i.e. higher IVs for out-of-the-money calls). This can be interpreted as a “risk premium” paid by option buyersfor securing supply (cf. [15, 36]). Finally, we show that the two-factor model isable to reproduce in a satisfactory way the different levels and shapes of the IVprofiles displayed by all the deliveries traded in the market, by outperforming boththe Black and the one-factor model.This paper is organized as follows. In the forthcoming Section 2 we state ourmodeling framework and discuss the no-overlapping-arbitrage conditions in this set-ting. Building on this, we address options on futures as well as its pricing by fastFourier methods in the following section. Then, Section 4 is devoted to the calibra-tion procedure when considering overlapping delivery periods and Section 5 to theempirical study. Finally, Section 6 concludes.4 Additive multifactor models for futures contracts inan overlapping arbitrage free framework
Let F ( t, T , T n ) denote the price at a given day t of a futures contract which deliversa fixed intensity of electricity over the period [ T , T n ]. This period is divided intosubperiods [ T i , T i +1 ], for i = 1 , · · · , n −
1. Thus, the periods are overlapping. Forexample, the delivery period [ T , T n ] can be one quarter long, and is divided into fourmonthly periods, or one year and it is divided into four quarters. Since contractsexpire right before delivery starts, we have that t ≤ T . We will state the multifactormodel in this framework. Nevertheless, when it is sufficient to consider one periodonly, we will use the period [ T , T ] representative for a single arbitrary period.We introduce a general framework for multifactor additive futures prices (as, forexample, in [13]) and, from this, we focus on a two-factor model inspired by [30],that will be of interest for application. More in detail, we introduce a stochasticevolution, parametrized by the delivery period (i.e. depending, in addition to thetrading day t , also on T , T ), driven by independent L´evy factors. We also computethe corresponding characteristic functions, that will constitute the main ingredientsin the computation of option prices (see Section 4.3). We will assume throughoutthis work that the risk-free interest rate is zero (similar arguments apply in the caseof deterministic flat interest rate after discounting, see [30]). In the rest of the paper, we will consider European options written on futurescontracts written on the delivery period [ T , T ] and we will denote by T < T theexercise date of these options. Therefore, for convenience we express the futuresprice at time T . We assume that, for any time t before the exercise of the option,i.e. 0 ≤ t < T , the futures prices F ( T, T , T ) are given by the following stochasticdifferential equation (here given in integral form) F ( T, T , T ) = F ( t, T , T ) + p X k =1 Z Tt Σ k ( u, T , T ) dW k ( u ) + m X j =1 Z Tt Γ j ( u, T , T ) dJ j ( u )(1)where W k are independent Brownian motions for k = 1 , . . . , p and J j ( u ) = R u R R y e N j ( dy, dv )are independent, pure-jump, centered L´evy processes such that R | y | > y ν j ( dy ) < ∞ for each j = 1 , . . . , m . This integrability assumption on the L´evy measure impliesthat J j are square-integrable martingales with zero expectation (for background onL´evy processes see e.g. [22]). We assume that all the stochastic factors are indepen-dent, so that, in particular, the m Poisson random measures are independent of the p Brownian components. The dynamics above are described under a risk-neutralmeasure Q . The absence of the drift follows from the fact that, by no-arbitrage,futures prices must be martingales under Q (see e.g. [14]).It is often the case, in power markets, that futures written on overlapping periodsare traded simultaneously: this originates from the so-called cascade mechanism , by5hich calendar and quarterly futures are split in contracts spanning smaller periods,see e.g. [10] or [30] for a graphical illustration. In particular, it may happen that, fora period [ T , T n ], the futures F ( t, T , T n ) and ( F ( t, T i , T i +1 )) i =1 ,...,n − are all quotedat the same time t < T . In this case, a naturally arising condition on the value V ( t, K ; T , T n ) of a futures contract with price K over the period [ T , T n ] and thevalue V ( t, K ; T i , T i +1 ) with the same price K is V ( t, K ; T , T n ) = n − X i =1 V ( t, K ; T i , T i +1 ) . (2)That is, the value of receiving electricity for the fixed price K over the period[ T , T n ] has to be the same as receiving the electricity for the same price over allpartial periods. Now, the forward price K = F ( t, T , T n ) is defined to be the pricethat makes the value V ( t, K ; T , T n ) of the contract being equal to zero. Thus, from(2) it follows the following no-overlapping-arbitrage (NOA) condition on the forwardprices: F ( t, T , T n ) = 1 T n − T n − X i =1 ( T i +1 − T i ) F ( t, T i , T i +1 ) . (3)For more of this topic, see [30] or [17]. Considering the dynamics of (1) and (3), theNOA condition on the futures price reads in our frameworkΣ k ( t, T , T n ) = 1 T n − T n − X i =1 ( T i +1 − T i )Σ k ( t, T i , T i +1 ) (4)Γ j ( t, T , T n ) = 1 T n − T n − X i =1 ( T i +1 − T i )Γ j ( t, T i , T i +1 ) (5)for all t < T , k = 1 , . . . , p and j = 1 , . . . , m , see e.g. [10].In order to apply the Fourier transform approach to option pricing, we needthe characteristic function of the underlying process. Consider now a generic pe-riod [ T , T ]. By introducing Z ( t, T, T , T ) := F ( T, T , T ) − F ( t, T , T ) and itscharacteristic function (as a function of v ∈ R )Ψ( t, T, T , T , v ) = E h e ivZ ( t,T,T ,T ) |F t i , we have that (see e.g. [13])log Ψ( t, T, T , T , v ) = − v p X k =1 Z Tt Σ k ( u, T , T ) du + n X j =1 ψ j ( t, T ; v Γ j ( · , T , T )) . (6)The function ψ j ( t, T ; θ ( · )) denotes ψ j ( t, T ; θ ( · )) = Z Tt e ψ j ( θ ( u )) du = Z Tt Z R ( e iθ ( u ) z − − iθ ( u ) z ) ν j ( dz ) du, e ψ j ( θ ) is the cumulant of the L´evy process J j computed in θ ∈ R , i.e. e ψ j ( θ ) = log E h e iθJ j (1) i and ν j is the L´evy measure of J j . In this section we consider a two-factor model of the type (1) based on the NormalInverse Gaussian (NIG) distribution. This is motivated by the framework [10] andarises as a natural generalization of [30]. The model assumes that futures prices arestirred by two stochastic factors built on Normal Inverse Gaussian L´evy processesmodulated by deterministic coefficients. The first factor has a delivery-averagedexponential behavior, meant to reproduce the so-called Samuelson effect. This isan observed feature of prices volatilities, common to many commodity markets,consisting of increasing volatility of prices as time approaches maturity [27]. For ananalysis of the impact of the Samuelson effect on option pricing, see [33] and [11].The second factor is independent of time, but varies for contracts with differentdelivery in a seasonal and no-arbitrage way (see [10, 30]) and accounts for a finerreproduction of the term structure of futures volatilities. We remark that, since ourmodel is additive, by “volatility” we mean the parameter (or function of parameters)that determines the variability of prices and not of log-prices as in geometric models.Building upon (1), we assume that the stochastic evolution of a generic futureprice F ( · , T , T ), delivering in the period [ T , T ], from t to T is described by F ( T, T , T ) = F ( t, T , T ) + Z Tt Γ ( u, T , T ) dJ ( u ) + Γ ( T , T )( J ( T ) − J ( t )) , (7)whereΓ ( u, T , T ) := 1 T − T Z T T γ e − µ ( τ − u ) dτ = γ ( e − µ ( T − u ) − e − µ ( T − u ) ) µ ( T − T ) , (8)Γ ( T , T ) := 1 T − T Z T T γ ( τ ) dτ. (9)This model is in the spirit of [10, 30]. Thus, the special form of the coefficientsarises from the implicitly underlying assumption that F can be written as the av-erage over an underlying artificial futures price with instantaneous delivery. Byfocusing our attention to the first component, the parameter γ ∈ R + models the base volatility , that is the volatility of contracts resulting from the first componentwith distant delivery i.e. T − t → ∞ . This coefficient has an exponential rate givenby µ ∈ R + , which determines an increase of volatility as time approaches deliveryi.e. T − t →
0. This effect is averaged over the delivery period for no-arbitrage argu-ments (as explained in [30]) and mimics the Samuelson effect. Regarding the second7omponent, the function γ : [0 , ∞ ) → R + models the general seasonal behaviorof volatility (it takes high values in periods of high volatility and low values forperiods of low volatility). See [18] for a time change model of seasonal volatility inelectricity markets. The seasonal function γ can be specified either in a parametricor nonparametric fashion. For instance, one can have γ ( τ ) := γ + bτ + m X j =1 (cid:0) a j cos( ωjτ ) + a j +1 sin( ωjτ ) (cid:1) , with ω = 2 π/ m ∈ N , b ∈ R (for capturing possible deterministic linear trends),( γ , a , . . . , a m +1 ) ∈ R m +1 .Note also that the driving L´evy processes are the same for all contracts. This isalso a result of the assumption of an underlying instantaneous futures dynamics(see [10], [30]). Note that we have fully specified the form of Γ ( · , T · , T · ), where thecoefficients γ and µ are independent of the delivery period. With this specification,the NOA condition (5) for Γ is naturally fulfilled. On the other hand, if the NOAcondition has to be satisfied for the chosen form of Γ ( u, T i , T i +1 ) as in (8), it isindeed not possible to choose γ and µ different for each delivery period. Though,we leave unspecified the coefficient of the second factor, Γ , so to have more modelingflexibility that is able to account for a finer reproduction of the term structure. Inthe empirical analysis to come, we will use a non-parametric approach and estimateone value for each delivery period. Thus, the following condition has to be fulfilledΓ ( T , T n ) = 1 T n − T n − X i =1 ( T i +1 − T i )Γ ( T i , T i +1 ) . (10)The coefficients in (8) and (9) modulate the variability of the two stochastic pro-cesses J and J , which are both defined as centered versions of NIG L´evy processes.Let us recall that a L´evy process is called Normal Inverse Gaussian with parameters( α, β, δ, µ ) if its characteristic triplet is ( χ, , ν ) with χ = m + 2 αδπ Z sinh( βx ) K ( αx ) dx, (11) ν ( dy ) = αδπ | y | K ( α | y | ) e βy dy, (12)where K is the modified Bessel function of the third kind with index 1 (in the ter-minology of [1, Section 9.6]), 0 ≤ | β | < α and δ > m ∈ R (see [5]). Given a NIGprocess L , the random variable L (1) is NIG distributed with parameters ( α, β, δ, µ ).The NIG distribution is a subclass of a very flexible family, the Generalized Hyper-bolic distributions, and it can accomodate heavy-tails and skewness. The parameter α rules the tail heaviness of the distribution, β determines the skewness, δ is a scaleparameter and µ indicates the location of the distribution.In general, a NIG process is not centered. Therefore, in order to define J and J we subtract from two general NIG processes L , L the corresponding expected value8multiplied by time). For j = 1 ,
2, let L j be a NIG L´evy process with parameters( α j , β j , δ j , m j ) and characteristic triplet ( χ j , , ν j ) and set J j ( t ) = L j ( t ) − E [ L j ( t )] = L j ( t ) − t χ j + Z | y |≥ y ν j ( dy ) ! . Then, J j is a centered NIG process. In particular, it can be easily shown (cf. thecharacteristic function of J j in (32) to see this) that J j is a NIG process withparameters α j , β j , δ j , − δ j β j q α j − β j ! .As already mentioned, in order to compute the option prices, we will need thecharacteristic function of F ( T, T , T ). This in turn reduces to finding the charac-teristic function of Z ( t, T, T , T ) := F ( T, T , T ) − F ( t, T , T ). We compute it inthe appendix in order not to make the presentation unnecessarily heavy. We consider the pricing of European vanilla options written on futures contracts ofthe type introduced in the previous section. We discuss the method for call options,being the case of puts completely analogous. Let C ( t ; T, K, T , T ) be the price attime t (the observation date) of a call option, with strike price K and exercise time T , that is written on a futures contract with delivery period [ T , T ]. By no-arbitrage(see, for instance, [14]), C ( t ; T, K ) = E [( F ( T ) − K ) + |F t ] , (13)where F ( T ) is the price of the underlying futures at the option exercise as in (1)and F t represents the filtration at time t , i.e. the information flow up to time t . In order to ease the notation, sometimes we will not write the dependence onthe delivery period, which does not come into play in this discussion. We recallthat the expectation is taken under the risk-neutral measure Q , even though it isnot explicitly indicated, and throughout this work Q will be the only probabilitymeasure we will deal with. By applying the definition of conditional expectation,we can write C ( t ; T, K ) = Z + ∞ K ( s − K ) q t,T ( s ) ds, (14)where q t,T is the (risk-neutral) density function of F ( T ) conditioned up to time t . This formula yields an expression for the option value, for instance, when thedistribution of the underlying F allows for an explicit formula for the density, that,moreover, can be integrated against the payoff function in a tractable way (as it isthe case in the Black model [16], where the underlying follows a geometric Brownianmotion).We follow the alternative approach of [20], which consists, roughly speaking, incomputing the Fourier transform of (13) as a function of K (after proper manipula-tions) so to recover the option value by inverse transform. This has been studied by9everal authors and it has been shown to be a very convenient way to compute optionprices in the case that the characteristic function of the underlying is known explic-itly, while the density is not. The starting point of the above mentioned approachis the observation that, as K goes to −∞ , C ( t ; T, K ) → ∞ , so that in particular C ( t ; T, K ) is not integrable as a function of K for large negative values. This meansthat the option value C ( t ; T, K ) does not satisfy the assumptions required for com-puting its Fourier transform. In order to overcome this, we follow [20], who suggesttwo approaches that we here recall and apply to the case of additive models.The first approach requires to modify the option value with a damping term e aK , where a > E [ e aZ ( t,T,T ,T ) ] < + ∞ . Whilethis does not pose problems in Gaussian models, where any a > a would deliver the sameresult, it is a well-known result that “extreme” values for a (i.e., too close to 0 or tothe upper bound) give numerical instabilities: thus, one should also optimize withrespect to a in order to have good numerical results.For the arguments above, we instead choose to use the second approach, whichconsists in substracting the time value of the option . Following [20] (see also [21]),define the modified time value of the option (as a function of the strike K ∈ R ) by z MTt,T ( K ) = C ( t ; T, K ) − ( F ( t ) − K ) + , (15)and, if it is square-integrable, compute its Fourier transform ξ MTt,T ( v ) = Z + ∞−∞ e ivK z MTt,T ( K ) dK. (16)The modified option price (15) is then recovered by Fourier inversion after integrat-ing by parts: this is derived in detail in the appendix. As we have derived in the previous section semi-analytical expressions (analyticalup to integration) for the option prices, we now move to discussing a calibrationmethodology that we are going to apply in our empirical study. First, we discretizethe option value that is given in integral form in (30) (see the appendix) in thedomain of integration. Then, we select a finite grid of strike prices, that consistsin practice of the listed options available in the market for a given underlying.This procedure reduces the valuation problem to the computation of a finite sumof vectors, where each component is the option price for a given strike. Then weintroduce a least squares problem designed to find the parameters that minimize anerror function (for a discussion about the choice of the error function, we refer to[17, Section 5.4]). We can compute it for two possible quantities, the model error on For the sake of completeness, we performed the calibration of Section 5 also with the firstapproach, which led to results analogous to those that we present here, but with the additionalcomplications mentioned above.
The quantity that we have to discretize in (30) takes the following form: z t,T ( K ) = 1 π Z + ∞ Re (cid:0) e − iKv ξ t,T ( v ) (cid:1) dv. (17)First, we choose an upper limit A ∈ R + in the above integration (see [20, Section3.1] for a discussion on how to do this optimally). Then, if we apply a simple Eulerrule to the truncated integral, we find an expression of the form z t,T ( K ) ≈ π N X j =1 Re (cid:0) e − iv j K ξ t,T ( v j ) (cid:1) η, (18)where N ∈ N , η := A/N and v j := η ( j −
1) is the integration step. Given M ∈ N different strikes with granularity κ >
0, we compute the function in (18) for thefollowing values of K : K u := K + κ ( u − , for u = 1 , . . . , M, being K the lowest strike price traded. By plugging this in (18), we get for u =1 , . . . , Mz t,T ( K u ) ≈ w t,T ( K u ) := 1 π N X j =1 Re (cid:16) e − iκη ( j − u − e − iKη ( j − ξ t,T ( η ( j − (cid:17) η. (19)As pointed out in [20], this formula is suitable for the application of the fastFourier transform (FFT) algorithm. In order to do this, one must impose that thenumber of strike prices considered is equal to the number of integration nodes, i.e. M = N (which is typically chosen as a power of 2). Moreover, it must hold that κη = 2 πN , which consists of a trade-off between the grid for the integration and the granularityof strike prices. In particular, since in practice the granularity κ is given, thisequality univocally defines the integration grid as a function of N . However, in ourapplication we will not make use of the FFT algorithm and so, in particular, wewill not impose the above restrictions on N, M, κ, η . This is motivated by the factthat we do not have a significant advantage in the computational complexity of theproblem, being the number of strikes consistently lower than N . Furthermore, thefocus of our work is not on the speed-up of the calibration procedure, but ratheron empirical results, so that we do not exclude the possibility to apply the FFTalgorithm, being still possible from a theoretical point of view.11 .2 Parameters estimation Since we have at disposal formulas in discrete form that can be readily implemented,we can now discuss how to fit the two-factor model to market data. The calibrationhappens statically, taking a snapshot of the market, in the sense that we fix a tradingdate and observe the market option prices for that date. We then find the parametersthat fit a certain distance from theoretical to model prices best. As mentioned atthe beginning of this section, one can alternatively use IVs instead of option prices.The futures model under consideration is defined in a such a way that there isno possibility of arbitrages, also in the case of overlapping delivery periods. No-arbitrage implies certain relations on the coefficients, that translate into parameterconstraints. In order to highlight this, we present the objective function and theparameters set for the following three cases: single underlying, many underlyingsbut non-overlapping delivery periods, and the general case of possibly overlappingdelivery periods.
The first and most used in practice model for option prices written on futures isthe Black model [16]. It assumes that the underlying follows a geometric Brownianmotion as in the Black-Scholes formula. Also, the expression for the call price is verysimilar to the Black-Scholes one, with the futures price replacing the stock price,but with different discounting. Since we are assuming that the risk-free rate is zero,the Black-Scholes and the option price given by the Black model are actually thesame in our case. We recall the Black formula here because, in addition to use itas a benchmark in the upcoming empirical application, we use it to compute theimplied volatility of both market and model prices: C BS ( t ; T, K ) = F ( t, T ) N ( d ) − KN ( d ) , (20)where N denotes the cumulative distribution function of a standard Normal randomvariable and d = log F ( t ) K + σ ( T − t ) σ √ T − t , d = d − σ √ T − t. The implied volatility of a given option with price P and strike K is defined as theonly (positive) number σ such that the Black formula for a strike K and volatility σ (all other quantities being equal) gives the price P . The two-factor model aims toreproduce the IV profile of market option prices, that is the plot of σ with respectto K . Black’s formula yields constant implied volatility with respect to K , whilereal option prices usually display smiles or smirks (i.e. the IV is not constant andshows a certain convexity). Assume that we observe at a certain date t < T the prices c ∗ ( t, T, T , T , K u ) of acall option for M different strike prices K u , with exercise at time T written on a12utures with delivery over [ T , T ]. To fit the model to the observed prices we searchfor the parameters that minimize a certain error function. Specifically, we introducethe following least-squares problem: b θ := arg min θ ∈ Θ M X u =1 | c ( t, T, T , T , K u ) − c ∗ ( t, T, T , T , K u ) | , (21)where the set Θ contains all the parameters appearing in the approximation formulaof the model option price C ( t, T, T , T , K u ) ≈ c ( t, T, T , T , K u ) where c ( t, T, T , T , K u ) := w t,T ( K u ) + ( F ( t, T , T ) − K u ) + (22)where w t,T ( K u ) for u = 1 , . . . , M is the discrete function in (19).If we assume that F ( T, T , T ) is given by a two-factor pure-jump model ofthe form (7) where Γ i are defined in (8) and (9) and J i is a centered NIG L´evyprocess with parameters ( α i , β i , { θ = ( α , α , β , β , µ, γ , Γ ( T , T )) ∈ ( R + ) × ( R +0 ) × R + : 0 ≤ | β j | < α j } . We do not need to indicate the parameter m of the original NIG L´evy process L i (see Section 3.2) since it does not appear in thecentered version. Also, the parameter δ is assumed to be 1 without here because ofthe presence of the multiplying factors Γ ( u, T , T ) and Γ ( T , T ), thanks to theproperty that, given a NIG( α, β, δ ) distributed random variable X , for a constant γ > γX is distributed as a NIG( α/γ, β/γ, δγ ). This means that letting δ varyin the parameters set would result in an overparametrization of the minimizationproblem. Let us assume that we are at time t and observe I call options written on futurescontracts with non-overlapping delivery periods [ T i, , T i, ] (for example Jan/YY,Feb/YY, Mar/YY), exercise at T i and M i strike prices K iu := K i + κ i ( u − , for u = 1 , . . . , M i and i = 1 , . . . , I . In the case of more than one contract, weintroduce the following least-squares problem: b θ := arg min θ ∈ Θ I X i =1 M i X u =1 | c ( t, T i , T i, , T i, , K iu ) − c ∗ ( t, T i , T i, , T i, , K iu ) | . (23)where Θ is now the set of parameters { θ = ( α , α , β , β , µ, γ , { Γ i } Ii =1 ) ∈ ( R + ) × ( R +0 ) × ( R + ) I : 0 ≤ | β j | < α j } . For I = 1 we recover exactly the previous case. In general options traded in power markets are written on I futures contracts withdifferent delivery length and possibly overlapping. For example, one can haveApr/YY, May/YY, Jun/YY, Jul/YY, Q2/YY, Q3/YY, Cal-YY simultaneously13raded. As explained in [30], in order to estimate the parameters in a consistent,i.e. arbitrage-free, way, we have to take into account additional constraints on theparameters. From Equation (9), it can be shown that the parameters Γ i satisfy thefollowing constraintsΓ i = Γ ( T i, , T i, ) = n X j =1 T j, − T j, T i, − T i, Γ ( T j, , T j, ) = n X j =1 T j, − T j, T i, − T i, Γ j (24)whenever [ T i, , T i, ] is the union of disjoint intervals [ T j, , T j, ] for j = 1 , . . . , n , i.e.for all the forwards with overlapping delivery, see Section 2 for details. Let us call atomic , the contracts whose delivery period can not be partitioned by the deliveryperiods of other futures. In other words, we suppose that m forwards F , . . . , F m have non-overlapping delivery periods [ T , , T , ] , . . . , [ T m, , T m, ] and such that thedelivery periods of the other contracts traded in the market can be expressed as unionof the former. For example, assume that we observe in the chosen calibration windowthe option prices for futures contracts delivering over Apr/YY, May/YY, Jun/YY,Jul/YY, Q2/YY, Q3/YY, Cal-YY. On one hand, Q2/YY is not atomic, since itcan be “splitted” into Apr/YY, May/YY and Jun/YY. On the other hand, Q3/YYturns out to be atomic, even if Jul/YY is already traded, as Aug/YY and Sep/YYare not observed. For the same reason, Cal-YY is considered atomic as well. Then,in this example, if Γ , . . . , Γ denote the corresponding parameters, we have thatΓ , Γ , Γ , Γ , Γ , Γ are free parameters, as they refer to atomic contracts, whereasto determine Γ we use Equation (24). Consequently, we define the same statisticsas in (23) but where now the vector of parameters is subject to the additionalconstraints given by Equation (24). For example, with the convention that Γ Q2 / YY denotes the parameter Γ i corresponding to the contract Q2/YY, thenΓ Q2 / YY = u Apr / YY Γ Apr / YY + u May / YY Γ May / YY + u Jun / YY Γ Jun / YY , (25)where the weights u i are defined according to the number of days in the month/quarter(e.g. for Apr/YY we have u Apr / YY = 30 / b θ := arg min θ ∈ Θ I X i =1 M i X u =1 | c ( t, T i , T i, , T i, , K iu ) − c ∗ ( t, T i , T i, , T i, , K iu ) | . (26)where Θ is now the set of parameters θ = ( α , α , β , β , µ, γ , { Γ i } Ii =1 ) ∈ ( R + ) × ( R +0 ) × ( R + ) I such that 0 ≤ | β j | < α j , Γ i = n X j =1 T j, − T j, T i, − T i, Γ j . Empirical study
In this section we describe our dataset, compute and discuss the empirically observedvolatility implied by settlement prices. We compare the performance of three models:a purely Gaussian model, a one-factor model (being a special case of the two-factorone) and the general two-factor model.
The contracts that we consider in our application are European-styled call optionstraded at the EEX Power Derivatives market. The underlying assets are futures con-tracts that prescribe the delivery of 1 MWh for each hour of each day of a month,a quarter or a year. More specifically, we will consider call options written on thePhelix Base index of the German/Austrian area. These options are called PhelixBase Month/Quarter/Year Options and the EEX official product codes are O1BM,O1BQ, O1BY. The term
Base refers to
Base Load , because the delivery of electricitytakes place for each hour of the day, in contrast to
Peak Load contracts, that insteadprescribe the delivery only for the hours from 8 to 20. Usually, the exercise of theoptions under consideration is few trading days before the start of delivery. Sincerecently, yearly options are available for four different exercise dates. We consider inour dataset only the ones expiring few days before delivery, in analogy to quarterlyand monthly contracts. Since there is not enough liquidity in the market in order toextract information on the IV surface from traded market quotes, we consider thesettlement prices, that are available for a sufficiently large range of maturities andstrikes: though settlement prices do not represent trades that really take place in themarket, they contain information on market expectations. We observe the marketfor a representative day: Monday, March 5, 2018. For each option, we consider thestrike prices in the range 90%–110% of the underlying current price (as, for exam-ple, in [3, 4]). At this date the listed options with available settlement prices arethe ones written on five monthly (Apr/18, May/18, Jun/18, Jul/18, Aug/18), sixquarterly (Q2/18, Q3/18, Q4/18, Q1/19, Q2/19, Q3/19) and three yearly (Cal-19,Cal-20, Cal-21) futures.Before moving on with the calibration procedure, let us comment on the empir-ical market-implied volatilities first, which we have calculated by inverting Black’sformula. They are displayed in Figure 4. We observe a very pronounced smile forthe contracts with forthcoming delivery, i.e. the contracts with delivery in April 18and in the second quarter. Furthermore we observe that the farther the beginningof the delivery period is in the future, the less pronounced is the smile: for thenext beginning periods in May, June and July the smile is present but less and lesspronounced, and already in August and the Q3 contract there is still a tendencyto smile, but not pronounced at all. Generally we observe a forward skew (i.e.higher IVs for out-of-the-money calls). This can be interpreted as a “risk premium”paid by option buyers for securing supply. Remember however, that we consideronly strikes in a range of about 90% − ( u, T , T ) defined in Equation (8), also referred to as the Samuelson effect; thesecond factor, modeled by Γ ( u, T , T ) as in (9), captures seasonal behavior depend-ing on both the delivery period and the no-overlapping-arbitrage condition (25).From the observation of the market implied volatilities, we have some expecta-tions on what the estimated parameters of our futures model should be. First of all,we should get better estimates incorporating the Samuelson-factor, as we do observea different smile behavior close to delivery and far away from delivery. We observe avery pronounced smile shortly before delivery, and a less pronounced smile far awayfrom delivery. Consequently, the estimated NIG parameters of the Samuelson factorshould be such that the resulting distribution is far away from a Normal distribu-tion. Furthermore, the estimated NIG parameters of the seasonal factor should besuch that they are closer to a normal distribution. We perform the calibration procedure described in Section 3.4.3, first for the one-factor model derived by (7) by setting the first coefficient to 0, i.e. Γ ( u, T , T ) ≡ A = 10 and computed by an adaptive Simpson quadrature rule alreadyimplemented in MATLAB (as in [32]), which takes into account the oscillatorybehavior of the integrand (cf. [20]).We find that, in general, the optimization routine falls into local minima. How-ever, by selecting a starting condition that is “sufficiently close” to the observed IVs,the minimization converges (cf. [22] for well-posedness of this kind of problems). Asa by-product from the estimation of the one-factor model, we derive that, under therisk-neutral measure, futures prices are leptokurtic and have significantly positiveskewness. This is reflected also into the shape of empirical IVs, which, as we alreadymentioned, display a forward skew (i.e. higher IVs for out-of-the-money calls). TheIVs of market and model prices are plotted in Figures 5–9. We first address the estimated parameters for the futures model, before we commenton the resulting fit of market implied volatilities. We have estimated the distributionof the driving NIG process. Remember though that we consider centered L´evy16 .0080.010.0120.0140.0160.0180.02 I m p li ed v o l a t ili t y Apr/18: 32.15 Eur/MWh
29 30 31 32 33 34 35
Strike marketBlack1 factor2 factors I m p li ed v o l a t ili t y May/18: 29.31 Eur/MWh
26 27 28 29 30 31 32
Strike marketBlack1 factor2 factors
Figure 1: Selection of monthly delivery periods. Implied volatility for the Black, one-factor and two-factor model compared to the empirical implied volatilities of optionslisted at March 5, 2018 (the corresponding underlying current price is indicatedabove each plot).processes, thus µ = 0 and δ = 1 due to the scaling property of the NIG distribution(see Section 4.2.2). The estimated parameters are reported in Table 1 and 2 for theone factor and two factor model respectively. α β Γ (1) Γ (2) Γ (3) Γ (4) Γ (5) Γ (6)0.0059 0.0019 0.0464 0.0327 0.0315 0.0311 0.0293 0.0368Γ (7) Γ (8) Γ (9) Γ (10) Γ (11) Γ (12) Γ (13) Γ (14)0.0284 0.0310 0.0304 0.0211 0.0209 0.0271 0.0255 0.0244 Table 1: Calibrated parameters for the one factor model at 2018, March 5th. α β α β γ µ Γ (1) Γ (2) Γ (3) Γ (4)0.1890 0.0586 0.0005 0.0002 0.1656 0.0044 0.0129 0.0054 0.0060 0.0068Γ (5) Γ (6) Γ (7) Γ (8) Γ (9) Γ (10) Γ (11) Γ (12) Γ (13) Γ (14)0.0064 0.0081 0.0066 0.0091 0.0093 0.0055 0.0057 0.0093 0.0084 0.0078 Table 2: Calibrated parameters for the two factor model at 2018, March 5th.In Figure 2 we have plotted the theoretical NIG density together with the corre-sponding normal density. We see that the Samuelson factor is far away from beingGaussian, exhibiting rather fat tails. Also the estimated density of the seasonal com-ponent clearly deviates from the Gaussian shape, but in a less pronounced way. Thisconfirms our expectations and justifies our choice of pure jump L´evy processes. Theseasonal component has a tendency to model more “normal” movements, while theSamuelson component is able to account for rare, bigger movements. This analysis iscomplemented by Table 3, where we have presented the estimation of the momentsof the NIG driving factors. Note that these are the (risk-neutral) moments of thejump drivers J (1) and J (1) only, i.e. without considering the coefficients and forfixed time equal to 1. By accounting for the coefficients and integrating in time, one17ould get the (risk-neutral) moments of the futures prices. -8 -6 -4 -2 0 2 4 6 800.050.10.150.20.250.30.350.4 NIGNormal -80 -60 -40 -20 0 20 40 60 8000.050.10.150.20.250.30.35
NIGNormal
Figure 2: Densities of the NIG factors J (left), J (right) compared to the corre-sponding Gaussian densities (i.e. Gaussian density with same mean and variance). component mean variance skewness excess kurtosis J J Table 3: Moments of the NIG driving factors J (top) and J (bottom).We now start to discuss the remaining estimates, that are reported in Table1 and 2 for the one factor and two factor model respectively. We would like todraw the attention to the estimates Γ i of the overlapping delivery periods that arerestricted through the NOA-condition in (24). For our dataset, this is Apr/18,May/18, Jun/18, as well as Q2/18. The corresponding estimates are illustrated inFigure 3.After addressing the properties of our estimated futures model, we move on todiscussing the fit of volatility smiles. We divide the discussion depending on thetime to maturity – forthcoming delivery periods (Apr/18, Q2/18), delivery periodsthat are not immediately forthcoming (May/18, Jun/18, Jul/18, Q3/18, Q4/18)and delivery periods that are far away (in terms of the number of forthcomingperiods). The implied volatility smile of a selection of the delivery periods is plottedin Figure 1, both for the 1-factor and 2-factor model. For comparison, we have addedthe calibrated volatility resulting from Black’s model, whose values are reported inTable 4. While the selection in Figure 1 includes the two forthcoming months, theplots of all smiles that we consider can be found in the Appendix B together withthe corresponding option prices, see Figures 5 - 9. Here, we have also displayedthe smiling volatility surface of monthly, quarterly and yearly delivery periods (seeFigure 4).For the forthcoming contracts, both the 1-factor and 2-factor model underes-timate the implied volatility at the money, and overestimate it in and out of themoney. Nevertheless, the two factor model does well, and captures especially theshape out of the money. It is our winner here. On the contrary, the one factor18 .0250.030.0350.040.0450.05 Γ XX/YY for overlapping contracts
Apr/18May/18Jun/18Q2/18
Figure 3: Γ for overlapping monthly delivery (solid line) and quarterly delivery(dashed lines). σ (1) σ (2) σ (3) σ (4) σ (5) σ (6) σ (7)0.0156 0.0132 0.0123 0.0121 0.0120 0.0137 0.0109 σ (8) σ (9) σ (10) σ (11) σ (12) σ (13) σ (14)0.0106 0.0106 0.0094 0.0094 0.0106 0.0104 0.0100 Table 4: Calibrated parameters for the purely Gaussian model on 2018, March 5th.model mismatches so much that it is hard to call it satisfactory. Consider now thecontracts that are not immediately forthcoming. Here we find an almost perfect fitof the two factor model, which is very satisfactory; however, the one factor modelseems to perform better than before, too. Finally, for the contracts with deliveryperiod in the far future, both factor models give good results. Nevertheless, as thesmile is quite flat, also the purely Gaussian model seems to be a reasonable choice.As we want to have one model with one single set of parameters, that capturesthe implied volatility features of all contracts, forthcoming and not, we can concludethat the two factor model performs overall very well and is a dignified winner of thisempirical study.
In this paper, we present the theory of arbitrage possibilities when trading in futuresand option contracts with overlapping delivery periods (see also [10, 17, 30]). Inthis setting, we discuss the necessary no-overlapping-arbitrage (NOA) conditionsof a NIG-driven two factor model with Samuelson and seasonal factor. Our mainpurpose is the calibration of this model to power option prices via Fourier transformmethods in an overlapping-arbitrage free way. This leads to an additional NOA19estriction in our optimization problem, see Section 4.2.4. Looking at the marketimplied volatilities at our chosen trading day, we observe a pronounced forward-skewed volatility smile for the forthcoming delivery periods. The smile flattens outwhen time to maturity increases. Thus, it is not surprising that our two factormodel — that accounts exactly for these short-term, long-term variations — is ableto fit the smile very well.
Acknowledgments
T. Vargiolu acknowledges financial support from the research grant BIRD172407-2017 of the University of Padova ”New perspectives in stochastic methods for fi-nance and energy markets”. This work has been initiated, and then revised, whileM. Schmeck was visiting the Department of Mathematics of the University of Padovathanks to the fundings provided by the “ACRI Young Investigator Training Pro-gram” (YITP-EFI2018) and the YITP 2019 connected to the Energy Finance ItaliaConference 2019 (UA.MB.D10-Dipartimento di Statistica e Metode Quantitativi,Universita’ degli Studi di Milano - Bicocca), and while M. Piccirilli visited the Uni-versity of Bielefeld.
A Appendix
A.1 Inversion of the Fourier transform: Carr-Madan approach
We here recall the Carr-Madan approach for recovering the option value from itsFourier transform (see Section 3). If ξ MTt,T , defined in (16), is integrable, by theFourier inversion theorem, the modified time value of the option can be recoveredby inverting again the last equation: z MTt,T ( K ) = 12 π Z + ∞−∞ e − iKv ξ MTt,T ( v ) dv. (27)Now, observe that, by the martingale property of F ( · ),( F ( t ) − K ) F ( t ) >K = E [( F ( T ) − K ) |F t ] F ( t ) >K . Then, we can write from (14) z MTt,T ( K ) = Z + ∞−∞ ( s − K ) q t,T ( s ) ( s>K − s ( t ) >K ) ds where, for any t ∈ [0 , T ], s ( t ) := F ( t ) and q t,T is the density function of F ( T )conditioned up to time t . As a consequence, if the interchange of integrals holds20see Proposition A.1), (16) can be written as ξ MTt,T ( v ) = Z + ∞−∞ e ivK Z + ∞−∞ ( s − K ) q t,T ( s ) ( s>K − s ( t ) >K ) ds dK = Z + ∞−∞ q t,T ( s ) Z ss ( t ) e ivK ( s − K ) dK ds (28)= Z + ∞−∞ q t,T ( s ) ( e ivK siv (cid:12)(cid:12)(cid:12)(cid:12) sK = s ( t ) − Z ss ( t ) K e ivK dK ) ds = Z + ∞−∞ q t,T ( s ) ( − e ivs ( t ) siv + e ivs ( t ) s ( t ) iv − e ivs v + e ivs ( t ) v ) ds = − v Z + ∞−∞ q t,T ( s ) e ivs ds + e ivs ( t ) v − e ivs ( t ) iv Z + ∞−∞ q t,T ( s )( s − s ( t )) ds. Since F is a martingale, Z + ∞−∞ q t,T ( s )( s − s ( t )) ds = E [( F ( T ) − F ( t )) |F t ] = 0and, by definition of characteristic function (see also (6)), Z + ∞−∞ q t,T ( s ) e ivs ds = E h e ivF ( T ) |F t i = e ivs ( t ) Ψ( t, T, v ) , where Ψ( t, T, v ) is the characteristic function of Z ( t, T ). Then, we finally arrive to ξ MTt,T ( v ) = e ivF ( t ) − Ψ( t, T, v ) v . (29)In analogy to the modified option approach, we have that ξ MTt,T has even real partand odd imaginary part and so z MTt,T ( K ) = 1 π Z + ∞ Re (cid:0) e − iKv ξ MTt,T ( v ) (cid:1) dv. (30)Finally, by recalling (15), the option value is computed from (30) by C ( t ; T, K ) = z MTt,T ( K ) + ( F ( t ) − K ) + . In order to apply this formula, we need to justify the interchange of integralsoperated in (28) under the following integrability assumptions on the futures priceprocess.
Proposition A.1 If Z ( T, T , T ) is square-integrable, then ξ t,T ( v ) = Z + ∞−∞ e ivK Z + ∞−∞ ( s − K ) q t,T ( s ) ( s>K − s ( t ) >K ) ds dK = Z + ∞−∞ q t,T ( s ) Z ss ( t ) e ivK ( s − K ) dK ds. roof. It is enough to show that the integral of the absolute value of the integrandwith respect to K , i.e. Z + ∞−∞ | s − K | | s>K − s ( t ) >K | dK can be integrated in s against the density q t,T ( s ). Since this is equal to Z ss ( t ) ( s − K ) dK = 12 ( s − s ( t )) , the integrability with respect to the density is equivalent to the existence of thesecond moment of Z ( t, T, T , T ). (cid:3) A.2 Characteristic function of the two-factor NIG model
As already stated in Section 4.2 for general multifactor additive models, the char-acteristic function of two-factor NIG model is defined as a function of v ∈ R byΨ( t, T, T , T , v ) = E h e ivZ ( t,T,T ,T ) |F t i , where Z ( t, T, T , T ) := F ( T, T , T ) − F ( t, T , T ), and it can be shown to be equalto log Ψ( t, T, T , T , v ) = ψ ( t, T ; y Γ ( · , T , T )) + ψ ( t, T ; y Γ ( T , T )) . (31)First, we need to compute the cumulant function e ψ j ( θ ) of J j ( j = 1 , e ψ L j ( θ ) of L j as follows. From thedefinition of J j we know that e ψ j ( θ ) = e ψ L j ( θ ) − iθ χ j + Z | y |≥ y ν j ( dy ) ! , where e ψ L j ( θ ) = δ (cid:16)q α j − β j − q α j − ( β j + iθ ) (cid:17) + iθm j . Moreover, since L j (1) is a NIG distributed random variable, its expected value is E [ L j (1)] = χ j + Z | y |≥ y ν j ( dy ) ! = m j + δ j β j q α j − β j , so that, by replacing it in the expression above, we find that e ψ j ( θ ) = δ j q α j − β j − q α j − ( β j + iθ ) − iθ β j q α j − β j . (32)22n particular, we observe that J j (1) is a NIG distributed random variable withparameters ( α j , β j , δ j , − δ j β j q α j − β j ). Now, we can compute the two jump componentsof the cumulant function in (31), by inserting the corresponding expressions of e ψ j ( θ )from (32) and Γ j as in (8)–(9) for j = 1 ,
2. The second coefficient can be directlycomputed as ψ ( t, T ; y Γ ( T , T )) = ( T − t ) δ (cid:16)q α − β − q α − ( β + iy Γ ( T , T )) − iy Γ ( T , T ) β p α − β (cid:17) , (33)while the first requires an integration in time: ψ ( t, T ; y Γ ( · , T , T )) = ( T − t ) δ q α − β − δ Z Tt q α − ( β + iy Γ ( u, T , T )) du − iy δ β p α − β Z Tt Γ ( u, T , T ) du. (34)Let us recall that, if ζ := Γ ( T , T ) is positive, by the properties of the NIG distri-bution (see e.g. [5]), ζ · J (1) is a NIG distributed random variable with parameters( α /ζ, β /ζ, δ ζ, − δ ζβ √ α − β ). Consequently, it is easy to see that we can assumewithout loss of generality that δ = 1, so that ψ ( t, T ; y Γ ( T , T )) = ( T − t ) (cid:16)q α − β − q α − ( β + iy Γ ( T , T )) − iy Γ ( T , T ) β p α − β (cid:17) , (35)For the same reason we can assume that δ = 1. By replacing Γ ( u, T , T ) = e µu γ ( e − µT − e − µT ) µ ( T − T ) := e µu e Γ ( T , T ) in (34) and integrating, we get ψ ( t, T ; y Γ ( · , T , T )) = ( T − t ) q α − β − η ( c ( y, T )) + η ( c ( y, t )) − iy e Γ ( T , T ) β ( e µT − e µt ) µ p α − β , (36)where η ( w ) := 1 µ (cid:18)q α + w − iβ arcsinh wα − q α − β log 2 α (cid:16) α − iβ w + p α − β p α + w (cid:17) ( w + iβ )( α − β ) , (37) c ( v, u ) := v e Γ ( T , T ) e µu − iβ . (38)23inally, the cumulant function of Z is explicitly given bylog Ψ( t, T, T , T , y ) = ψ ( t, T ; y Γ ( · , T , T )) + ψ ( t, T ; y Γ ( T , T ))= ( T − t ) (q α − β + q α − β − iy e Γ ( T , T ) β µ p α − β ( e µT − e µt ) T − t + Γ ( T , T ) β p α − β ! − η ( c ( y, T )) − η ( c ( y, t )) T − t − q α − ( β + iy Γ ( T , T )) ) . with η ( · ) and c ( · , · ) as in (37)–(38). B Figures
90 92 94 0.008 96 Q2/18 98 Q3/18 1000.01 102Q4/18 104Q1/19 1060.012 Q2/19 108110Q3/190.0140.016 90 92 94 96 Cal-190.01 98 1001020.0105 Cal-20 1041061080.011 110Cal-210.0115
Figure 4: Surface plot of implied volatilites for montly, quarterly and yearly deliveryperiods. Empirical values (dotted, blue) and two factor model (red).25 .0080.010.0120.0140.0160.0180.02 I m p li ed v o l a t ili t y Apr/18: 32.15 Eur/MWh
29 30 31 32 33 34 35
Strike marketBlack1 factor2 factors O p t i on p r i c e Apr/18: 32.15 Eur/MWh
29 30 31 32 33 34 35
Strike marketBlack1 factor2 factors I m p li ed v o l a t ili t y May/18: 29.31 Eur/MWh
26 27 28 29 30 31 32
Strike marketBlack1 factor2 factors O p t i on p r i c e May/18: 29.31 Eur/MWh
26 27 28 29 30 31 32
Strike marketBlack1 factor2 factors I m p li ed v o l a t ili t y Jun/18: 31.44 Eur/MWh
28 29 30 31 32 33 34 35
Strike marketBlack1 factor2 factors O p t i on p r i c e Jun/18: 31.44 Eur/MWh
28 29 30 31 32 33 34 35
Strike marketBlack1 factor2 factors
Figure 5: Implied volatility for the Black, one-factor and two-factor model comparedto the empirical implied volatilities of all the options listed at March 5, 2018 (left).The corresponding underlying current price is indicated above each plot. On theright hand side the corresponding prices are shown. Monthly delivery periods fromApril 18 - June 18. 26 .0080.010.0120.0140.0160.0180.02 I m p li ed v o l a t ili t y Jul/18: 32.33 Eur/MWh
29 30 31 32 33 34 35 36
Strike marketBlack1 factor2 factors O p t i on p r i c e Jul/18: 32.33 Eur/MWh
29 30 31 32 33 34 35 36
Strike marketBlack1 factor2 factors I m p li ed v o l a t ili t y Aug/18: 31.19 Eur/MWh
28 29 30 31 32 33 34
Strike marketBlack1 factor2 factors O p t i on p r i c e Aug/18: 31.19 Eur/MWh
28 29 30 31 32 33 34
Strike marketBlack1 factor2 factors
Figure 6: Implied volatility for the Black, one-factor and two-factor model comparedto the empirical implied volatilities of all the options listed at March 5, 2018 (left).The corresponding underlying current price is indicated above each plot. On theright hand side the corresponding prices are shown. Monthly delivery periods, Julyand August 2018. 27 .0080.010.0120.0140.0160.0180.02 I m p li ed v o l a t ili t y Q2/18: 30.95 Eur/MWh
28 29 30 31 32 33 34
Strike marketBlack1 factor2 factors O p t i on p r i c e Q2/18: 30.95 Eur/MWh
28 29 30 31 32 33 34
Strike marketBlack1 factor2 factors I m p li ed v o l a t ili t y Q3/18: 32.68 Eur/MWh
29 30 31 32 33 34 35 36
Strike marketBlack1 factor2 factors O p t i on p r i c e Q3/18: 32.68 Eur/MWh
29 30 31 32 33 34 35 36
Strike marketBlack1 factor2 factors I m p li ed v o l a t ili t y Q4/18: 38.03 Eur/MWh
34 35 36 37 38 39 40 41 42
Strike marketBlack1 factor2 factors O p t i on p r i c e Q4/18: 38.03 Eur/MWh
34 35 36 37 38 39 40 41 42
Strike marketBlack1 factor2 factors
Figure 7: Implied volatility for the Black, one-factor and two-factor model comparedto the empirical implied volatilities of all the options listed at March 5, 2018 (left).The corresponding underlying current price is indicated above each plot. On theright hand side the corresponding prices are shown. Quarterly delivery periods in2018. 28 .0080.010.0120.0140.0160.0180.02 I m p li ed v o l a t ili t y Q1/19: 38.05 Eur/MWh
34 35 36 37 38 39 40 41 42
Strike marketBlack1 factor2 factors O p t i on p r i c e Q1/19: 38.05 Eur/MWh
34 35 36 37 38 39 40 41 42
Strike marketBlack1 factor2 factors I m p li ed v o l a t ili t y Q2/19: 30.04 Eur/MWh
27 28 29 30 31 32 33
Strike marketBlack1 factor2 factors O p t i on p r i c e Q2/19: 30.04 Eur/MWh
27 28 29 30 31 32 33
Strike marketBlack1 factor2 factors I m p li ed v o l a t ili t y Q3/19: 30.27 Eur/MWh
27 28 29 30 31 32 33
Strike marketBlack1 factor2 factors O p t i on p r i c e Q3/19: 30.27 Eur/MWh
27 28 29 30 31 32 33
Strike marketBlack1 factor2 factors
Figure 8: Implied volatility for the Black, one-factor and two-factor model comparedto the empirical implied volatilities of all the options listed at March 5, 2018 (left).The corresponding underlying current price is indicated above each plot. On theright hand side the corresponding prices are shown. Quarterly delivery periods in2019. 29 .0080.010.0120.0140.0160.0180.02 I m p li ed v o l a t ili t y Cal-19: 33.9 Eur/MWh
31 32 33 34 35 36 37
Strike marketBlack1 factor2 factors O p t i on p r i c e Cal-19: 33.9 Eur/MWh
31 32 33 34 35 36 37
Strike marketBlack1 factor2 factors I m p li ed v o l a t ili t y Cal-20: 33.61 Eur/MWh
30 31 32 33 34 35 36 37
Strike marketBlack1 factor2 factors O p t i on p r i c e Cal-20: 33.61 Eur/MWh
30 31 32 33 34 35 36 37
Strike marketBlack1 factor2 factors I m p li ed v o l a t ili t y Cal-21: 33.99 Eur/MWh
31 32 33 34 35 36 37
Strike marketBlack1 factor2 factors O p t i on p r i c e Cal-21: 33.99 Eur/MWh
31 32 33 34 35 36 37
Strike marketBlack1 factor2 factors
Figure 9: Implied volatility for the Black, one-factor and two-factor model comparedto the empirical implied volatilities of all the options listed at March 5, 2018 (left).The corresponding underlying current price is indicated above each plot. On theright hand side the corresponding prices are shown. Yearly delivery periods.30 eferences [1] M. Abramowitz and I. A. Stegun.
Handbook of mathematical functions withformulas, graphs, and mathematical tables , volume 55 of
National Bureau ofStandards Applied Mathematics Series . For sale by the Superintendent of Doc-uments, U.S. Government Printing Office, Washington, D.C., 1964.[2] A. Andresen, S. Koekebakker, and S. Westgaard. Modeling electricity forwardprices using the multivariate normal inverse gaussian distribution.
The Journalof Energy Markets , 3(3):3, 2010.[3] J. C. Arismendi, J. Back, M. Prokopczuk, R. Paschke, and M. Rudolf. Sea-sonal Stochastic Volatility: Implications for the pricing of commodity options.
Journal of Banking & Finance , 66(C):53–65, 2016.[4] J. Back, M. Prokopczuk, and M. Rudolf. Seasonality and the valuation ofcommodity options.
Journal of Banking & Finance , 37(2):273–290, 2013.[5] O. E. Barndorff-Nielsen. Normal Inverse Gaussian distributions and stochasticvolatility modelling.
Scand. J. Statist. , 24(1):1–13, 1997.[6] O. E. Barndorff-Nielsen. Processes of Normal Inverse Gaussian type.
Financeand Stochastics , 2(1):41–68, 1997.[7] F. E. Benth and J. altyt Benth. The normal inverse gaussian distribution andspot price modelling in energy markets.
International Journal of Theoreticaland Applied Finance , 07(02):177–192, 2004.[8] F. E. Benth, J. Kallsen, and T. Meyer-Brandis. A non-Gaussian Ornstein-Uhlenbeck process for electricity spot price modeling and derivatives pricing.
Appl. Math. Finance , 14(2):153–169, 2007.[9] F. E. Benth and R. Kufakunesu. Pricing of exotic energy derivatives basedon arithmetic spot models.
International Journal of Theoretical and AppliedFinance , 12(04):491–506, 2009.[10] F. E. Benth, M. Piccirilli, and T. Vargiolu. Mean-reverting additive energy for-ward curves in a Heath-Jarrow-Morton framework.
Mathematics and FinancialEconomics. In press. , 2019.[11] F. E. Benth and M. D. Schmeck. Pricing and hedging options in energy marketsusing Black-76.
Journal of Energy Markets , 7(2):35–69, 2014.[12] F. E. Benth and M. D. Schmeck.
Pricing Futures and Options in ElectricityMarkets , pages 233–260. Springer Berlin Heidelberg, Berlin, Heidelberg, 2014.[13] F. E. Benth, J. ˇSaltyt˙e Benth, and S. Koekebakker.
Stochastic Modelling ofElectricity and Related Markets . Advanced Series on Statistical Science & Ap-plied Probability, 11. World Scientific Publishing Co. Pte. Ltd., Hackensack,NJ, 2008. 3114] N. H. Bingham and R. Kiesel.
Risk-neutral valuation . Springer Finance.Springer-Verlag London, Ltd., London, second edition, 2004. Pricing and hedg-ing of financial derivatives.[15] O. H. Birkelund, E. Haugom, P. Molnr, M. Opdal, and S. Westgaard. A com-parison of implied and realized volatility in the Nordic power forward market.
Energy Economics , 48:288 – 294, 2015.[16] F. Black. The pricing of commodity contracts.
Journal of Financial Economics ,3(1):167 – 179, 1976.[17] R. B¨orger. Energy-related commodity futures: statistics, models and deriva-tives. 2007. Ph.D. Dissertation.[18] S. Borovkova and M. D. Schmeck. Electricity price modeling with stochastictime change.
Energy Economics , 63:51–65, 2017.[19] D. Brigo and F. Mercurio.
Interest rate models—theory and practice . SpringerFinance. Springer-Verlag, Berlin, second edition, 2006. With smile, inflationand credit.[20] P. Carr and D. B. Madan. Option valuation using the fast fourier transform.
Journal of Computational Finance , 2:61–73, 1999.[21] R. Cont and P. Tankov.
Financial Modelling with Jump Processes . Chapman &Hall/CRC Financial Mathematics Series. Chapman & Hall/CRC, Boca Raton,FL, 2004.[22] R. Cont and P. Tankov. Nonparametric calibration of jump-diffusion optionpricing models.
Journal of Computational Finance , 7:1–49, 2004.[23] E. Edoli, M. Gallana, and T. Vargiolu. Optimal intra-day power trading witha Gaussian additive process.
Journal of Energy Markets , 10(4):23–42, 2017.[24] EEX. Trading Volume in Power Derivatives, June 2018. .[25] V. Fanelli, L. Maddalena, and S. Musti. Modelling electricity futures pricesusing seasonal path-dependent volatility.
Applied Energy , 173:92 – 102, 2016.[26] V. Fanelli and M. D. Schmeck. On the seasonality in the implied volatility ofelectricity options.
Quantitative Finance , 19(8):1321–1337, 2019.[27] E. Jaeck and D. Lautier. Volatility in electricity derivative markets: TheSamuelson effect revisited.
Energy Economics , 59:300 – 313, 2016.[28] R. Kiesel and F. Paraschiv. Econometric analysis of 15-minute intraday elec-tricity prices.
Energy Economics , 64:77 – 90, 2017.[29] R. Kiesel, G. Schindlmayr, and R. H. Brger. A two-factor model for the elec-tricity forward market.
Quantitative Finance , 9(3):279–287, 2009.3230] L. Latini, M. Piccirilli, and T. Vargiolu. Mean-reverting no-arbitrage additivemodels for forward curves in energy markets.
Energy Economics , 79:157–170,2019.[31] E. Nastasi, A. Pallavicini, and G. Sartorelli. Smile modelling in commoditymarkets. Preprint, 2018.[32] N. K. Nomikos and O. A. Soldatos. Analysis of model implied volatility forjump diffusion models: Empirical evidence from the Nordpool market.
EnergyEconomics , 32(2):302–312, 2010.[33] M. D. Schmeck. Pricing options on forwards in energy markets: the role of meanreversion’s speed.
International Journal of Theoretical and Applied Finance ,19(8):1650053, 2016.[34] L. Schneider and B. Tavin. From the samuelson volatility effect to a samuelsoncorrelation effect: An analysis of crude oil calendar spread options.
Journal ofBanking & Finance , 2016.[35] L. Schneider and B. Tavin. The Samuelson Effect and Seasonal StochasticVolatility in Agricultural Futures Markets. Preprint, 2018.[36] S. J. Taylor.