# The heat modulated infinite dimensional Heston model and its numerical approximation

The HEat modulated Infinite DImensional Heston (HEIDIH) model and its numerical approximation are introduced and analyzed. This model falls into the general framework of infinite dimensional Heston stochastic volatility models of (F.E. Benth, I.C. Simonsen '18), introduced for the pricing of forward contracts. The HEIDIH model consists of a one-dimensional stochastic advection equation coupled with a stochastic volatility process, defined as a Cholesky-type decomposition of the tensor product of a Hilbert-space valued Ornstein-Uhlenbeck process, the mild solution to the stochastic heat equation on the real half-line. The advection and heat equations are driven by independent space-time Gaussian processes which are white in time and colored in space, with the latter covariance structure expressed by two different kernels. In the first part of the paper, regularity results for the HEIDIH model in fractional Sobolev spaces are formulated. These are achieved under smoothness conditions on the covariance kernels, which in particular allow for weighted Matérn kernels. In the second part, numerical approximation of the model is considered. An error decomposition formula, pointwise in space and time, for a semi-explicit finite-difference scheme is proven. For a special case, essentially sharp convergence rates are obtained when this is combined with a fully discrete finite element approximation of the stochastic heat equation. The analysis takes into account a localization error, a pointwise-in-space finite element discretization error and an error stemming from the noise being sampled pointwise in space. The rates obtained in the analysis are higher than what would be obtained using a standard Sobolev embedding technique. Numerical simulations illustrate the results.

• 7 publications
• 4 publications
• 4 publications
• 6 publications
12/18/2020

### Strong Rates of Convergence for Space-Time Discretization of the Backward Stochastic Heat Equation, and of a Linear-Quadratic Control Problem for the Stochastic Heat Equation

We introduce a time-implicit, finite-element based space-time discretiza...
10/31/2019

### A stochastic transport problem with Lévy noise: Fully discrete numerical approximation

Semilinear hyperbolic stochastic partial differential equations have var...
12/21/2021

### SPDE bridges with observation noise and their spatial approximation

This paper introduces SPDE bridges with observation noise and contains a...
08/23/2022

### Improved rates for a space-time FOSLS of parabolic PDEs

We consider the first-order system space-time formulation of the heat eq...
03/19/2020

### Parameter estimation for discretely sampled stochastic heat equation driven by space-only noise revised

The main goal of this paper is to build consistent and asymptotically no...
09/10/2019

### Weak convergence of fully discrete finite element approximations of semilinear hyperbolic SPDE with additive noise

We consider the numerical approximation of the mild solution to a semili...
10/15/2020

### Approximation of BSDE with Hidden Forward Equation and Unknown Volatility

In the present paper the problem of approximating the solution of BSDE i...

## 1. Introduction

The use of SPDEs (stochastic partial differential equations) for the modeling of the risk-neutral dynamics of forward prices in commodity markets is well-established, see, e.g, [2, 7, 8, 10]. By a forward price we refer to the price at of some contract on a commodity (oil, gas, an instantaneous delivery of electricity) delivering at a maturity time . The basic idea comes from the Heath, Jarrow and Morton approach [23]

. A heuristic is the following: assume that the forward price is modeled by an Itô process via

. Here the initial value as well as the standard Brownian motion and the stochastic volatility depend on the maturity time . Let denote the time to maturity, consider a space-time Gaussian process and let , . Then we define a space-time process by

 (1) X(t,x):=f(0,t+x)+∫t0σ(s,x+t)dB(s,x)=S(t)X(0)(x)+∫t0(S(t−s)g(s))(x)dB(s,x)

for . Here is the left-shift semigroup defined by for functions , having the infinitesimal generator . This is made rigorous by letting be the mild solution at time to a stochastic advection equation on driven by a Wiener process in a Hilbert space following the setting of [16]. Since the term structure of forward prices is believed to be continuous and eventually constant, is typically chosen to be a Hilbert space of functions displaying these two features. This is a simple model that nevertheless captures non-Gaussian effects (by letting the volatility term be stochastic) as well as correlation structures in time to maturity.

In this paper, we present, analyze and study the numerical approximation of the HEIDIH model, which is a special case of the SPDE approach to forward price modeling. Compared to other such models, our work stands out for three reasons.

1. [label=()]

2. In the context of forward pricing, the properties of are often discussed in terms of its incremental covariance kernel (or function) , see, e.g., [2, 10]. It is, however, rare to find any discussion on what properties a given kernel should have in order for to be well-defined and thus allow for the existence of in a given Hilbert space . In this work, we give concrete conditions on the covariance kernels involved that allow us to make precise statements not only on the existence of but on its regularity in terms of fractional Sobolev spaces.

3. We derive convergence rates for numerical approximations of . While numerical simulations of SPDEs for forward pricing abound, rigorous derivation of sharp convergence rates for the algorithms considered are rare. We mention [5, 6] but note that these only consider one-dimensional noise (i.e., for all ). Such convergence results are needed to guarantee the reliable simulation of forward prices of contracts and options written on them. Furthermore, in contrast to many other numerical analyses of SPDEs, our work focuses on pointwise-in-space errors for . This is practically relevant since evaluated at the point corresponds to the forward price of a contract with maturity time .

4. The model we construct is parsimonious but sufficiently flexible to allow practitioners to match features of forward prices encountered in real-world data. We return to this point below.

The HEat modulated Infinite DImensional Heston (HEIDIH) model is a special case of the infinite dimensional stochastic volatility model introduced in [4, 10]. Therein, a price process and a volatility process are coupled via another process in the system

 (2)

on a fixed interval . Here is a Hilbert space of functions on on which the left-shift semigroup is strongly continuous. With being another strongly continuous semigroup, we have written and for the generators of and , respectively. By and we denote two independent Wiener processes in . The mapping is not function-valued as in (1) but operator-valued and is defined by . The stochastic process is -valued and takes values on the unit ball in and denote the Hilbert tensor product. The process is an intrinsic part of the model and each choice yields different dynamics. One example is the setting for all , where is a deterministic function with unit norm. Another is obtained by letting for and otherwise. In any case , where is the adjoint of , so that

is a square root of the stochastic variance process

[10, Proposition 7]. This is the reason the authors of [10] refer to this system as a Heston stochastic volatility model in Hilbert space.

The HEIDIH model is obtained by specifying in (2), which is the Laplacian with either Dirichlet or Neumann boundary conditions at coupled with a diffusivity constant . For the state space we adopt the setting of [18] and choose where is a fractional Sobolev space and . An example of a realization of and for this model when for all and fixed is shown in Figures 1(a) and 1(b), see Remark 4.6 for precise parameter choices.

The generator has been chosen to illustrate how precise regularity results and approximation convergence rates can be obtained for Heston stochastic volatility models in Hilbert spaces under suitable conditions on the covariance kernels of and . We believe that our results on kernels, regularity in Sobolev spaces and pointwise-in-space approximation results are of interest also for other SPDE models of forward prices. Moreover, the model itself is flexible enough to capture real-world features of forward prices. Specifically:

1. [label=()]

2. The covariance kernels of and of can be chosen to match the regularity (i.e., for which the process takes values in

) of the prices being modeled. First steps towards a theory for direct estimation of the covariances structures of infinite-dimensional stochastic integrals are taken in

[9].

3. A practitioner can influence the behavior of the volatility process by specifying the diffusivity constant . The process can be seen to be Gaussian in due to the noise being additive, and is therefore completely determined by its covariance function . By an argument that is made rigorous below, the left-shift property of the semigroup implies that

 X(t,x)=X(0,x+t)+∫t0Y(s,x+s)⟨Z(s),⋅⟩HdW(s).

By properties of the heat kernel, as (see Figure 1(c)) and the rate of decrease can be increased by letting . The degree to which the dependence of the volatility with regards to time to maturity is local can therefore be controlled. The local dependence is a realistic feature of the HEIDIH model since contracts close in maturity affect each other.

4. The model directly displays the Samuelson effect. This says that the volatility should be decreasing in time to maturity (see, e.g., [3] for a discussion of this in a commodity pricing context) i.e., that as the time to maturity . This is obtained from the fact that , see also Figure 1(a).

We now outline our results in more detail. In Section 2, we describe some necessary mathematical preliminaries, starting with notation and operator theory. After this, we reiterate properties of reproducing kernel Hilbert spaces (RKHSs) and fractional Sobolev spaces that we need, and define our state space .

In Section 3 we introduce cylindrical Wiener processes with covariance kernels and prove some preliminary results of a general nature. First, we completely characterize the kernels that yield classical -Wiener processes in . Then, we derive a class of kernels (which include weighted Matérn kernels) that allow for precise regularity and numerical convergence results in the remainder of the paper.

Section 4 concerns existence, uniqueness and regularity of and . First, we show that under suitable conditions on and , classical results guarantee existence, uniqueness and Sobolev regularity of . We prove a temporal regularity result for the process with fixed that is needed for numerical approximation, and note that it is better than what would be obtained from a Sobolev embedding technique. These results are then used when the HEIDIH model is introduced. Under further assumptions on , and , we show that is well-defined in for suitable .

In Section 5 we turn our attention to numerical approximation of the HEIDIH model. In Section 5.1 we derive an error decomposition formula for a finite-difference approximation of . This applies to any Heston stochastic volatility model in a Hilbert space as considered in [10], as long the Hilbert space in question has a reproducing kernel. In Section 5.2 we focus on the approximation of in the special case of Dirichlet boundary conditions and . The analysis is complicated by the fact that the errors are considered pointwise in space. This is rarely considered for finite element approximations, which we use for the approximation of . We also need to take into account localization errors, meaning that we approximate the process on a fixed domain instead of . Since we wish to employ rapid noise simulation methods such as the circulant embedding method, we consider pointwise sampling of rather than the typical case of sampling , where denotes the -projection onto a finite element grid. This leads to an additional error. We take all of these error sources into account and show in numerical simulations that our results are sharp. Finally, in Section 5.3 we put these results together to derive convergence rates for the full HEIDIH model and discuss how our regularity results influence practitioners’ choice of algorithms, specifically what space-time grids to use for the components and . Section 6 finishes the paper with a discussion of how our results relate to existing results for numerical approximations of SPDEs and an outline of future work.

## 2. Preliminaries

This section briefly presents the mathematical machinery necessary for our results.

### 2.1. Notation and operator theory

Let and be Banach spaces. All Banach spaces in this paper are taken over unless otherwise stated. We denote by the space of bounded linear operators from to equipped with the usual operator norm. If and are separable Hilbert spaces we write for the space Hilbert–Schmidt operators. This is a separable Hilbert space with an inner product, for an arbitrary ONB (orthonormal basis) of , given by

 ⟨Γ1,Γ2⟩L2(U,V):=∞∑j=1⟨Γ1ej,Γ2ej⟩V for Γ1,Γ2∈L2(U,V).

We have if and only if and their norms coincide. We use the shorthand notations and when and we write for the class of positive semidefinite operators when is a Hilbert space. For operators , we say that they are trace class if for one, equivalently all, ONBs of .

For Hilbert spaces and , the tensor is regarded as an element of by the relation for and .

We need the concepts of approximation and entropy numbers from [33]. Given a linear operator between Banach spaces and , its th approximation number is given by

 ψn(Γ):=inf{∥Γ−~Γ∥:~Γ∈L(U,V),rank(~Γ)

, while its th entropy number is

 ϵn(Γ):=inf{δ≥0:Γ(BU)⊂∪mj=1{vj+δBV} for some v1,…,vm∈V and m≤2n−1}.

Here and denote the unit spheres in and , respectively. By [33, Theorems 11.2.3, 12.1.3], both types of numbers satisfy

 (3) ϵ1(Γ)=ψ1(Γ)=∥Γ∥L(U,V).

Both numbers are also multiplicative in that, if and for some additional Banach space , then

 (4) ϵm+n−1(Γ1Γ2)≤ϵm(Γ1)ϵn(Γ2) and ψm+n−1(Γ1Γ2)≤ψm(Γ1)ψn(Γ2)

for all [33, Theorems 11.9.2 and 12.1.5]. If and are Hilbert spaces, then

 (5) ψn(Γ)≤2ϵn(Γ)

for all and

 (6) ∥Γ∥L2(U,V)=∥(ψn(Γ))∞n=1∥ℓ2,

where denotes the usual space of -summable sequences [33, Theorems 11.3.4, 12.3.1 and 15.5.5]. Taking (3) and (4) into account, this means in particular that is an operator ideal.

For two Banach spaces and , we write if and for some constant and all , i.e., the embedding operator . For Hilbert spaces, we write as shorthand for

Throughout this paper, we adopt the notion of generic constants, which may vary from occurrence to occurrence and are independent of any parameter of interest, such as spatial or temporal step sizes. By , for , we denote the existence of a generic constant such that .

### 2.2. Reproducing kernel Hilbert spaces

We make heavy use of the theory of reproducing kernel Hilbert spaces (RKHSs). The properties that we are going to need are listed here, we refer to, e.g., [11, 40] for further details. Throughout the paper we consider symmetric positive semidefinite kernels on a non-empty subset , , but usually just refer to as a kernel on when there is no risk of confusion. A Hilbert space is said to be the RKHS of a kernel if it is a Hilbert space of real-valued functions on a non-empty index set , , such that the conditions

1. [label=()]

2. for all and

3. for all ,

are satisfied. The property 2 is referred to as the reproducing kernel property of the space . For each kernel there is one and only one Hilbert space of functions on with these two properties. Moreover, a Hilbert space of functions on is a RKHS if and only if the evaluation functional , defined by , is continuous for all [40, Theorem 10.2]. We write for when it is clear from the context what index set we have in mind. Since is positive semidefinite, for all .

If is separable with an ONB , [11, Theorem 14] yields the kernel decomposition , the sum being convergent in . This implies that with convergence in . For index sets that are domains in , continuity of is a sufficient condition for separability [11, Theorem 15].

Suppose that we are given two RKHSs on an index set with kernels and . Then if and only if there is some constant such that is a positive semidefinite kernel [1, Theorem I.13.IV, Corollary I.13.IV.2]. We express this by writing .

For any Hilbert spaces and operator , we may interpret the range as a Hilbert space when equipped with the inner product , where is the pseudoinverse of [34, Proposition C.0.3]. The norm on may also be represented by

 ∥u∥Γ(U)=minv∈UΓv=u∥v∥U,

see [34, Remark C.0.2]. Applying this to a RKHS and the operator that restricts functions on to functions on a subset , it follows from [11, Theorem 6] that maps to and with equal norms. Moreover, by [40, Theorem 10.46] there is a linear extension operator of functions on to functions on such that .

### 2.3. Fractional Sobolev spaces

For a domain , , we write , , for the usual Banach space of -integrable equivalence classes of functions with respect to the Lebesgue measure. We omit when it is clear from the context which domain we refer to. The fractional Sobolev space , , is a separable Hilbert space and a generalization of the usual Sobolev space , . It consists of all such that

 ∥u∥2Hr(\amsmathbbRd):=(2π)−d/2∫\amsmathbbRd|^u(ξ)|2(1+|ξ|2)rdξ<∞,

where

is the Fourier transform of

. For , the fractional Sobolev space , , consists of all for which for some . It is equipped with the norm

 ∥u∥Hr(E):=minv∈Hr(\amsmathbbRd)v|E=u∥v∥Hr(\amsmathbbRd).

When is sufficiently regular (such as when it is a possibly unbounded interval in the case of ), has a bounded extension operator. It therefore coincides with the usual Sobolev space with equivalent norms. We refer to [36, Chapter VI] for more details on this. In terms of the restriction operator , we have with equal norms. We again omit from when it is clear from the context which domain we refer to.

If , then , , is a RKHS with a stationary (that is, for ) kernel given by

 mr(x):=21−r(Γ(r))−1|x|r−d/2Kr−d/2(|x|),

see [40, Chapter 10]. Here denotes the modified Bessel function of the second kind and order . Note that is continuous and bounded on . By the reproducing property, we obtain that all functions are continuous and bounded, which is a special case of the classical Sobolev embedding theorem.

We also make use of the weighted spaces where . For these we introduce the weight function given by for . We set . These are special cases of the weighted spaces considered in [22], which we refer to for further details, see also [38, Theorems 2.3.9 and 2.5.6]. They are Hilbert spaces equipped with the inner product . For , with the embedding being compact if and only if both inequalities are strict and [22, Theorem 2.3]. For , we note that for all and that for ,

 wα(x)f(x)=⟨wαf,mr(x−⋅)⟩Hr=⟨f,w−αmr(x−⋅)⟩Hr,α,x∈\amsmathbbRd.

From this, we obtain that is a RKHS with kernel .

The state space that we will work with for the HEIDIH model is given by , . We define it for . It consists of all where and . It is equipped with the norm

 ∥f∥2Hr:=∥~f∥2Hr(\amsmathbbR+)+b2.

When , it is a RKHS with a kernel given by . We consider as a subspace of by identifying it with .

## 3. Cylindrical Wiener processes in fractional Sobolev spaces

In this section, we recall the concept of cylindrical Wiener processes with a focus on processes that have a covariance kernel. We prove a result that completely characterizes -Wiener processes in the state space . Then, we outline the assumptions on the kernels and that are made in the forthcoming sections. Finally, we construct a class of kernels that fulfill all the assumptions made.

Consider

, a filtered probability space

fulfilling the usual conditions and a Hilbert space . We follow [35] and define a (strongly) cylindrical Wiener process as a process such that is a real-valued Wiener process for all . This definition agrees with what is called a generalized Wiener process in [16].

Given a separable RKHS with ONB and a sequence of real-valued Wiener processes, we construct a cylindrical Wiener process in by

 (7) W(t)u:=∞∑j=1βj(t)⟨u,ej⟩H,t∈[0,T],u∈H,

which converges in , see [35, Remark 7.3]. Note that for . The representation does not depend on in the sense that if is another ONB of and we define a real-valued Wiener processes by , then

 W(t)u=∞∑j=1~βj(t)⟨u,fj⟩H,t∈[0,T],u∈H,\amsmathbbP-a.s.

If is another Hilbert space such that , we obtain a cylindrical Wiener process in by replacing in (7) with . In , will have an incremental covariance operator given by in the sense that for all . Similarly, has covariance in . These operators have unique positive square roots. By [16, Corollary B.3], . In this sense, the distribution of does not depend on the choice of in (7). Therefore, we are justified in calling a cylindrical Wiener process with kernel and need not specify in which Hilbert space we consider it.

The SPDEs we consider in the paper are built on stochastic integrals with respect to cylindrical Wiener processes . Consider a predictable process . Itô integrals taking values in are well-defined with

 \amsmathbbE[∥∥∥∫t0Ψ(s)dW(s)∥∥∥2H]=∫t0\amsmathbbE[∥Ψ(s)∥2L2(Hq,H)]ds,

for , provided that the integral on the right hand side is finite. We refer to [35] and [16, Section 4.2] for more details on Wiener processes in Hilbert spaces and the Itô integral.

If , the sum

 W(t):=∞∑j=1βj(t)ej,t∈[0,T],

converges in . It is then called a -Wiener process, and its covariance operator is trace-class. It induces a cylindrical Wiener process by and we do not make a notational difference between the two concepts. Many papers dealing with SPDE models for forward prices assume that is a -Wiener process, see, e.g., [4, 10]. It is therefore important to clarify when this is the case in our setting with the state space . The following theorem completely characterizes the kernels that satisfy , so that is a -Wiener process in .

###### Theorem 3.1.

There is a separable RKHS such that if and only if there is a separable RKHS such that , a function and a constant with

 q(x,y)=⟨~q(x,⋅)+f,~q(y,⋅)+f⟩H~q+c=~q(x,y)+f(x)+f(y)+∥f∥2~Hq+c

for all .

If this holds for , then the cylindrical Wiener process with kernel is an -valued -Wiener process and fulfills

 \amsmathbbE[⟨W(1),u⟩H0⟨W(1),v⟩H0] =∫\amsmathbbR+×\amsmathbbR+~q(x,y)~u(x)~v(y)dxdy+∫\amsmathbbR+au~v(x)f(x)+av~u(x)f(x)dx+auav(∥f∥2H~q+c)

for all .

If this holds for , then is a mean square continuous and bounded random field for all and for all .

Note that if for some , without the property of the embedding, there might not be a -Wiener process with covariance function . However, as long as , this theorem shows that we can still interpret as the covariance function of in a weaker sense. One might also consider embeddings in negatively weighted spaces. We do not pursue this direction but focus instead on constructing kernels such that .

###### Proof of Theorem 3.1.

First, suppose that there is a separable RKHS such that . Let be an ONB of and write with . Since

 (8) ∞∑j=1∥fj∥2Hr+c2j=∥IHq↪Hr∥L2(Hq,Hr)<∞

we have . Moreover, due to the fact that , the evaluation operator is well-defined on . Since

 ∞∑j=1|fj(x)|2=∞∑j=1|ej(x)−cj|2≤2∞∑j=1c2j+2∞∑j=1|ej(x)|2=2∥(cj)∞j=1∥2ℓ2+2q(x,x)<∞,

we may define a kernel by

 ~q(x,y)=∞∑j=1fj(x)fj(y)

for . The kernel is symmetric and positive semidefinite. We define a Hilbert space by

 H:={v=∞∑j=1vjfj,(vj)∞j=1⊂\amsmathbbR such that ∥v∥2:=∞∑j=1v2j<∞}.

An ONB of this space is . Moreover, is a RKHS: for and ,

 |v(x)|2≤(∞∑j=1|vj||fj(x)|)2≤∞∑j=1|vj|2∞∑j=1|fj(x)|2=∥v∥2H~q(x,x).

Since also and , we find that is the RKHS of . The fact that follows directly from (8). Let . Then

 q(x,y) =∞∑j=1(fj(x)+cj)(fj(y)+cj) =∑j∈Jfj(x)fj(y)+∑j∈Jcjfj(x)+∑j∈Jcjfj(y)+∑j∈Jc2j+∑j∈Jcc2j,

where the split is justified by and for all . By setting and we obtain one direction of the first claim of Theorem 3.1. The other direction is obtained by analogous arguments.

For the second claim, consider the same setting as above. Since ,

 (9) \amsmathbbE[⟨W(1),u⟩Hr⟨W(1),v⟩Hr]=⟨Qu,v⟩Hr=⟨I∗Hq↪Hru,I∗Hq↪Hrv⟩Hq=∞∑j=1⟨I∗Hq↪Hru,ej⟩Hq⟨I∗Hq↪Hrv,ej⟩Hq=∞∑j=1⟨u,ej⟩Hr⟨v,ej⟩Hr.

Applying this identity with yields

 \amsmathbbE[⟨W(1),u⟩H0⟨W(1),v⟩H0] =∞∑j=1(⟨~u,fj⟩L2+aucj)(⟨~v,fj⟩L2+avcj) =∑j∈J⟨~u,fj⟩L2⟨~v,fj⟩L2+∑j∈Javcj⟨~u,fj⟩L2+∑j∈Jaucj⟨~v,fj⟩L2+∞∑j=1c2j =∫\amsmathbbR+×\amsmathbbR+~q(x,y)~u(x)~v(y)dxdy+∫\amsmathbbR