DIMENSIONAL MODELING AND SIMULATION
ANDREA BARTH AND FRED ESPEN BENTH
Abstract. In this paper an infinite dimensional approach to model energy forward mar- kets is introduced. Similar to the Heath-Jarrow-Morton framework in interest rate mod- eling, a first order hyperbolic stochastic partial differential equation models the dynamics of the forward price curves. These equations are analyzed, and in particular regularity and no-arbitrage conditions in the general situation of stochastic partial differential equa- tions driven by an infinite dimensional martingale process are studied. Both arithmetic and geometric forward price dynamics are studied, as well as accounting for the delivery period of electricity forward contracts. A stable and convergent numerical approximation in form of a finite element method for hyperbolic stochastic partial differential equations is introduced and applied to some examples with relevance to energy markets.
1. Introduction
The Heath-Jarrow-Morton (HJM) framework was originally introduced to model prod- ucts in fixed income markets. Rather than to base the analysis of price processes on the spot rate(s), Heath, Jarrow and Morton [34] proposed to model the forward rates directly.
This approach, however, is not limited to interest-rate derivatives: forward markets rep- resent a considerable portion of the trading on energy and commodities. Therefore, it is natural to apply the HJM approach to such settings (see Benth et al. [16], Koekebakker and Ollmar [37], etc). To introduce models for the forward price dynamics directly rather than via spot prices circumvents the difficulty of linking spot and forward prices. In a market like electricity, the spot-forward relationship is far from straightforward.
The HJM modeling paradigm is naturally embedded into an infinite dimensional frame- work. In fixed-income theory, there is already an extensively developed theory in this direction (see Filipovi´c [29], Carmona and Tehranchi [20], and the references therein). The purpose of this paper is to analyze a class of HJM models for forwards in energy markets within an infinite dimensional framework. Further, we develop convergent approximations
Date: November 15, 2013.
Key words and phrases. Energy markets, forward price, random fields, infinite dimensional analysis, stochastic partial differential equation.
It is a pleasure to thank Claude J. Gittelson, Annika Lang, Santiago Moreno Bromberg, J¨urgen Pot- thoff and Oleg Reichmann for all the fruitful discussions and helpful comments. We are grateful for the comments and critics of two anonymuous referees. Fred Espen Benth acknowledges financial support from the Norwegian Research Council under the eVita project 205328 “Energy Markets: modeling, optimization and simulation” (Emmos).
1
of the resulting infinite dimensional equations with the purpose of simulating sample paths or functionals of the solution, such as moments, via Monte Carlo methods.
Energy markets have several distinct features that qualify them for a separate study within the HJM approach. First, in markets such as gas and electricity, contracts do not deliver at fixed maturities. By the very nature of these commodities, delivery must take place over an agreed period of time. The holder of such a commodity receives a flow of gas/electricity for some prescribed period. This introduces a new parameter: the length- of-delivery. The presence of such a crucial element further complexifies the existing theory for forward rates in money markets, which only present a fixed maturity as parameter. We shall dubflow forwards those contracts that deliver over the span of some time period, and present an infinite dimensional framework for both such forwards and those delivering at a fixed maturity time. As it turns out, the ”classical” forward contracts with fixed time of delivery of the underlying is an important modeling tool for the price dynamics of flow forwards.
A precise model of electricity markets requires an infinite number of stochastic pro- cesses driving the forward curve. In interest rate theory, as suggested in Carmona and Tehranchi [20], three driving noise sources (Wiener processes) and their linear combina- tions cover more than 99% of the variation in the data. In these markets there seem to be a strong connection among various maturities, which may explain why low-dimensional noise processes are sufficient to model the yield curve. In contrast, the picture in electricity markets is very different. In Koekebakker and Ollmar [37] the authors show, using princi- pal component analysis, that more than 12 factors are necessary to cover only 95% of the variation in the data. This suggests that there is a sizable number of noise sources in the forward market, where each contract has a large degree of idiosyncratic risk associated to it. Hence, an infinite dimensional approach to modeling the noise in the forward curves of energy markets is in order.
There are other arguments in favor of using infinite dimensional processes as the driving noise of the forward curves in energy markets, in particular markets where the underlying commodity is not storable (such as electricity). As it is explained in Benth and Meyer- Brandis [15], in such markets forward-looking information may be relevant to the forward markets, but not to the spot market. In electricity markets, say, there is a vast stream of information that impacts forward prices, including weather prediction, power plant shut- downs, production plans, reservoir situations, etc. This translates into a large number of noise sources that affect the forward curve at different maturities, and therefore create a high-dimensional noise pattern. Using a completely different approach to forward price modeling, namely the certainty equivalence principle, Benth, Cartea and Kiesel [11] show that no-arbitrage models can only be obtained by introducing idiosyncratic risk to each traded forward. If one accepts such a pricing principle, then one is forced to deal with as many sources of noise as maturity times, which would as well support the choice of infinitely many sources of noise in our model.
In Benthet al. [16] and Frestadet al. [32], it is shown that the (log-)returns of flow for- ward prices in electricity markets are non-normally distributed, independent of the delivery period. The normal inverse Gaussian (NIG) distribution has been used successfully to fit
the observed return data. This points towards using L´evy processes to model the noise driving the forward curve, which, together with the above observation of high dimensional noise in these markets, suggests the use of an infinite dimensional L´evy processes as the driving noise source in the forward curve dynamics.
In this paper, we separate into two classes of forward models, additive and multiplicative.
Inspired by the existing models in the literature and the “statistical” modeling view, we take simple, infinite dimensional dynamics for the forward curve evolution specified under the Musiela parametrization. In fact, our infinite dimensional dynamics is specified in terms of a drift and volatility term structure, as well as a noise modeling the dependency structure over delivery times (spatial dependency, that is). This resembles the “statistical”
modeling view of asset prices, first advocated by Eberlein and Keller [26] for stock prices using an exponential L´evy process and recently applied to electricity forward pricing by Andresenet al.[1]. As forward prices are positive, we also treat exponential models in this paper, which we refer to as multiplicative. The infinite dimensional process driving the forward curve will be explicitly representable, and not derived as a solution to a stochastic differential equation (in infinite dimensions). For a general introduction to Hilbert space valued stochastic differential equations with Brownian noise we refer the reader to Da Prato and Zabczyk [22], Chow [21], Pr´evˆot and R¨ockner [41] and to Peszat and Zabczyk [40] for the case of L´evy processes.
Additive forward price models have been used with great success in energy markets, see for example Benth et al. [12]. Similar approaches in interest-rate theory are presented in Ekeland and Taflin [27], Peszat and Zabczyk [40], Filipovi´c [29] and Carmona and Tehranchi [20]. Here we extend these approaches to energy markets. We cover as well sev- eral other models, for example the multiplicative ones by Audetet al.[2] and Schwartz [47].
Particularly, the model presented in Audetet al.[2] is further investigated as it is the basis for some numerical examples. We also note that the ambit fields suggested as a model for energy forwards in Barndorff-Nielsen et al. [3] is nicely corresponding to our framework.
Audet et al.[2] applied their model to price European options on electricity forwards. Re- cently, this model was generalized to L´evy random fields by Hepperger [35], who develops efficient numerical schemes for solving the integro-partial differential equations describing the price of options on forwards.
Flow forwards are linked to forwards with fixed delivery time through a simple integral.
Both classes of contracts must satisfy a martingale condition in order to avoid an arbitrage dynamics for the prices. We derive this explicitly for both the additive and multiplicative models in the case of fixed delivery forwards, and for the additive case for flow forwards.
Furthermore, flow forwards must satisfy an additional consistency condition with respect to the delivery period in order to prevent arbitrage opportunities between contracts with overlapping delivery. We characterize the volatility term structures such that the flow forward dynamics satisfies the arbitrage consistency condition. A particular example is discussed in detail, leading back to a term structure for the volatility incorporating the Samuelson effect for forwards with a delivery period.
One of the main contributions of the paper is to provide a numerical method for simulat- ing the forward curves. We face the problem of numerically solving a first-order hyperbolic
partial differential equation, which we approximate via a Finite Element method. This method provides great flexibility in the choice of the nodal points of the approximation, which will be the various times to maturity of the forward contracts. For first-order hy- perbolic problems with non-smooth data, a standard Galerkin approximation will lead to instabilities, thus we chose a Petrov–Galerkin method (SUPG method). Convergence results for such schemes for deterministic equations and a thorough introduction to Fi- nite Element methods for deterministic first-order hyperbolic equations can be found, e.g., in [43, 36, 48, 24]. In the stochastic case, L2-convergence was proved in Barth [5]. To calculate expectations of functionals of the forward prices, such as moments or payoffs, a fast multilevel Monte Carlo algorithm was introduced in Barth and Lang [9] for hyper- bolic stochastic partial differential equations. There the authors showed convergence of the multilevel Monte Carlo algorithm and its advantages over a singlelevel Monte Carlo method. The method discussed and applied here can also be used in fixed-income theory for simulating forward rates. We have included several numerical examples with relevance to energy forward pricing. In particular, we look at a forward price model which resembles the one suggested by Audet et al.[2], but where we also consider L´evy-noise driven fields.
We choose L´evy fields which are of normal inverse Gaussian type, relating our cases back to the statistical analysis and modeling in Frestad et al. [32] and Andresen et al. [1]. In particular, the implied spot price dynamics is simulated from the fields and we discuss the effect of spatial correlation on the paths.
The results are presented as follows: We describe our infinite dimensional framework for forward contracts delivering at a fixed time of maturity in Section 2, and present regularity results on the smoothness of the model. We analyze the features of additive and multiplicative models. Further, we briefly discuss an extension to these popular models.
In Section 3 we develop an infinite dimensional model directly for flow forwards. This approach satisfies the no-arbitrage condition presented in Benth et al. [16]. Finally, in Section 4 we develop a numerical scheme and present some simulation results.
2. Infinite dimensional modeling of forward prices
Here we present a general infinite dimensional model for the forward price f(t, T) at time t ≥ 0 of a contract delivering the underlying commodity at time T ≥ t, where t, T ∈ [0, τ] and τ < +∞. Much of the theory developed for such price dynamics can be transferred to flow forwards, which we shall study in more detail in the next section. We provide conditions for an arbitrage-free price dynamics for both additive and multiplicative models.
As it is standard in interest rate theory, we consider the so-calledMusiela parametrization (see [20]), where we set x =T −t, the time to maturity, and consider the price g(t, x) of a forward contract at time t ≥ 0 with x ≥ 0 time left to maturity. As we see, g(t, x) is a function of “time and space” (t, x) ∈ [0, τ]×[0, τ], and it is therefore natural to look for a stochastic dynamics in time t 7→ g(t,·), where g(t,·) takes values in some appropriate function space. We place ourselves within the theory of Hilbert space valued stochastic processes when modeling g, which is done under the risk-neutral measure. Thus, we must
choose dynamics for g such that the process
t 7→f(t, T) = g(t, T −t) (t∈[0, T])
is a martingale for every maturityT. This leads tomartingale conditionson the parameters in the specification ofg(see [20], for example, for the similar situation in interest-rate theory and [23] on the fundamental theorem of asset pricing).
It is not obvious whether the model for g should be additive or multiplicative, that is, if we should model the dynamics of g directly rather than ln(g). In Benth et al. [12]
additive spot-price models were introduced and justified for energy markets. This leads to additive forward (and, specifically flow-forward) price dynamics. Furthermore, in Bern- hardt, Kl¨uppelberg and Meyer-Brandis [17], a statistical analysis of Singapore electricity spot data showed evidence in favor of additive models, in particular continuous-time au- toregressive moving-average processes. In view of these results, one may be tempted to model g directly. On the other hand, forward prices are naturally positive, thus it may be advantageous to model ln(g) instead. This would also be in line with some classical examples of forward price models like the one by Audet et al.[2]. In this paper we discuss both alternatives.
In the mathematical specifications of our model we follow the general set-up presented by Peszat and Zabczyk [40] (in particular Chapters 8 and 9). Let (Ω,F,Q) be a probability space equipped with a filtration F = (Ft, t ≥ 0), satisfying “the usual condition”. The process M is a (possibly discontinuous) square-integrable martingale taking values in a separable Hilbert space (U,(·,·)U). This class of processes is denoted by M2(U). Observe that we denote the probability measure on (Ω,F) by Q, indicating that we model the forward curve dynamics directly in the risk-neutral world.
The process M is used as an integrator in the random evolution of the forward curve, and for this purpose we restrict ourselves to the following subclass ofM2(U) of martingales with bounded stochastic covariance:
M2b(U) :={M ∈ M2(U) :∃Q∈L+1(U) s. t. ∀t≥s≥0, hhM, Miit− hhM, Miis≤(t−s)Q}, where L1+(U) is the space of nuclear, symmetric, and positive definite linear operators on U. For every Q ∈ L+1(U), there exists an orthonormal basis (en, n ∈ N) of U consisting of the eigenvectors of Q. This leads to the spectral representation Qen =γnen, where γn
is the eigenvalue of Q corresponding to en. We use the notation H =Q1/2(U), where the (Hilbert-Schmidt) operator Q1/2 is the pseudo square root of Q. Note that TrQ < +∞
since Q is nuclear.
We define the following stochastic partial differential equation with additive noise, i.e.
a Hilbert-space-valued stochastic differential equation for t∈[0, τ] (2.1) dX(t) = (AX(t) +a(t))dt+b(t)dM(t)
with initial condition X(0) =X0 ∈L2(Ω, H), where H denotes a separable Hilbert space, not necessarily equal to U. Later, we shall introduce conditions on the Hilbert spaces which are natural in our context. The stochastic differential equation is defined on a finite time interval 0≤t ≤ τ, where τ <+∞. We assume further, that the linear (differential) operator A generates a C0-semigroup S of contractions on H. When we link X to g, we
show subsequently that A is a first order differential operator withS being the right-shift semigroup. We suppose that a is a mapping from [0, τ] into H which has Q-a.s. Bochner- integrable trajectories in the finite interval [0, τ]. The process b : [0, τ] → L(H, H) is for every t ∈[0, τ] a linear operator from H to H.
Further, we denote by k · kLHS(H,H) the Hilbert–Schmidt norm on the space of Hilbert–
Schmidt operators fromHtoH. Equation (2.1) is only well defined ifX0isF0-measurable.
It is well known (see Peszat and Zabczyk [40], Filipovic, Tappe and Teichmann [30] and Marinelli [39]) that (under appropriate assumptions on the coefficients) a unique mild solutionX inH to Equation (2.1) exists, and that the predictable processX : [0, τ]×Ω→ H is given by
(2.2) X(t) = S(t)X0+ Z t
0
S(t−s)a(s)ds+ Z t
0
S(t−s)b(s)dM(s). Next, we introduce the forward price dynamics based on the process X.
As we discussed earlier, both additive and multiplicative models are of interest in energy forward markets. We introduce both classes here, based on the general infinite dimensional stochastic process X(t) presented in Equation (2.1) (or equivalently, Equation (2.2)).
We remark that the model of X in Equation (2.1) is not the most general one from the mathematical point of view, with no dependency on the states in the drift a and volatilityb. In fact, we have taken the “statistical” modeling view here where one specifies directly the drift and volatility, as well as the noise structure. The crucial points from a statistical viewpoint is to capture the dependency between contracts as well as the volatility term structure. Typically, the existing models in energy can be traced back to the form Equation (2.2) (see, e.g., Barndorff-Nielsen et al. [3]). For example, this is the case with the model of Audet et al. [2] and Andresen et al. [1]. The former can be represented as the exponential of a Wiener field modulated by a volatility structure, while the latter corresponds to the exponential of a L´evy field of normal inverse Gaussian type. We shall come back to this with more details below.
So far, we have formulatedX as a stochastic process with values in some abstract Hilbert space. Since our interest is to model the dynamics of curves t 7→ g(t, x), where g(t, x) is the forward price at time t of a contract with time-to-maturity x≥ 0, we must introduce a function space H. A minimum requirement (see, e.g., Carmona and Tehranchi [20] and Filipovi´c [29]) is the following assumption:
Assumption 2.1. H is a space of real-valued continuous functions in x∈[0, τ]⊂R+, on which the evaluation functional δx is continuous, that is, belongs to the dual H∗.
Further, the Musiela parametrization tells us that f(t, T) = g(t, T −t), where f(t, T) is the forward price at time t for a contract maturing at time T. Thus, looking at the differential off introduces naturally a derivative ofgwith respect toxin the drift. Thus, we only consider an operatorAbeing the first-order differential operatorA=∂/∂x. Obviously, under Assumption 2.1 withH being a space of continuous functions on a subset Dof R+, the domain ofA,D(A), is in general not equal toH. However,D(A) is a Hilbert space with
norm kfkD(A) = (kfk2H +kAfk2H)12. Generated by A is the shift operator S(t) which acts on elements f ∈H as S(t)f(x) = f(x+t), for x∈D. We have,the following assumption:
Assumption 2.2. The operator A is the generator of the C0-semigroup of contractions S(t), t≥0.
Remark 2.3. Assumption 2.2 is quite restrictive, but in the previously described setting sufficient. Equation (2.2) is also defined for the more general classes of semigroups.
In the rest of this Section we suppose that the separable Hilbert spaceHand the operator A satisfy Assumptions 2.1 and 2.2.
Next, we give examples of Hilbert spaces H which satisfy our conditions. The Hilbert space H is a function space over D ⊂ R+. If we assume that H = L2(D), then for all Sobolev spaces Hs(D) with s > 12 the evaluation functional δ is an element of the dual space (Hs(D))∗, also called H−s(D). For integers s the norm in Hs(D) is given by:
kfk2Hs(D) = Z
D
f2 + Xs
i=1
dif dxi
2! dx .
The right shift semigroup S is well defined on Hs(D), since Hs(R) and H−s are invariant under the shift translation. Here,H−s denotes the subset ofHs(R) containing all functions f with support (−∞,0], such that f(t) = 0 for all t >0.
Note that for s > 12 all functions in Hs(D) converge to zero at τ → ∞. In an additive model this would mean that in the far end of the market, forward prices are essentially zero, which seems like an unreasonable property. Most spot models considered for energy markets predict forward prices being constant in the far end of the forward curve, however, not being equal to zero (see e.g. Ornstein-Uhlenbeck models discussed in Benthet al.[16]).
To deal with this issue, Ekeland and Taflin introduced in [27] a special class of Sobolev spaces allowing for non-zero asymptotic behaviour. We define
Us(D) := {h | h= ˆh+ ˆα, ˆh∈Hs(D) and ˆα∈R}, equipped with the norm
khk2Us(D) =khkˆ 2Hs(D)+ ˆα2.
The decompositionh= ˆh+ ˆαis unique, soUs(D) has the Hilbert space structureHs(D)⊕R, and its dual is given by (Us(D))∗ = (Hs(D))∗⊕R. The shift semigroupS can be extended toUs(D) by defining for ˆα∈R, S(t)ˆα= ˆα.
We have the following regularity result for the mild solutionX, given in Equation (2.2):
Lemma 2.4. If EkX0k2H1(D)<∞ and for a constant C < ∞, it holds sup
t∈[0,τ]
(ka(t)k2H1(D)+kb(t)k2LHS(H;H1(D)))≤C ,
then, the unique solution X to Equation (2.2) satisfies supt∈[0,τ]EkX(t)k2H1(D) <∞.
Proof. By the definition of the mild solution we may write EkX(t)k2H1(D) =EkS(t)X0+
Z t 0
S(t−s)a(s)ds+ Z t
0
S(t−s)b(s)dM(s)k2H1(D). With the Itˆo isometry (see [40]) for square integrable (noncontinuous) martingales and since S is a C0-semigroup it follows that
EkX(t)k2H1(D)
≤3(EkX0k2H1(D)+Ek Z t
0
S(t−s)a(s)dsk2H1(D)+Ek Z t
0
S(t−s)b(s)dM(s)k2H1(D))
≤C(τ)(EkX0k2H1(D)+ Z t
0
ka(s)k2H1(D)ds+ Z t
0
kb(s)k2LHS(H,H1(D))ds)
<∞,
where the constant C only depends on τ.
We remark that the regularity result is valid for more general Hilbert spaces and, in particular, also for U1(D) resp. U2(D).
Based on the H-valued stochastic process X we define the additive forward price dy- namics as
(2.3) ga(t, x) = δx(X(t)) = X(t, x). The multiplicative forward dynamics is assumed to follow
(2.4) gm(t, x) = exp(δx(X(t))) = exp(X(t, x)).
We recall that we model the dynamics under the risk-neutral probability directly, which is the common approach also in fixed-income markets. Notice that we obtain the spot price dynamics E(t) in the market by consideringx= 0, that is, E(t) :=Ea(t) =δ0(X(t)) for the additive model, and E(t) := Eg(t) = exp(δ0(X(t))) in the case of a multiplicative forward model. The spot dynamics is of interest, and we shall investigate it numerically in Section 4.
We must impose conditions on the parameters and the operators in the dynamics ofX(t) to guarantee that f(t, T) is a martingale. For the additive model we find the following:
Proposition 2.5. Let f(t, T) = ga(t, T −t) for t ≤ T. Then, under the assumptions of Lemma 2.4, t7→f(t, T) is a martingale if and only if a= 0 almost everywhere.
Proof. Since δx is a linear functional on H, we have ga(t, x) = δx(X(t))
=δxS(t)X(0) +δx Z t
0
S(t−s)a(s)ds+δx Z t
0
S(t−s)b(s)dM(s)
=δxS(t)X(0) + Z t
0
δxS(t−s)a(s)ds+ Z t
0
δxS(t−s)b(s)dM(s). From elementary considerations, we haveδx =δ0S(x). Thus, for every f ∈H we find
δxS(t−s)f =δ0(S(x)S(t−s))f =δ0S(x+t−s)f . Therefore, by lettingx=T −t
f(t, T) = δ0S(T)X(0) + Z t
0
δ0S(T −s)a(s)ds+ Z t
0
δ0S(T −s)b(s)dM(s).
Since δ0S(T −s)b(s) is a linear functional on U, and M is a U-valued martingale, it holds that the process
t 7→
Z t 0
δ0S(T −s)b(s)dM(s),
is a real-valued martingale (see Thm. 8.7(iii) in Peszat and Zabczyk [40]). Hence, t 7→
f(t, T) is a martingale if and only if the drift is zero. The proposition holds.
We now turn our attention to the more involved geometric case. First, by definition, the forward price dynamics is given as
f(t, T) = gm(t, T −t)
= exp δT−tS(t)X(0) + Z t
0
δT−tS(t−s)a(s)ds+ Z t
0
δT−tS(t−s)b(s)dM(s)
= exp δ0S(T)X(0) + Z t
0
δ0S(T −s)a(s)ds+ Z t
0
δ0S(T −s)b(s)dM(s) ,
since δxS(y) =δ0S(x+y). Letting lnf(0, T) =δ0S(T)X(0) and observing that a(t)(x) = δ0S(x)a(t), we find
(2.5) f(t, T) = f(0, T) exp Z t
0
a(s)(T −s)ds+ Z t
0
δ0S(T −s)b(s)dM(s)
.
we introduce the following short-hand notation: define the real-valued, continuous, finite variation process φ as
(2.6) φ(t) =
Z t 0
a(s)(T −s)ds , and the real-valued stochastic process ψ as
(2.7) ψ(t) =
Z t 0
δ0S(T −s)b(s)dM(s).
Note that by Thm. 8.7 (iii) in Peszat and Zabczyk [40], ψ is a martingale. We have the following martingale condition in the geometric case:
Proposition 2.6. The process t 7→f(t, T), for t≤T, is a (local) martingale on R if and only if
dφ(t) =−1
2d[ψ, ψ]c(t)−
e∆ψ(t)−1−∆ψ(t) ,
where [ψ, ψ]c is the continuous part of the bracket process of the martingale ψ (see Prot- ter [42]) and ∆ψ(t) = ψ(t)−ψ(t−) is the jump of ψ at time t.
Proof. We apply Itˆo’s Formula as formulated in Appendix D.1 in Peszat and Zabczyk [40]
to Equation (2.5), which gives f(t, T) = f(0, T) +
Z t 0
f(s−, T) (dφ(s) +dψ(s)) + 1
2 Z t
0
f(s−, T)d[ψ, ψ]c(s)
+X
s≤t
{∆f(s, T)−f(s−, T)∆(φ(s) +ψ(s))} .
Since φ(s) is a continuous process, we have ∆(φ(s) +ψ(s)) = ∆ψ(s). Furthermore,
∆f(s, T) =f(s, T)−f(s−, T) = f(s−, T) e∆ψ(s)−1 . Hence,
f(t, T) = f(0, T) + Z t
0
f(s−, T)dψ(s) + Z t
0
f(s−, T)dφ(s) + 1
2 Z t
0
f(s−, T)d[ψ, ψ]c(s)
+X
s≤t
f(s−, T)
e∆ψ(s)−1−∆ψ(s) .
The integral with respect to ψ above gives rise to a (local) martingale, which completes
the proof.
As an example, consider first the simple case of a real-valued noise processM =W,W being a Brownian motion. Hence, we study the equation
X(t) = S(t)X(0) + Z t
0
S(t−s)a(s)ds+ Z t
0
S(t−s)b(s)dW(s).
Then, φ(t) = Rt
0 S(t−s)a(s)ds is a continuous (deterministic) function, ifa is continuous.
Further, b(s), for s ∈ [0, τ], is the multiplication operator on a Hilbert space H. So we have
ψ(t) = Z t
0
b(s)(T −s)dW(s).
It holds that d[ψ, ψ]c(t) = b2(t)(T −t)dt, and we acquire the condition a(t)(x) =−1
2b2(t)(x),
as expected. Extending this to an infinite-dimensional Q-Wiener processM =W, we may write
ψ(t) = X∞
n=1
pλn
Z t 0
δ0S(T −s)b(s)endWn(s)
by appealing to the series representation of W, given by
(2.8) W(t) =
X∞
n=1
pλnWn(t)en.
Here, Wnare independent, real-valued Brownian motions, and λn are the eigenvalues of Q with eigenvectors en constituting a basis inU. By independence of Wn, we find
[ψ, ψ]c(t) = X∞
n=1
λn
Z t 0
(δ0S(T −s)b(s)en)2ds
= Z t
0
X∞
n=1
|δ0S(T −s)b(s)Q1/2en|2ds
= Z t
0
kδ0S(T −s)b(s)Q1/2k2LHS(U,R)ds ,
wherek · kLHS(U,R) denotes the Hilbert-Schmidt norm for operators from U toR. Thus, we get the martingale condition
a(t)(x) = −1
2kδ0S(x)b(t)Q1/2k2LHS(U,R).
We give an example relating our model to the geometric Brownian motion dynamics for forward contracts proposed by Audet et al. [2].
Example 2.7. Given the forward dynamics
(2.9) df(t, T) = f(t, T)e−α(T−t)σ(T)dBT(t) ∀t ∈[0, T] and T ∈[0, τ],
whereσ is a deterministic volatility function depending on maturity timeT only, and αis a positive constant such that exp(−α(T−t))σ(T)models the Samuelson effect in the forward market. This term may also be related to the spot price dynamics, with α as the speed of mean-reversion (see Benth et al. [16]). The noiseBT(t) is a Gaussian random field, being a Wiener process for each fixed maturity T, with t≤T. For each fixed timet ≤T, it has a covariance kernel given by the function q(T, T′) = exp(−κ|T −T′|), here κ >0 is a range parameter. Indeed, such a covariance structure was partly rejected by Frestad [31] and Andersen et al. [1], who detected additional seasonality effects in the correlation between swap contracts of different maturity in the Nordic electricity market NordPool.
With x = T −t, we define W(t, x) = Bt+x(t). W may be interpreted as a Gaussian process with values in the Hilbert space U = L2(R) and having a covariance kernel given by q(x, y) = exp(−κ|x−y|). Thus, we can express f(t, t+x) by
lnf(t, x+t) = lnf(0, x+t)−1 2
Z t 0
σ2(x+ (t−s) +s)e−2α(x+t−s)ds +
Z t 0
σ(x+ (t−s) +s)e−α(x+t−s)dW(s, x).
We therefore identify the volatility operator b(t)as the multiplication operator on U into a Hilbert space H of functions on R as
δx(b(t)(u)) = σ(x+t)e−αxu(x).
We must impose regularity conditions on σ in order to ensure that the operator b(t) maps into the chosen Hilbert space. For instance, we may choose H =L2(R), and require that y 7→ σ(y) is a bounded measurable function on R. This identifies the model of Audet et al. [2] in our infinite-dimensional framework. We note that the martingale condition also holds, giving a no-arbitrage dynamics for the forward process.
Equation (2.1) describes the risk-neutral dynamics of the log-forward price. In energy markets, forward contracts are applied for hedging purposes, and thus one is interested in the market dynamics of the forward prices. Also, if we want to estimate the parameters in the forward price model, the observed data is under the market probability and not the risk-neutral one. Obviously, if one has option prices on forwards available, one could use these to extract the “implied” parameters. However, sometimes, like in power, the option market may be rather illiquid and it is preferrable to use the forward prices observed in the market for estimation instead. To obtain the dynamics under the real-world measure P, we must change the measure fromQ toP. This is what we discuss in some more detail next:
In the language of energy markets, themarket price of risk, may be modeled as a process with values in the same space as the noise process. If U is a function space over x ∈ D, time-to-maturity, then we can have a market price of risk being a function of x. This is in contrast to the case with finite-dimensional noise, where the market price of risk must have the same dimension as the noise. In situations where one applies for example utility- based methods for pricing, one is forced to enlarge the dimensionality of the driving noise to match the number of forward contracts in the market to avoid arbitrage possibilities, this is the same as having one market price of risk per contract (see e.g. the discussion in Benth, Cartea and Kiesel [11]). By lifting the driving noise to infinite dimensions we may achieve a richer structure on the market price of risk. On the other hand, if the market price of risk has to be a function of time-to-maturity, then the market must admit an infinite dimensional noise in its dynamics.
In the case where a L´evy fieldM =Ldrives the dynamics ofX, we introduce an infinite- dimensional version of the Esscher transform (see Benth et al. [16] and Shiryaev [45] for an account on the finite-dimensional Esscher transform in finance). We suppose that L has the characteristic triplet (µ, Q, ℓ) (see Peszat and Zabczyk [40] for a definition), and denote by ϕ(x) the cumulant function of L for x ∈U, that is, from Thm. 4.27 in Peszat and Zabczyk [40],
ϕ(x) = lnE
ei(x,L(1))
= i(µ, x)U −1
2(Qx, x)U + Z
U
ei(x,y)U −1−i(x, y)U1(|y|U≤1)
ℓ(dy).
For a θ ∈ U, suppose that ϕ(−iθ) is well-defined, and define the real-valued stochastic process
Z(t) = exp ((θ, L(t))U−ϕ(−iθ)) ,
fort∈[0, τ]. Note thatZ(0) = 1 andZis a martingale sinceLhas independent increments.
Introduce the density process
dP dQ Ft
=Z(t), t∈[0, τ],
of the probability measure P. We show thatL preserves the L´evy process property under P. By Bayes’ Formula and the independent increment property of Lunder Q, we have for t≥s,
EP
ei(x,L(t)−L(s))U| Fs
=E
ei(x,L(t)−L(s))UZ(t) Z(s)| Fs
= exp (ϕ(x−iθ)−ϕ(−iθ)).
Thus, by definition of ϕ and the symmetry of the covariance operator Q, we find that the characteristic triplet of Lunder P is (µθ, Q, ℓθ), with
(µθ, x)U = (a, x)U + Z
U
(x, y)U1(|y|U≤1)ℓθ(dy), and
ℓθ(dy) = e(θ,y)Uℓ(dy).
This defines an infinite-dimensional version of the Esscher transform, and so we introduced a market price of risk θ ∈U also for the case of L´evy noise M(t) = L(t). Observe that we do not change the covariance operator Q similar to the Girsanov transform. It is simple to extend the Esscher transform to θ being time-dependent, deterministic functions with values in U, integrable with respect to L(t). This comes, however, at the expense of losing the L´evy property of L(t) under P. On the other hand, we preserve the property of independent increments.
We end this section with a brief discussion on the generalization of our forward price dynamics to include a multiplicative structure in the noise and the drift. To this end, we introduce the stochastic partial differential equation
(2.10) dX(t) = (AX(t) +a(t, X(t)))dt+b(t, X(t))dM(t)
with initial condition X(0) =X0 ∈ H. In this case we suppose that a is a mapping from [0, τ]×H into H which has Q-a.s. Bochner-integrable trajectories in the finite interval [0, τ]. The processb : [0, τ]×H→L(H, H) is for every t∈[0, τ] and element in H a linear operator from H to H. We assume that the operators a and b satisfy the following linear growth and Lipschitz conditions for u, v ∈H and t∈[0, τ]:
ka(t, u)k2H ≤C12(1 +kuk2H) ka(t, u)−a(t, v)kH ≤C1ku−vkH , and
kb(t, u)k2LHS(H;H)≤C22(1 +kuk2H) kb(t, u)−b(t, v)kLHS(H;H)≤C2ku−vkH,
for constants C1, C2 < +∞. Then. a unique mild solution X in H to Equation (2.10) exists (see [30] and [39]), and the predictable process X : [0, τ]×Ω→H is given by (2.11) X(t) =S(t)X0+
Z t 0
S(t−s)a(s, X(s))ds+ Z t
0
S(t−s)b(s, X(s))dM(s). A specification of the process X as the solution of the stochastic partial differential equa- tion (2.10) is clearly a generalization of the dynamics defined in Equation (2.1).
We have the following regularity result, similar to Lemma 2.4
Lemma 2.8. If EkX0k2H1(D) <∞, then the unique solutionX to Equation (2.11) satisfies supt∈[0,τ]EkX(t)k2H1(D) <∞.
Again, we could assume that the forward curve (under the Musiela parameterization) is either defined by ga(t, x) = δx(X(t)) or gm(t, x) = exp(δx(X(t))) for the additive and multiplication dynamics respectively. We leave out a further analysis of these more general models here, but remark that they can be approximated and simulated with the methods discussed in Section 4 below.
3. Modeling of electricity flow forwards
In this Section we consider flow forward prices and model them using infinite dimensional dynamics. Since flow forwards deliver over a time period, we enlarge the function space to dimension 2, since now we have both time-to-maturity and length-of-delivery in the Musiela parametrization. In this case one could argue in favour of adding noise in both variables. However, a no-arbitrage condition linked to different delivery periods restricts the range of possible noise specifications in the “length-of-delivery” dimension. We provide a discussion of functions satisfying this condition.
Suppose F(t, T1, T2) is the flow forward price at time t for a contract delivering the commodity (electricity or gas, typically) over the time period [T1, T2], wheret, T1, T2 ∈[0, τ] and t ≤ T1 ≤ T2. We know that this price dynamics must be a martingale in the risk- neutral framework, i.e.
(3.1) t7→F(t, T1, T2)
is aQ-martingale. One can show (see [16]) that the dynamics must also satisfy an additional no-arbitrage condition, namely
(3.2) F(t, T1, T2) = 1
T2−T1
Z T2
T1
F(t, u, u)du .
In many electricity markets, flow forward contracts may be traded over various delivery periods, some of which may overlap. For example, in the NordPool market, one may invest in flow forwards delivering electricity in the four quarters of the year, or one may trade in a flow forward delivering over the whole year. Similarly there are contracts delivering over the different months constituting a quarter. Hence, there must be specific no-arbitrage relations between these contracts with overlapping delivery periods. In Benth et al. [16]
such conditions are presented, and it is shown that they imply Equation (3.2) when con- sidering a market where all possible delivery periods are included. Thus, models for the flow-forward priceF(t, T1, T2),which must hold for all possible choices of T1 and T2, must satisfy Equation (3.2). In an infinite dimensional framework this is exactly the case.
Let us introduce the Musiela parametrization x=T1−t and y=T2−T1, and define G(t, x, y),F(t, t+x, t+x+y).
Our dynamics for the flow forwards F needs to satisfy the no-arbitrage condition (3.2).
Translated into a Musiela parametrization, it takes the form,
(3.3) yG(t, x, y) =
Z y 0
G(t, x+z,0)dz .
Thus, any specification of G(t, x, y) has to satisfy this condition in order to prevent ar- bitrage possibilities between flow forward contracts of different length and time to matu- rity. We observe that relation (3.3) immediately suggests a way to model the flow for- ward dynamics via the use of fixed-delivery contracts studied in the previous Section. Let G(t, x,0),g(t, x), with g equal toga orgm given in Equation (2.3), resp. Equation (2.4).
DefineG(t, x, y) by
yG(t, x, y) = Z y
0
g(t, x+z)dz .
Since g(t,·) takes values in a function space of regular functions, the integral makes sense for each y > 0, if we assume, for example, H = L2(R+) and the regularity result in Lemma 2.4.
In the following we investigate functions which fulfill the no-arbitrage condition (3.2).
Notice that F(t,·,·) is symmetric in T1, T2. We remark that the no-arbitrage condition only affects the start of delivery and the end of delivery, and not the present time t. If we just look at those variables for the time being we find that the values of the function on the diagonal, i.e F(t, u, u), characterize the function on every point (T1, T2) ∈[0, τ]×[0, τ] ⊂ R2+. If we are given an arbitrary integrable function ψ(u) with an antiderivative given by Ψ(u), i.e Ψ′(u) = ψ(u) an set ψ(u) =F(t, u, u) we get from Equation (3.2)
F(t, T1, T2) = 1 T2−T1
Z T2
T1
ψ(u)du= Ψ(T2)−Ψ(T1) T2−T1
. If we set, as an example, ψ(u) = (n+ 1)un, forn ∈N, we get that
F(t, T1, T2) = T2n+1−T1n+1 T2−T1
fulfills the no-arbitrage condition. The choice ψ(u) = eu gives F(t, T1, T2) = eT2 −eT1
T2−T1
.
It is possible to derive a similar characterization for functions fulfilling the parameterized version of the no arbitrage condition, Equation (3.3). We remark as above that this condition only affects the time to maturity x and the length of delivery y and not t.
Here we choose as above an arbitrary integrable function ψ with antiderivative Ψ and set ψ(u) =G(t, u,0), then we may write
G(t, x, y) = 1 y
Z y 0
ψ(x+z)dz = 1
y(Ψ(x+y)−Ψ(x)).
As an example: if ψ(u) = eu, it follows that G(t, x, y) = 1y(ex+y −ex) is fulfilling the no- arbitrage condition. It should be mentioned that we could also choose the antiderivative Ψ as an absolutely continuous function and derive the same results.
We proceed with the Hilbert-space-valued stochastic differential equation to model di- rectly the flow forwards. Recall theH-valued stochastic processX given by the solution of the stochastic differential equation (2.1), and with mild solution given in Equation (2.2).
We suppose H is a Hilbert space of functions in (x, y), and assume that M is U-valued, with U also being a Hilbert space (again possibly a space of functions in (x, y)). We con- sider the evaluation functional δ(x,y), and extend the assumptions on the forward models from the previous Section. The discussion of some possible choices of Hilbert spaces also extends naturally. Note, however, that in order for the evaluation functional to be con- tinuous, we must increase the order of the Sobolev spaces to s > 3/2, since our functions now are defined on R2+ rather than the positive real line (recall the Sobolev Embedding Theorem (Morrey’s inequality)).
As it is shown in Benth and Koekebakker [13], it is not possible to model F(t, T1, T2) using a geometric Brownian motion whose volatility depends on T1 and T2. Only volatility specifications depending on time t are valid. This essentially rules out all relevant models since naturally the volatility should have a dependency on the delivery period due to the well-known Samuelson effect (see [13] for more a detailed discussion). It is an established empirical fact that the volatility of forwards increases as time to maturity approaches zero (see, e.g., Frestad et al. [32]), and to have reasonable models one is forced to take this into account. We conclude with this that multiplicative models for flow forwards seem difficult to specify in the context of power markets. Admittedly, the result in Benth and Koekebakker [13] was proven for real-valued noise processesM, but the arguments may be generalized to infinite dimensional noises as well.
In the remaining of this Section we focus solely on additive models, since they provide the richest structures for flexible flow forward dynamics. To this end, define
(3.4) G(t, x, y) :=δ(x,y)(X(t)).
Since we have that
F(t, T1, T2) =G(t, T1 −t, T2−T1),
the operatorA is still given by A=∂/∂x. The martingale condition on X(t) is derived as in Prop. 2.5, and becomes a(t) = 0. From now on we suppose this to hold.
We move on by studying the possible range of models G based on the specification in Equation (3.4) that will satisfy the no-arbitrage condition (3.3). First, we remark that the condition can be interpreted as a side-constraint to the stochastic differential equation, (2.1), of X, since it implicitly imposes certain structural conditions on the evolution of X. Next, since G(t, x, y) is fully specified byG(t, z,0) via the no-arbitrage condition, the
variations inlength-of-delivery y is simply an integral of those intime-to-maturity. Hence, there cannot be any separate, independent noisy behaviour. As we have discussed in the Introduction, there is ample statistical evidence for noisy behaviour of forward prices along time-to-maturity x. We are not aware of any studies of noise in y, the length-of-delivery, and how this relates to x.
In order to study the range of volatility specifications b in more detail, we restrict ourselves to driving noise processes M which fulfill the Itˆo isometry:
Ek Z t
0
S(t−s)b(s)dM(s)k2H =E Z t
0
kS(t−s)b(s)k2LHS(H,H)ds .
We remark that for the class of square integrable martingales M2b(U), introduced in Sec- tion 2, this is in general not true. However, for many processes, such as the important class of L´evy processes the isometry holds.
Proposition 3.1. Let S2 be the shift operator in the second variable onH, and suppose the initial forward curve G(0, x, y) satisfies the no-arbitrage condition (3.3). Then G(t, x, y) fulfills the no-arbitrage condition (3.3) if and only if for every 0≤s ≤t≤τ andx, y ≥0, the volatility process b has the property
(3.5)
S2(y)− 1 y
Z y 0
S(z)dz
b(s)u∈ker(δ0,0S(t−s+x)), for all u∈U, where the integration of S is understood as a Bochner integral.
Proof. By the continuity ofδx,y, we have that G(t, x, y) =δx,yS(t)X(0) +
Z t 0
δx,yS(t−s)b(s)dM(s). It follows from the no-arbitrage condition (3.2) that
G(t, x, y) = 1 y
Z y 0
G(t, x+z,0)dz
= 1 y
Z y 0
δx+z,0S(t)X(0)dz +1 y
Z y 0
Z t 0
δx+z,0S(t−s)b(s)dM(s)dz . If we compare the first terms, we may write
δx,yS(t)X(0) = 1 y
Z y 0
δx+z,0S(t)X(0)dz , or, by definition,
G(0, x, y) = 1 y
Z y 0
G(0, x+z,0)dz . By assumption, this is supposed to hold.
Next, comparing the terms with the stochastic integral, we find by linearity of the operators
Z t 0
δx,yS(t−s)b(s)dM(s) = 1 y
Z y 0
Z t 0
δx+z,0S(t−s)b(s)dM(s)dz
= Z t
0
1 y
Z y 0
δx+z,0S(t−s)b(s)dz dM(s).
Hence, if the integral process with respect to M satisfies the Itˆo isometry, we find the operator equality
δx,yS(t−s)b(s) = 1 y
Z y 0
δx+z,0S(t−s)b(s)dz . It is easily verified that, for v ∈H
δx,yS(t−s)v =δ0,yS(t−s+x)v =δ0,0S(t−s+x)S2(y)v . Applying this and the semigroup property ofS, we get,
δ0,0S(t−s+x)S2(y)b(s) = 1 y
Z y 0
δx+z,oS(t−s)b(s)dz
= 1 y
Z y 0
δ0,0S(x+z)S(t−s)b(s)dz
= 1 y
Z y 0
δ0,0S(t−s+x+z)b(s)dz
=δ0,0
1 y
Z y 0
S(t−s+x+z)b(s)dz
=δ0,0S(t−s+x) 1
y Z y
0
S(z)dz
b(s). This is equivalent to
δ0,0S(t−s+x)
S2(y)b(s)− 1 y
Z y 0
S(z)dzb(s)
= 0.
This proves the Proposition.
We note that a sufficient condition for the volatility b to ensure no-arbitrage is that S2(y)b(t) = 1
y Z y
0
S(z)dzb(t).
Recall b is an operator mappingU to H, while the shift operators S and S2 maps H into itself. Thus, for h∈U we have that, after evaluating the function b(t)h in (u, v),
(b(t)h)(u, v+y) = 1 y
Z y 0
(b(t)h)(u+z, v)dz .
We now discuss this condition for the special case when the driving noise M of X(t) is one-dimensional, that is, M(t) is a real-valued square-integrable martingale. Further we need the integral with respect to M to fulfill the Itˆo isometry.
By definition of G, we have (assuming a drift-free equation) G(t, x, y) =X(0)(x+t, y) +
Z t 0
b(s)(x+t−s, y)dM(s).
Appealing to the stochastic Fubini Theorem (see Peszat and Zabczyk [40]), implicitly assuming that b and M are sufficiently regular for this theorem to hold, we find from Equation (3.2) that the volatility b must satisfy
(3.6)
Z t 0
b(s)(t−s+x, y)dM(s) = 1 y
Z y 0
Z t 0
b(s)(t−s+x+z,0)dM(s)dz .
Let us specialize to the case of volatilities which are additive and with a multiplicative structure between x and y, that is, we suppose that
b(t)(x, y) = α(t, x)β(t, y),
for two positive functions αand β. To avoid technicalities, we assume them to be smooth.
The no-arbitrage condition imposes (3.7)
Z t 0
β(s, y)
β(s,0)dM(s) = Z t
0
1 y
Z y 0
α(s, t−s+x+z)
α(s, t−s+x) dz dM(s).
The left-hand side of Equation (3.7) is independent ofx, which means that the right-hand side must be a function independent ofx as well. This holds in particular if there exists a positive smooth function k(s, t, z) such that
α(s, t−s+x+z)
α(s, t−s+x) =k(s, t, z).
We note that k(s, t,0) = 1. In light of the Samuelson effect, it is reasonable to have a decaying volatility in time-to-maturity, meaning that x 7→ b(t)(x, y) is decreasing. This implies that z 7→ k(s, t, z) decreases as well, if we impose the Samuelson effect into the model. We find,
α(s, t−s+x+z)−α(s, t−s+x)
z = k(s, t, z)−k(s, t,0)
z α(s, t−s+x), and passing to the limitz ↓0 we write (for s =t),
∂xα(t, x) =∂zk(t, t,0)α(t, x).
Here we introduce the notation ∂x =∂/∂x. Putting c(t) :=∂zk(t, t,0), we conclude that
(3.8) α(t, x) = α(t,0)ec(t)x.
If we take the Samuelson effect in the volatility structure into account, z 7→ k(s, t, z) is decreasing. This implies that c(t) is non-positive. With α as in Equation (3.8), we find (3.9)
Z t 0
β(s, y)dM(s) = Z t
0
β(s,0)
c(s)y ec(s)y−1
dM(s).
In conclusion, a sufficient class of volatility structures b(t, x, y) satisfying the no-arbitrage condition (3.6) is given by
b(t, x, y) = h(t)
c(t)y ec(t)(x+y)−ec(t)x ,
for smooth functionshand c. If the Samuelson effect holds, then candh are non-positive.
Using thatx=T1−t and y=T2−T1, we find the volatility structure of the forward price dynamics to be
σ(t, T1, T2) :=b(t, T1, T2−T1) = h(t)
c(t)(T2−T1) ec(t)(T2−t)−ec(t)(T1−t) .
We observe that for c(t) = −κ, this volatility structure is resulting from a one-factor spot price model being an Ornstein-Uhlenbeck process with speed of mean-reversionκand “spot volatility” h(t) (see Schwartz [47] and Benth et al. [16] for details).
We remark that this result is a concrete example for the structure of the forward dy- namics imposed by the no-arbitrage condition, which we derived in the beginning of this Section. Due to the additive structure of the model here we can drive the structure only for the volatility. For the model to be conform with the no-arbitrage condition the initial condition has to have a similar structure.
As we have mentioned before, reasonable arbitrage-free geometric models seem hard to construct in the framework of flow forwards. One approach to avoid the complicating no- arbitrage condition in Equation (3.3) is to model time to maturity and length of delivery simultaneously by letting x be the time to some point within the delivery period. Thus, one introduces a dynamics g(t, x) with x =τ −t, for a given τ ∈ [T1, T2]. This approach was implicitly suggested and used in Bjerksund et al. [18] and applied for example in Audet et al.[2]. As we see, this approach brings us back to the forward price dynamics for contracts with fixed time of delivery. It is important to notice that we only obtain a flow forward price dynamics up to time τ. In some markets, like the NordPool power market, one cannot trade the forwards within their delivery period, so in fact t ≤ T1 ≤ τ is the natural modelling horizon. In other markets, like the German power exchange EEX, the flow forwards can be traded within the delivery period, however the liquidity is very low.
As we recall from Section 2, we have conditions ensuring an arbitrage-free dynamics of the geometric model which allow for a rich class of possible processes. Moreover, by choosing τ > T1 we also incorporate the empirical fact that flow forward prices do not converge to the spot price as time to delivery goes to zero.
4. Simulations and numerics
In this Section we present a numerical method for the simulation of forward curves based on numerical solution of the stochastic partial differential equation (2.1). We restrict our attention to forward contracts with fixed maturity time. Both the additive ga(t, x) in Equation (2.3) and the multiplicative gm(t, x) in Equation (2.4) are covered by our numerical approach, since we solve for X(t). However, in the examples that we present we shall be concerned with gm(t, x). The method we propose does not have any restriction on the volatility structure as in Bjørk and Gombani [19] to allow for finite dimensional representations. We equilibrate the error for the finite dimensional realization of the infinite dimensional noise with the errors introduced by the space and time discretization of the stochastic partial differential equation. Moreover, we remark that it is not restricted to
energy markets, but provides an alternative approach to simulating forward rates in fixed- income markets for example, where the HJM approach has been applied.
First, we introduce a space and time approximation for the simulation of a sample path.
To calculate functions of the solution, such as moments, we refer to the multilevel Monte Carlo method as introduced in Barth and Lang [9]. The multilevel Monte Carlo method exhibits significant improvement over the Monte Carlo method in terms of computational complexity. This is especially important if we want to solve stochastic partial differential equations, since their computation is expensive. We have to solve for each sample a partial differential equation, which entails in solving a (possibly very large) system of equations in every time-step.
To simulate the dynamics of the first-order hyperbolic, stochastic partial differential equation (2.1) we need an approximation in space, leading to the semidiscretized dynam- ics, as well as an approximation in time (full discrete form). For the finite dimensional approximation in space we introduce a Petrov–Galerkin method (SUPG method), which is well known to yield reasonable results for deterministic, first-order hyperbolic differential equations. A Galerkin approximation would provide convergence as well, but oscillations might occur when the solution is not smooth enough. Since the solution to a stochastic partial differential equation is not necessarily smooth we have to introduce an artificial diffusion via a Petrov–Galerkin approximation. There are several publications on conver- gence of these methods for deterministic equations (c.f. [38] or [24]). In general, a finite element method is preferable since it allows for more freedom when choosing the nodal points of the approximation than, for instance, a method based on finite differences.
We are facing the initial boundary value problem ut=Au+ξ ,
with initial conditionu(0,·) =v in some convex domain D⊂Rd with Lipschitz boundary
∂D. Here ut denotes the partial derivative with respect to t. For a given vector field a, the first order differential operator A is defined as
Aφ(x) :=
Xd
i=1
ai(x)φxi(x). The inflow boundary is the set
∂D−:={x∈∂D:a(x)·n(x)>0},
where n(x) denotes the exterior normal to ∂D at x. To have a well posed problem we impose a boundary condition on the inflow boundary and the initial condition u(0,·) =v matching the boundary condition in∂D−. The source term ξ is in our case the sum of the drift and the random field. As in the deterministic case, we introduce a finite-dimensional finite element space Sh as an approximation for the space variable x ∈ D. In the case studied we have D = [0, τ] ⊂ R+. Sh is a family of spaces of piecewise linear finite element functions with respect to a partition of D, without any boundary conditions on the functions in Sh, thusSh ⊂H1(D). The expectation of the solution X(t), for t∈[0, τ] is, according to Lemma 2.4, an element ofH1(D) if we have certain regularity on the initial
condition. For first-order hyperbolic equations we just have essential boundary conditions at the inflow boundary ∂D−, which leads to the following definition for the finite element space:
Sh−={χ∈ Sh : χ=c at∂D−}
for some real constant c. The parameter h is a measurement of the fineness of the ap- proximation. It denotes the longest side of any element of the triangulation in the space approximation. Translated into the problem at hand, the inflow boundary is the long end of the market. The constantc is, in this context, the value at the long end of the market, i.e. the value contracts reach when time-to-delivery tends to infinity.
The deterministic, semidiscrete, first order hyperbolic differential equation converges with rate h3/2 for the SUPG method, see, e.g., Larsson and Thom´ee [38]. As indicated in Barth [5] and Barth and Lang [6], convergence of the stochastic partial differential equation only depends on the convergence of the corresponding partial differential equation and the compatibility of the scheme with the Itˆo isometry or a martingale inequality. A convergence result for this kind of equations was derived in Barth [5].
The main difference between SUPG and Galerkin methods is the test functions used. In a Galerkin approximation, one uses the same functions as test functions and trial functions.
In contrast, to implement our Petrov–Galerkin method we choose as test functionsχ+γχ,ˆ where χ is, for instance, the hat function basis and ˆχ denotes the (weak) derivative of the basis function. The trial functions remain χ. The parameter γ has to be chosen according to the problem, leading to a numerical scheme which is problem dependent. This could be circumvented with a Discontinuous Galerkin method (DG method). However, the parameter does not depend on the noise, but rather on the continuity of the initial and/or the boundary condition.
For the fully discrete problem we have to approximate the derivative int as well. There are several options for this approximation. We use a Crank–Nicolson approach, which is known to be unconditionally stable. As a discretization of [0, τ] we introduce the equidis- tant time stepping 0 =t0 < t1 < . . . < tN =τ, whereti = Nτi, fori= 1, . . . , N,N ∈N. We couple the number of time steps with the degrees of freedom hin the space discretization, to have an equilibrated error in space and time, as suggested in Barth and Lang [9].
For the simulations we assume the driving noise process to be a Gaussian random field W or a random field of correlated L´evy processes L. For the latter, according to Peszat and Zabczyk [40] we have
L(t)(x) = X∞
k=1
(L(t), ek)Hek(x)
where L(t) is a square integrable L´evy process in H and (ek, k ∈ N) is an orthonormal basis of H. The series on the right-hand side converges in L2. For the finite-dimensional approximation of the infinite dimensional noise we truncate the sum after K terms. The rate of convergence of this truncated sum should be higher than the convergence of the full approximation scheme. In the case of a Gaussian random field we truncate the expansion (2.8). Results on how to truncate these series can be found in [7]. There the authors show
