Dirichlet Forms and Finite Element Methods for the SABR Model
DDIRICHLET FORMS AND FINITE ELEMENT METHODS FOR THESABR MODEL
BLANKA HORVATH AND OLEG REICHMANN
Abstract.
We propose a deterministic numerical method for pricing vanilla options under theSABR stochastic volatility model, based on a finite element discretization of the Kolmogorovpricing equations via non-symmetric Dirichlet forms. Our pricing method is valid under mildassumptions on parameter configurations of the process both in moderate interest rate environ-ments and in near-zero interest rate regimes such as the currently prevalent ones. The parabolicKolmogorov pricing equations for the SABR model are degenerate at the origin, yielding non-standard partial differential equations, for which conventional pricing methods —designed fornon-degenerate parabolic equations— potentially break down. We derive here the appropriateanalytic setup to handle the degeneracy of the model at the origin. That is, we construct anevolution triple of suitably chosen Sobolev spaces with singular weights, consisting of the do-main of the SABR-Dirichlet form, its dual space, and the pivotal Hilbert space. In particular,we show well-posedness of the variational formulation of the SABR-pricing equations for vanillaand barrier options on this triple. Furthermore, we present a finite element discretization schemebased on a (weighted) multiresolution wavelet approximation in space and a θ -scheme in timeand provide an error analysis for this discretization. Introduction
The stochastic alpha beta rho (SABR) model introduced by Hagan et. al in [41, 43] is todayindustry standard in interest rate markets. The model with parameters ν > β ∈ [0 , ρ ∈ [ − , X t = Y t X βt d W t , X = x > , d Y t = νY t d Z t , Y = y > , d (cid:104) Z, W (cid:105) t = ρ d t, ≤ t ≤ T < ∞ , where W and Z are ρ -correlated Brownian motions on a filtered probability space (Ω , F , ( F t ) t ≥ , P ).The SABR process ( X, Y ) takes values in the state space D = [0 , ∞ ) × (0 , ∞ ), and describes thedynamics of a forward rate X with stochastic volatility Y and with the initial values x > y >
0. The constant elasticity of variance (CEV) parameter β determines the general shape ofthe volatility smile and the parameter ν (often denoted by α ) governs the volatility the stochasticvolatility. The first pricing formula for the SABR model (the so-called Hagan formula) proposedin [41, 43], is based on an expansion of the Black-Scholes implied volatility for an asset driven by(1.1). This tractable and easy-to-implement asymptotic expansion of the implied volatility madecalibration to market data easier. This, and the model’s ability to capture the shape and dynamics(when the current value X = x of the asset changes) of the volatility smile observed in the marketare virtues of the SABR model, which soon became a benchmark in interest rates derivatives mar-kets [5, 8, 10, 65]. The Hagan expansion is only accurate when the expansion parameter is smallrelative to the strike, that is when time to maturity or the volatility of volatility ν is sufficientlysmall. For low strike options such as in low interest rate and high volatility environments, much like Date : January 10, 2018.2010
Mathematics Subject Classification.
Key words and phrases.
SABR model, Finite Element Methods, Dirichlet Forms.BH is grateful for financial support from the Swiss National Science Foundation (SNSF Grant 165248). Theauthors are also indebted to the anonymous referees for their helpful comments and to Robbin Tops for severalinvaluable discussions. The views expressed in this article are those of the authors and do not necessarily representthose of the European Investment Bank. a r X i v : . [ q -f i n . M F ] J a n BLANKA HORVATH AND OLEG REICHMANN the ones we are facing today, this formula can yield a negative density function for the process X in (1.1), which leads to arbitrage opportunities. Therefore, as the problem of negative densitiesand arbitrage became more prevalent, it has been addressed for example in [6, 7, 10, 29, 42] bydifferent approaches, some suggesting modifications of the SABR model or its implied volatilityexpansion. The attempt of suggesting suitable modifications to the original model is an intricatechallenge since the Hagan expansion is deeply embedded in the market and fits market pricesclosely in moderate interest rate environments. Any model that deviates from its prices in thoseregimes may be deemed uncompetitive. This makes such pricing techniques desirable, which areapplicable to the original model in all market environments. There exist several refinements tothis asymptotic formula: in [62] a correction the leading order term is proposed, and [63] providesa second-order term. In the uncorrelated case ρ = 0 the exact density has been derived for theabsolutely continuous part of the distribution of X on (0 , ∞ ) in [6, 33, 48] and the correlated casewas approximated by a mimicking model. However, it seems that these refinements have not fullyresolved the arbitrage issue near the origin. Recent results [29, 38, 39] focus on the singular partof the distribution and suggest an explanation for the irregularities appearing at interest ratesnear zero; and [38, 39] provides a means to regularize Hagan’s asymptotic formula at low strikesfor specific parameter configurations, based on tail asymptotics derived in [26, 37].We propose here a numerical pricing method for the (original) SABR model (1.1) with rathermild assumptions on the parameters. It is consistently applicable in all market environments andallows for the derivation of convergence rates for the numerical approximations of option prices.The most popular numerical approximation methods which were considered so far for the SABRmodel (or closely related models) fall into the following classes: probabilistic methods— comprisingof path simulation of the process combined with suitable (quasi-) Monte Carlo approximation—were considered in [18] for the SABR model and [1, 2, 17, 21, 47] for related models. Furthermore in[18] some difficulties of Euler methods in the context of SABR are discussed. Splitting methods—where the infinitesimal generator of the process is decomposed into suitable operators for whichthe pricing equations can be computed more efficiently—provide a powerful tool in terms of com-putation efficiency for sufficiently regular processes. Such methods are considered in [12] for amodel closely related to SABR (see also [11]), and in [28] for a large class of models. However, theapplicability of corresponding convergence results to the SABR model itself is not fully resolved.Among fully deterministic PDE methods are most notably finite difference methods, which wereconsidered in [5, 53] for the modification [42] of the SABR model (1.1) and finite element meth-ods, which were described in the context of mathematical finance in [75]. In the recent textbook[44] finite element methods have been applied to a large class of financial models—including theclosely related process (1.3)—and provide a robust and flexible framework to handle the stochasticfinesses of these models. In spite of this, finite element approximation methods did not appearin the context of the SABR model so far in the corresponding literature. For a broad review ofsimulation schemes used for the SABR model (1.1) in specific parameter regimes, see also [57] andthe references therein.Standard theory provides convergence of the above methods if the considered model satisfiescertain (method-specific) regularity conditions. However in the case of the SABR model, obtainingconvergence rates is non-standard for these methods: The degeneracy of the SABR Kolmogorovequation at the origin violates the assumptions needed in conventional finite difference methodsand—for a range of parameters—also those of ad-hoc (i.e. unweighted) finite element methods.Path simulation of the SABR process also requires non-standard techniques due to the degeneracyof the diffusion (1.1) at zero. Nonstandard techniques often become necessary for the numericalsimulation of a stochastic differential equation, when the drift and diffusion ( b and σ ) do notsatisfy the global Lipschitz condition(1.2) | b ( x ) − b ( y ) | + | σ ( x ) − σ ( y ) | ≤ C | x − y | , for x, y ∈ R n and a constant C > x and y , cf. [68]. The degeneracy of theSABR model (1.1) at X = 0, originates from failure of condition (1.2) for the CEV process (cid:101) X , IRICHLET FORMS AND FINITE ELEMENT METHODS FOR THE SABR MODEL 3 described for parameters α >
0, and β ∈ [0 ,
1] by the equation(1.3) d (cid:101) X t = α (cid:101) X βt d W t (cid:101) X = (cid:101) x > , ≤ t ≤ T < ∞ . Although the exact distribution of the CEV process is available [55], simulation of the full SABRmodel based on it can in many cases become involved and expensive. In fact, exact formulasdecomposing the SABR-distribution into a CEV part and a volatility part are only available inrestricted parameter regimes, see [6, 32, 48] for the absolutely continuous part and [38, 39] for thesingular part of the distribution.A simple space transformation (see (1.5) below) makes some numerical approximation resultsfor the CIR model (the perhaps most well-understood degenerate diffusion) applicable to certainparameter regimes of the SABR process. The CIR process(1.4) d S t = ( δ − γS t ) d t + a √ S t d W t , S = s > , ≤ t ≤ T < ∞ , with a > , δ ≥ γ = 0, a = 2 to a squared Bessel process withdimension δ on the positive real line, the connection to CEV is then made via(1.5) ϕ : R ≥ −→ R ≥ s (cid:55)−→ − δ/ s − δ , where β = − δ − δ , for δ (cid:54) = 2 , that is, assuming absorbing boundary conditions at zero, the law of S in (1.4) (for the parameters γ = 0, a = 2) under the space transformation ϕ in (1.5) coincides with the law of (cid:101) X in (1.3).Recent results in this direction, exploring probabilistic approximation methods for diffusions wherethe global Lipschitz continuity (1.2) is violated, can be found for example in [17, 21] and [47], seealso [3, 27] and the references therein. Establishing strong convergence rates in the case 2 δ < a where the boundary is accessible as in [47], is of particular difficulty, as this case renders coefficientsof the SDE (1.4) neither globally, nor locally Lipschitz continuous on the state space. Yet, theseconvergence results do not cover the parameter range of the SABR model. Further approximationschemes are presented in [1, 2] which apply to CIR processes with accessible boundary and bothstrong and weak convergence for the approximation are studied. The weak error analysis of Talayand Tubaro [74] yields second order convergence of the schemes proposed in [1, 2], which coversthe parameters 0 < β < but the results do not directly carry over to the case ≤ β < well posedness of the latter. We use the weighted multiresolution(wavelet) Galerkin discretization of [13] in the state space to approximate variational solutions ofthe SABR-Kolmogorov pricing equations for financial contracts (vanilla and barrier options). Forthe time discretization of the semigroup generated by the process (1.1) we propose a θ -scheme.We derive approximation estimates tailored to our weighted setup, measuring the error betweenthe true solution of the pricing equations and their projection to the discretization spaces. Basedon these, we conclude error estimates akin to [64] for our fully discrete scheme. Under appropriateregularity assumptions on the payoff we obtain the full convergence rate for our finite element ap-proximation. The advantage of the presented method is that it allows for a consistent pricing withvery mild parameter assumptions on the SABR process and it is robustly applicable for moderateinterest rate environments as well as in the current low interest rate regimes. Furthermore, theproposed discretization can be applied to compute prices of compound options or multi-periodcontracts without substantial modifications of the numerical methodology, cf. [68].The article is organized as follows. Section 2 is devoted to the formulation of the SABR pricingproblem in the appropriate analytic setting and an outline of the general idea of the variational BLANKA HORVATH AND OLEG REICHMANN analysis underlying the proposed finite element method for the SABR model. In Section 2 weintroduce some notations. In Section 2.1 we introduce the SABR Dirichlet form and cast thevariational formulation of the SABR pricing equation in a suitable setting. We then propose aGelfand triplet of spaces for our finite element discretization, consisting of a space V of admissiblefunctions (the domain of the SABR-Dirichlet form), its dual space V ∗ and a pivotal Hilbert space H , containing V . In Section 2.2 we briefly recall some relevant existing results to prove well-posedness of the SABR-pricing problem on the triplet V ⊂ H ⊂ V ∗ , and conclude the existence ofa unique weak solution to the variational formulation of the Kolmogorov partial differential equa-tions on these spaces. We furthermore derive in this section a non-symmetric Dirichlet form forthe SABR model, thereby extending the results of [22] on Dirichlet forms on SABR-type modelsto the non-symmetric case. In Section 3 we present the finite element discretization of the weaksolution of the equation examined in the previous sections. Section 3.1 is devoted to the spacediscretization, which is carried out through a spline wavelet discretization of spaces V and H . Wereview the multiresolution spline wavelet analysis of [44, 64] (in the unweighted case) to discretizethe volatility dimension. The forward dimension (the CEV part) is more delicate, due to its de-generacy at zero. In this case we apply the weighted multiresolution norm equivalences, provenin [13] which are suitable to this degeneracy. We pass from the univariate case to the bivariatecase by constructing tensor products of the discretized spaces in each dimension as outlined in[44]. Finally, we specify the mass- and stiffness matrices involved in the space discretization. InSection 3.2 we present the fully discrete scheme by applying a θ -scheme in the time-stepping. Wefollow [64], to conclude that the stability of the θ -scheme continues to hold in the present setting ofweighted Sobolev spaces. In Section 4 we derive error estimates for our finite element discretiza-tion. In Section 4.1 approximation estimates of the projection to our discretization spaces areestablished based on multiresolution (weighted) norm equivalences. We cast our estimates underspecification of different regularity assumptions on the solution of our pricing equations. We usethese estimates in Section 4.2 to derive convergence rates for our finite element discretization andconclude that under some regularity assumptions on the payoff, examined in the previous section,we obtain the full convergence rate. We remark here, that the approximation- and error estimatespresented in Section 4 for the SABR model readily yield the corresponding approximation- anderror estimates for the CEV model as a direct corollary. The well-posedness of the variationalformulation of the CEV pricing equations has been studied in [44, 68], however to the best of ourknowledge, a presentation of the full error analysis thereof was not available in the correspondingliterature so far. 2. Preliminaries and problem formulation
Preliminaries and Notations: there exists a unique weak solution of the system (1.1) whichwill be established via the associated martingale problem [50, Theorem 21.7], and by (pathwise)uniqueness of (1.1), via [50, Theorem 23.3]. Furthermore, the process X in (1.1) is a martingalewhenever β < β = 1 it is a martingale if and only if ρ ≤ V the notation ||·|| V ≈ ||·|| V indicates that the norms are equivalent on V . Function spaces of bivariate functionswill be denoted in italic ( V , H , . . . ) and spaces of univariate functions by ( V, H, . . . ) accordingly.For a domain G ⊂ R (resp. interval I ⊂ R ) we denote by L loc ( G ) (resp. L loc ( I )) the locallyintegrable functions on G (resp. I ), and by C ∞ ( G ) (resp. C ∞ ( I )) the smooth functions withcompact support. Derivatives with respect to time will be denoted by ˙ u , ¨ u , . . . accordingly to easenotation.2.1. Analytic setting for the SABR model.
In this section we establish the triplet
V ⊂ H ⊂V ∗ of spaces, tailored to the SABR process on which we cast the Kolmogorov pricing equations invariational form. We then proceed to show the well-posedness of these pricing equations on thistriplet, i.e. that inequalities (2.17) and (2.18) are fulfilled. We conclude the section by provingthat the bilinear form resulting from the weak formulation of the pricing equations is indeed anon-symmetric Dirichlet form corresponding to the (unique) law of the SABR process. Fix a time IRICHLET FORMS AND FINITE ELEMENT METHODS FOR THE SABR MODEL 5 horizon [0 , T ] and set (cid:101) Y = log( Y ), so that the SDE (1.1) becomes(2.1) dX t = X βt exp( (cid:101) Y t ) dW t X = x > ,d ˜ Y t = νdZ t − ν dt ˜ Y = log y , y > d (cid:104) W, Z (cid:105) t = ρdt. where we impose absorbing boundary conditions at X = 0 to ensure martingality of the process.The solution to (2.3) (and (1.1)) is then uniquely defined. Indeed, for the parameters β ∈ [0 . , β ∈ [0 , . t ∈ J = (0 , T ) of a European-type contract on (2.1) with payoff u is (2.2) u ( t, z ) = E [ u ( Z zτ )] , t ∈ J where τ := ( T − t ) and Z zτ := ( X τ , (cid:101) Y τ ) is the process started at z := ( x, ˜ y ) ∈ R ≥ × R , with( x, ˜ y ) = ( X t , ˜ Y t ) P -a.s. Then for u ∈ C , ( J ; R ≥ × R ) ∩ C ( ¯ J ; R ≥ × R ) the Kolmogorov pricingequation to (2.1) is(2.3) ˙ u ( t, z ) − Au ( t, z ) = g ( t, z ) in J × R ≥ × R ,u (0 , z ) = u ( z ) in R ≥ × R where g ∈ L ( J, V ) ∩ H ( J, H ) denotes a general forcing term (see [44, Section 4]) with V ∗ as in(2.12) and where the infinitesimal generator A of (2.1) reads(2.4) Af = x β e y ∂ xx f + ρνx β e ˜ y ∂ x ˜ y f + ν ∂ ˜ y ˜ y f − ν ∂ ˜ y f for f ∈ C ( D ) ⊂ D ( A ) . From now on we drop the tilde in the logarithmic volatility for notational convenience. Theoperator A in (2.4) is a linear second order operator, degenerate (i.e. non-elliptic) at the boundary { ( x, y ) : x = 0 , y > } . The domain D ( A ) of the operator A is equipped with a norm || · || V (specified in Definition 2.4 below) and the completion of D ( A ) under this norm will be denoted by V . Furthermore, we denote by H (specified in Definition 2.2 below) a separable Hilbert space—henceforth the pivot space —containing V , such that V (cid:44) → H is a dense embedding. The innerproduct ( · , · ) H of H is extended to a duality pairing ( · , · ) V ∗ ×V on V ∗ × V , where V ∗ denotes thedual space of V , equipped with the dual norm || · || V ∗ . Identifying the Hilbert space H with itsdual H ∗ we obtain the corresponding Gelfand-triplet(2.5) V (cid:44) → H ∼ = H ∗ (cid:44) → V ∗ , where (cid:44) → denotes a dense embedding.With view to the discretization, it is customary (cf. [44, 64, 68]) to localize the spatial domainof the PDE ( R ≥ × R in (2.3)) at this point to a bounded domain G ⊂ R ≥ × R with Lipschitzboundary. In what follows, all localization domains are rectanguar G = [0 , R x ) × ( − R y , R y )denoting the range of admissible values which can be taken by the price (and volatility) process. Remark 2.1.
For call and put options the error made by truncating the domain to G ⊂ R × R ≥ corresponds to approximating the option prices by a knock-out barrier options, up to the firsthitting time of the boundary ∂G , see [44, Sections 5.3 and 6]. The estimate in [44, Theorem5.3.1] can be directly applied to the volatility dimension Y , and yields that the truncated problemconverges to the original problem exponentially fast. Furthermore, for the CEV process in the assetprice the estimate serves as an upper bound by comparison between the CEV ( β ∈ [0 , β = 1. The probabilistic argument to estimate the localizationerror by a knock-out barrier option was suggested by Cont and Voltchkova in [19, 20] even in amore general setting of L´evy models. Indeed, for the SABR model, the probability that the first For notational simplicity we consider zero risk-free interest rates here. Furthermore, we make some standardtechnical assumptions on the payoff function u : In accordance with [44, Equation (5.10) page 47] we assume u to satisfy u (0) = 0 and a polynomial growth condition (satisfied by vanilla contracts), and that u ∈ H (seeDefinition 2.2) in accordance with [64, Equations (2.10) and (2.12)]. BLANKA HORVATH AND OLEG REICHMANN hitting time of R x resp. R y occurs before T converges to zero as R x , R y → ∞ . The lower boundaryhowever cannot be truncated to any domain ( (cid:15), R x ) for a positive (cid:15) without possibly introducing asignificant localization error. This is due to the fact that the SABR process accumulates a positivemass at zero for every T > β <
1. For details see [38, 39], where the mass at zero iscalculated for several relevant parameter configurations.
Definition 2.2.
Let G := [0 , R x ) × ( − R y , R y ) ⊂ R + × R , R x , R y > G we define the weighted space(2.6) H := L ( G, x µ/ ) = { u : G → R measurable | || u || L ( G,x µ/ ) < ∞} with µ ∈ (cid:40) [ − β,
0] for β ∈ [0 , )[ − , − β ] for β ∈ [ , , (2.7)where || u || L ( G,x µ/ ) := ( u, u ) H for the bilinear form(2.8) ( u, v ) H := (cid:90) G u ( x, y ) v ( x, y ) x µ dxdy, u, v ∈ V. Remark 2.3.
Note that ( H , ( · , · ) H ) is a Hilbert space, see Appendix A.1, Lemma A.5 and A.3.A possible choice for the weight is µ = − β for any β ∈ [0 , β ∈ [ ,
1) and β ∈ [0 , ) and choose µ = − β for β ∈ [ , µ = 0 for β ∈ [0 , ) ∪ { } .Distinguishing the above parameter regimes has the advantage that we preserve the classical settingof an unweighted L ( G )-space for the parameters β ∈ [0 , ) ∪ { } . The latter choice furthermorehighlights that our analytic setting consistently extends the univariate CEV case to the bivariateSABR case: see [44, Section 4.5] and Remark 2.8 for the analytic setting for CEV. Definition 2.4.
Set G := [0 , R x ) × ( − R y , R y ) ⊂ R + × R , R x , R y >
0, the domain of interest. Forthe first coordinate we consider L ([0 , R x ) , x µ/ ) = { u ∈ L ([0 , R x )) with || x µ/ u || L < ∞} . On L ([0 , R x ) , x µ/ ) we define the space V x := C ∞ ([0 , R x )) ||·|| Vx where || u || V x := || x β + µ/ ∂ x u ) || L (0 ,R x ) + || x µ/ u || L (0 ,R x ) , u ∈ C ∞ ([0 , R x )) . For the second coordinate we consider on L ( − R y , R y ) the space V y := H ( − R y , R y ) , where || u || V y := || ∂ y u || L ( − R y ,R y ) + || u || L ( − R y ,R y ) , u ∈ H ( − R y , R y ) . We define for the bivariate case V := (cid:0) V x ⊗ L ( − R y , R y ) (cid:1) (cid:92) (cid:16) L ([0 , R x ) , x µ/ ) ⊗ V y (cid:17) . (2.9)The dual space will be denoted by V ∗ and equipped with the usual dual norm(2.10) || v || V ∗ = sup u ∈V ( v, u ) V ∗ ×V || u || V , v ∈ V ∗ . Remark 2.5.
Note that the norm on the space (2.9) is by construction equivalent to(2.11) || u || V ≈ || x β + µ/ ∂ x u || L ( G ) + || x µ/ ∂ y u || L ( G ) + || x µ/ u || L ( G ) , for u ∈ V . Lemma 2.6.
For any µ in (2.7) the spaces H , V and V ∗ in Definitions 2.2 and 2.4 form a Gelfandtriplet. In particular, the inclusion maps are dense embeddings (2.12) V (cid:44) → H (cid:44) → V ∗ , where V and H are specified in (2.6) and (2.9) . The scalar product ( · , · ) H can hence (by Lemma 2.6) be extended to a dual pairing ( · , · ) V×V ∗ . See Appendix A.2.1 for details. The analogous statement in the one-dimensional ( V x ) part was presented in the CEV-analysis in [44, Equation(5.33), page 56]. IRICHLET FORMS AND FINITE ELEMENT METHODS FOR THE SABR MODEL 7
Proof.
The first inclusion follows by construction (cf Definition 2.4) and the second by directapproximation in the weighted space H = L ( G, x µ/ ) by smooth functions with compact support. (cid:3) We extend the inner product ( · , · ) H of the Hilbert space to the dual pairing ( · , · ) V ∗ ×V as describedabove . Recalling that A ∈ L ( V ; V ∗ ) we apply the duality pairing ( · , · ) V ∗ ×V : A bilinear form a ( · , · ) : V × V → R is thus associated with the operator A in (2.4) by setting(2.13) a ( u, v ) := − ( Au, v ) V ∗ ×V , u, v ∈ V , where the operator A acts on V in the weak sense, see in (2.14) below. Hence, we define theSABR-bilinear form by the relation (2.13), which therefore reads a ( u, v ) = (cid:90) (cid:90) G x β + µ e y ∂ x u ∂ x v dxdy + β + µ (cid:90) (cid:90) G x β + µ − e y ∂ x u v dxdy + ρν (cid:90) (cid:90) G x β + µ e y ∂ x u ∂ y v dxdy + ρν (cid:90) (cid:90) G x β + µ e y ∂ x u v dxdy + ν (cid:90) (cid:90) G x µ ∂ y u ∂ y v dxdy + ν (cid:90) (cid:90) G x µ ∂ y u v dxdy, u, v ∈ V , (2.14)which is obtained from (2.13)—by the divergence theorem together with V ⊂ L loc ( G )—when Au , u ∈ V are interpreted as weak derivatives . With the spaces V ⊂ H ⊂ V ∗ and the bilinear form(2.13) at hand, we can formulate the variational (or weak) framework corresponding to (2.3). Themotivation for passing to the weak formulation is that for degenerate equations (such as (2.3)) itis often not possible to find a classical solution u ∈ C , ( J, R ) ∩ C ( ¯ J, R ) to the original problem(2.3). The variational reformulation (2.15) problem may then still permit a (weak) solution u withless regularity. Whenever a (weak) solution of the variational problem is sufficiently smooth, itcoincides with the solution of the corresponding original problem. Definition 2.7 (Variational formulation of the SABR pricing equation) . Let V , V ∗ and H (resp. L ( G, x µ/ )) be as in Definitions 2.4 and 2.2, and let the bilinear form a ( · , · ) on V be as in (2.14).Furthermore, let u ∈ H (resp. u ∈ L ( G, x µ/ )) and consider for a T > J = (0 , T ). Then the variational formulation of the SABR pricing problem reads as follows: Find u ∈ L ( J ; V ) ∩ H ( J ; V ∗ ), such that u (0) = u , and for every v ∈ V , ϕ ∈ C ∞ ( J )(2.15) − (cid:90) J ( u ( t ) , v ) H ˙ ϕ ( t ) dt + (cid:90) J a ( u ( t ) , v ) ϕ ( t ) dt = (cid:90) J ( g ( t ) , v ) V ∗ ×V ϕ ( t ) dt. The time derivative of a function u in the appropriate Bochner space is understood in the weaksense: For u ∈ L ( J ; V ), its weak derivative in ˙ u ∈ L ( J, V ∗ ) ∩ H ( J ; V ∗ ) is defined by the relation(2.16) (cid:90) J ( ˙ u ( t ) , v ) V ∗ ×V ϕ ( t ) dt = − (cid:90) J ( u ( t ) , v ) V ∗ ×V ˙ ϕ ( t ) dt, see [44, Sections 2.1 and 3.1] for definitions and properties of Bochner spaces. Remark 2.8 (The CEV case: ν = 0) . In [44, 68] a corresponding analytic setup is studied forthe univariate case: For the CEV model the Gelfand triplet V ⊂ H ∼ = H ∗ ⊂ V in [44, 68] consistsof the weighed spaces H := L ((0 , R ) , x µ/ )and V := C ∞ ([0 , R )) ||·|| V with || u || V := || x β + µ/ ∂ x u ) || L (0 ,R ) + || x µ/ u || L (0 ,R ) such as its dual V ∗ , where the parameter µ ∈ [max {− , − β } , − β ] is chosen as in Definitions2.2, and 2.4. Indeed, setting ν = 0, the SABR process (1.1) (resp. (2.1)) with trivial volatilityprocess reduces to the CEV model, and the state space G reduces to [0 , R x ) ⊂ R . Accordingly, for ν = 0 the spaces H , V , V ∗ in Definitions 2.2 and 2.4 coincide with the spaces V , H and V ∗ above.Also the SABR bilinear form (2.14) reduces to the corresponding CEV bilinear form. Hence, Note that for this the operator A need not be self-adjoint nor needs the associated bilinear form be symmetric. Multiplying − Au with v ∈ C ∞ and integrating gives (2.13), partial integration (cf. (2.16)) then yields (2.14). BLANKA HORVATH AND OLEG REICHMANN our analytic setup extends the univariate setup of the CEV model consistently to the bivariateSABR-case. See [68, Equation (21)],[44, Equations (4.30) and (4.33)] such as [44, page 62] for thedefinitions of V, H, V ∗ and [44, Equation (4.28)] for the CEV bilinear form.2.2. Well-posedness of the variational pricing equations and SABR Dirichlet forms.
In this section we show that the triplet of spaces
V ⊂ H ∼ = H ∗ ⊂ V ∗ (Definitions 2.2 and 2.4) istailored to the degeneracy of the infinitesimal generator (2.4) at zero: in the sense that in thissetting (as a consequence of Theorem 2.11 and well-posedness cf. Theorem 2.13), the variationalformulation of the SABR pricing equation has a unique solution in V and hence the family ofoption prices P t = E [ u ( Z t )], t ≥ H , cf. Theorems 2.12 and 2.13.Key result in this section is the well-posedness of the SABR variational equations, established inTheorem 2.13. Furthermore, we show that the SABR-bilinear form (2.14) is a (non-symmetric)Dirichlet form with domain V on the Hilbert space H , cf. Theorem 2.19. The latter extends theresults of [22] on SABR-Dirichlet forms. We start by briefly recalling some concepts and resultsused in this section. Definition 2.9 (Continuity) . The form a ( · , · ) is called continuous on V if there exists a 0 < C < ∞ such that(2.17) ∀ u, v ∈ V : | a ( u, v ) | ≤ C || u || V || v || V . Definition 2.10 (Coercivity, G˚arding inequality) . The form (2.13) is said to satisfy the G˚ardinginequality on V , if there exists a constant C ≥ ∀ u ∈ V : a ( u, u ) ≥ C || u || V − C || u || H . If (2.18) holds with C = 0, the form a ( · , · ) is coercive and the equivalence a ( · , · ) ≈ || · || V holds. Theorem 2.11.
Let V and H be separable Hilbert spaces with a continuous dense embedding V (cid:44) → H . Furthermore, let a : V × V → R be a bilinear form satisfying the inequalities (2.17) and (2.18) . Then the corresponding variational parabolic problem (recall Definition 2.7) has a uniquesolution in L ( J ; V ) ∩ H ( J ; V ∗ ) .Proof. See [54, Theorem 4.1] for a proof. (cid:3)
Theorem 2.12.
Consider a bilinear form a ( · , · ) : V × V → R associated with an A ∈ L ( V , V ∗ ) via (2.13) . If a ( · , · ) satisfies the properties (2.17) and (2.18) , then − A is the infinitesimal generator ofa bounded analytic C -semigroup ( P t ) t ≥ in V ∗ . In this case, for given u ∈ H and g ∈ L ( J ; V ∗ ) ,the unique variational solution of the corresponding equation can be represented as (2.19) u ( t ) = P t u + (cid:90) t P t − s g ( s ) ds. Proof.
See [59, Theorem 2.3.], [64, Section 2. Equation (2.13) and Remark 2.1] and also [54]. (cid:3)
We now formulate the main theorem in this section:
Theorem 2.13 (Well-posedness of the SABR pricing equation) . For every configuration ( β, | ρ | , ν ) ∈ [0 , × [0 , × R + of the SABR parameters, which satisfy the condition | ρ | ν < and for any x , y > , the variational formulation (2.15) of the pricing equation (2.3) corresponding tothe SABR model (2.3) admits a unique solution u ∈ L ( J, V ) ∩ H ( J, V ∗ ) for any forcing term Note that
V, H, V ∗ are presented here in a form which is adjusted to our current notation. It is standard (see for example[64, Remark 2.4]) that the (weaker) G˚arding inequality can be reduced to thecoercivity property by the substitution v := e − C t u . In case (2.18) is fulfilled for the operator A at u , then (2.18)with C = 0 is fulfilled at v . Then the operator A + C I is coercive and solves the related problem˙ v ( t, z ) + ( A + C I ) v ( t, z ) = e − C t g ( t, z ) in t ∈ J, z ∈ R . By Theorem 2.11 above.
IRICHLET FORMS AND FINITE ELEMENT METHODS FOR THE SABR MODEL 9 g ∈ L ( J, V ∗ ) and any u in H . The unique variational solution of the pricing equation can berepresented as (2.20) u ( t, z ) = P t u ( z ) + (cid:90) t P t − s g ( s ) ds, t ≥ , z ∈ G for a strongly continuous semigroup ( P t ) t ≥ on H with the infinitesimal generator A in (2.4) .Proof. This is a direct consequence of Theorem 2.11 applied to Lemmas 2.14 and 2.16 below,which establish continuity (2.17) and the G˚arding inequality (2.18) for the SABR Dirichlet form(2.14) on this triplet. (cid:3)
Lemma 2.14.
The bilinear form (2.14) is continuous: There exists a constant C > such that | a ( u, v ) | ≤ C || u || V || v || V , for all u, v ∈ V . (2.21) Proof.
The continuity statement (2.21) is a direct consequence of the following six estimates, eachof which corresponds to a component of a ( · , · ) in (2.14):(1) (cid:82) (cid:82) G x β + µ e y ( ∂ x u )( ∂ x v ) dxdy ≤ (cid:16) || x β + µ/ e y ∂ x u || L ( G ) + || x β + µ/ e y ∂ x v || L ( G ) (cid:17) , (2) By the Cauchy-Schwarz inequality, β + µ (cid:82) (cid:82) G x β + µ − e y ( ∂ x u ) vdxdy ≤ β + µ (cid:0)(cid:82) (cid:82) G x β + µ − e y v dxdy (cid:1) / (cid:0)(cid:82) (cid:82) G x β + µ e y ( ∂ x u ) dxdy (cid:1) / , and an upper bound for the latter is derived from Hardy’s inequality [44, (5.56) p. 54] : ≤ β + µ | β + µ − | || x β + µ/ e y ∂ x v || L ( G ) || x β + µ/ e y ∂ x u || L ( G ) ≤ β + µ | β + µ − | (cid:16) || x β + µ/ e y ∂ x v || L ( G ) + || x β + µ/ e y ∂ x u || L ( G ) (cid:17) (3) ρν (cid:82) (cid:82) G x β + µ e y ( ∂ x u )( ∂ y v ) dxdy ≤ | ρν | (cid:16) || x β + µ/ e y ∂ x u || L ( G ) + ν || x µ/ ∂ y v || L ( G ) (cid:17) ,(4) ρν (cid:82) (cid:82) G x β + µ e y ( ∂ x u ) vdxdy ≤ | ρν | (cid:16) || x β + µ/ e y ∂ x u || L ( G ) + ν || x µ/ v || L ( G ) (cid:17) , (5) ν (cid:82) (cid:82) G x µ ( ∂ y u )( ∂ y v ) dxdy ≤ ν (cid:16) || x µ/ ∂ y u || L ( G ) + || x µ/ ∂ y v || L ( G ) (cid:17) ,(6) ν (cid:82) (cid:82) G x µ ( ∂ y u ) vdxdy ≤ ≤ ν (cid:16) || x µ/ ∂ y u || L ( G ) + || x µ/ v || L ( G ) (cid:17) . (cid:3) Remark 2.15.
The proof of Lemma 2.14 reveals that analogous estimates are valid if the domain G is the whole (not-truncated) state space R × R ≥ . In the next lemma, an analogous statementholds true if in estimate (2) and (6) integration by parts is valid with vanishing boundary terms. Lemma 2.16.
The bilinear form (2.14) satisfies the G˚arding inequality, i.e. there exist constants C > and C ≥ such that a ( u, u ) ≥ C || u || V − C || u || H , for all u ∈ V . (2.22) Proof.
The G˚arding inequality (2.22) is obtained from the following estimates:(1) (cid:82) (cid:82) G x β + µ e y ∂ x u∂ x udxdy = || x β + µ/ e y ∂ x u || L ( G )9 We use Hardy’s inequality as in [44, (5.56) p. 54], setting ε = 2 β + µ − C = β + µ : C (cid:82) (cid:82) x ε − v ( x, y ) dx e y dy ≤ C (cid:82) | ε − | (cid:82) x ε | ∂ x v ( x, y ) | dxe y dy . (2) Using ( ∂ x u ) u = ∂ x ( u ) and integration by parts in the second term of (2.14) yields β + µ (cid:82) (cid:82) G x β + µ − e y ∂ x ( u ) dxdy = − β + µ (2 β + µ − (cid:82) (cid:82) G x β + µ − e y u dxdy The last term is non-negative if and only if µ ∈ [ − β, − β ], which is satisfied by (2.7).(3) ρν (cid:82) (cid:82) G x β + µ e y ( ∂ x u )( ∂ y u ) dxdy ≥ −| ρ | ν (cid:16) δ || x β + µ/ e y ∂ x u || L ( G ) + δ || x µ/ ∂ y u || L ( G ) (cid:17) , for a constant δ > ρν (cid:82) (cid:82) G x β + µ e y ( ∂ x u )( ∂ y u ) dxdy ≥ −| ρ | ν (cid:16) (cid:15) || x β + µ/ e y ∂ x u || L ( G ) + (cid:15) || x µ/ u || L ( G ) (cid:17) , for a constant (cid:15) > ν (cid:82) (cid:82) G x µ ∂ y u∂ y udxdy = ν || x µ/ ∂ y u || L ( G ) ,(6) ν (cid:82) (cid:82) G x µ ( ∂ y u ) udxdy = ν (cid:82) (cid:0)(cid:82) u ( ∂ y u ) dy (cid:1) x µ dx = ν (cid:82) (cid:0) − (cid:82) u ( ∂ y u ) dy (cid:1) x µ dx = 0 by inte-gration by parts, the (Dirichlet) boundary conditions for u ∈ C ∞ ( G ) and by the densityof C ∞ ( G ) in V .Hence, a ( u, u ) ≥ (cid:16) − | ρ | ν δ − | ρ | ν (cid:15) (cid:17) || x β + µ/ e y ∂ x u || L + (cid:16) ν − | ρ | ν δ (cid:17) || x µ/ ∂ y u || L − | ρ | ν (cid:15) || x µ/ u || L , which yields the inequality (2.22) a ( u, u ) ≥ C (cid:16) || x β + µ/ e y ∂ x u || L ( G ) + || x µ/ ∂ y u || L ( G ) + || x µ/ u || L ( G ) (cid:17) − C || x µ/ u || L ( G ) = C || u || V − C || u || H , with C := min { ν − | ρ | ν δ , − | ρ | ν δ − | ρ | ν (cid:15) } and C := C + | ρ | ν (cid:15) . (cid:3) It remains to verify that the constants δ and (cid:15) can be chosen accordingly such that C > Lemma 2.17.
For every configuration ( β, | ρ | , ν ) ∈ [0 , × [0 , × R + of the SABR parameters,which satisfy the condition | ρ | ν < there exist constants (cid:15) > and δ > such that C > . Remark 2.18 (Discussion of the parameter restrictions) . Note that the case ρ = 0 , ν > C > ν = 0 yields the CEV model, as discussed in Remark 2.8. Furthermore,the condition on the parameters in Lemma 2.17 is satisfied for any | ρ | ∈ (0 ,
1] for example if0 < ν < √
2. The latter condition on ν is fulfilled in most practical scenarios observed in themarket as usual values of this parameter are well below √
2: The volatility of volatility typicallycalibrates to values around ν = 0 .
2; 0 .
4; 0 .
6, see for example [41, 62, 66].
Proof of Lemma 2.17.
For any parameter configuration with | ρ | ν < δ in such a way that(2.23) | ρ | ν > δ > , and the constants (cid:15) accordingly, such that(2.24) ν | ρ | − δ > (cid:15) > . If the inequalities (2.23) and (2.24) are satisfied then C > .
24) poses no contradiction to (2 . δ are(2.25) δ ∈ (cid:16) | ρ | ν , | ρ | ν (cid:17) . Indeed if | ρ | ν <
2, then the set in (2.25) is nonempty, and there exists an (cid:15) satisfying (2.24). (cid:3)
Theorem 2.19 (The non-symmetric SABR Dirichlet form) . Let the bilinear form a ( · , · ) be as in (2.14) and its domain V as in (2.9) . Then the pair ( a ( · , · ) , V ) is a (non-symmetric) Dirichlet formon the Hilbert space ( H , ( · , · ) H ) in (2.6) , for every parameter configuration ( β, | ρ | , ν ) ∈ [0 , × [0 , × R + with | ρ | ν < , and for any µ as in (2.7) . IRICHLET FORMS AND FINITE ELEMENT METHODS FOR THE SABR MODEL 11
Proof.
The crucial statement in the above theorem is that the pair ( a, V ) is a coercive closed form on the Hilbert spaces H and L ( G, x µ/ ) for any µ as in (2.7), where a denotes SABR-bilinearform (2.14) and V is the space (2.9). For this we first define an auxiliary symmetric bilinear formrelated to (2.14) as:(2.26) E ( u, v ) := (cid:82) (cid:82) G x β + µ e y ∂ x ( u ) ∂ x ( v ) dx dy + (cid:82) (cid:82) x µ ∂ y ( u ) ∂ y ( v ) dx dy D ( E ) := C ∞ ( G ) . It follows directly from [71, Corollary 3.5] or [14, Proposition 1] that for any open G ⊂ R ≥ × R ,the bilinear form (2.26) is closable on the Hilbert space H = L ( G, x µ/ ). Then, closability of( a, C ∞ ( G )) in H and L ( G, x µ/ ) is inherited from the closability of the auxiliary symmetric form E in (2.26) via the equivalence of norms a ( · , · ) ≈ ||·|| V , see [70, Section 3, Proposition 3.5] . The latterequivalence is a direct consequence of the continuity property (2.21) and G˚arding inequality (2.22),which were proven for the form a ( · , · ) on the triplets V ⊂ H ⊂ V ∗ such as V ⊂ L ( G, x µ/ ) ⊂ V ∗ inSection 2.2. Hence, in the inequality (2.21) the norm || · || V on the right hand side can be replacedby ( a ( · , · )) / , yielding the strong (resp. weak) sector condition (see (B.1) and Remark B.3) forthe SABR-bilinear form: There exists a C >
0, such that for all u, v ∈ V| a ( u, v ) | ≤ C ( a ( u, u )) / ( a ( v, v )) / resp. | a ( u, v ) | ≤ C ( a ( u, u )) / ( a ( v, v ))) / for a ( u, v ) := a ( u, v ) + ( u, v ) H , u, v ∈ V . cf. [70, Equations (2.4) resp. (2.3)]. Hence, ( a ( · , · ) , V ) isa coercive closed form on ( H , ( · , · ) H ) and on ( L ( G, x µ/ ) , ( · , · ) H )), in the sense of [70, Definition2.4], see Section B.2 below.Now for the coercive closed form ( a ( · , · ) , V ) to be a Dirichlet form (cf. Definition B.4) it remainsto show that it is sub-Markovian (i.e. it satisfies contraction properties (2.27)). These contractionproperties follow via Theorem B.7 from the respective contraction properties of the unique (cf.Theorem B.5) semigroups on ( H , ( · , · ) H ) and ( L ( G, x µ/ ) , ( · , · ) H )) associated to ( a ( · , · ) , V ).The contraction properties (sub-Markovianity) in Definition B.2 for the bilinear form (2.14) canalso be shown directly: For any u ∈ V it holds that u + ∧ ∈ V (since derivatives are taken in theweak sense) and the the contraction properties are by non-negativity of the form equivalent to(2.27) a ( u + u + ∧ , u − u + ∧ ≥ a ( u + ∧ , u − u + ∧ ≥ ,a ( u − u + ∧ , u + u + ∧ ≥ a ( u − u + ∧ , u + ∧ ≥ . Since the functions u + ∧ u − u + ∧ a ( u + ∧ , u − u + ∧
1) =0 follows directly from the construction (2.14) of the bilinear form, yielding sub-Markovianity.Hence, the form a ( u + ∧ , u − u + ∧
1) is a non-symmetric (as a ( u, v ) (cid:54) = a ( v, u ) in general) Dirichletform. (cid:3) For the sake of completeness, we included in Appendix B the properties of non-symmetric Dirichletforms which are involved in Theorem 2.19. For a comprehensive treatment of symmetric and non-symmetric Dirichlet forms, see the monographs [15, 34] and [70] respectively.
Remark 2.20.
To extend the above proof to the untruncated problem, see Remark 2.15 above.We remark moreover that the form (2.14) is the (unique) Dirichlet form corresponding to the(SABR) Markov process (2.3). Recall, that the law of the process (2.14) is unique, by pathwiseuniqueness of the solutions of (1.1) when imposing Dirichlet boundary conditions at zero, and thecorresponding martingale problem (see [50, Theorem 21.7], and [50, Lemma 21.17]) is well-posed.3.
Discretization
In this section we derive a suitable discretization in space and time for the variational formula-tion of the SABR pricing equations. We propose a multiresolution approximation inspired by the(unweighted) wavelet discretization in [64, Section 3.4]. To accommodate to the current setting,we shall rely on the weighted multiresolution analysis established in [13, Sections 5.2 and 5.3], forfurther reference see also [68].
Space discretization and the semidiscrete problem.
Given u ∈ V and g ∈ L ( J ; V ∗ ),we first choose an approximation u (0 ,L ) ∈ V L of the initial data u , where V L ⊂ V is a finitedimensional subspace. Then the semi-discrete problem reads as follows: Find u L ∈ H ( J ; V L ),such that u L (0) = u (0 ,L ) (3.1) ( ddt u L , v L ) H + a ( u L , v L ) = ( g ( t ) , v L ) V ∗ ×V , ∀ v L ∈ V L . (3.2)We assume u (0 ,L ) = P L u , were P L : V → V L , u (cid:55)−→ u L is an appropriate projector (see equation(3.17) below). The semi-discrete problem is an initial value problem for N = dim V L ordinarydifferential equations(3.3) M ddt u ( t ) + A u ( t ) = g ( t ) , u (0) = u , where u ( t ) denotes the coefficient vector of u L ( t ) for t ≥
0, and M and A denote the mass andstiffness matrix respectively with respect to some basis of the discretization space V L , which weconstruct in the subsequent sections.3.1.1. Discretization spaces.
Recall that the bivariate pricing equation (2.3) is parabolic withdegenerate elliptic operator A . To accommodate to the degeneracy of A in (2.4), we introducedthe weighted spaces H and V with weights which are singular at the boundary x = 0 of thedomain G (cf. Definitions 2.2 and 2.4) to establish well-posedness of the variational formulationof the pricing equations. To construct finite element approximation spaces for V and H we firstestablish univariate approximation in each dimension separately. The approximation spaces tothe (weighted) univariate spaces V x , V y and H x , H y are then assembled to obtain the bivariateapproximation spaces to V and H . From now on we restrict ourselves without loss of generality tothe unit interval I = [0 ,
1] if not stated otherwise. Our multiresolution analysis on I ⊂ R consistsof a nested family of spaces(3.4) V ⊂ V ⊂ . . . ⊂ V L ⊂ . . . ⊂ L ( I, ω ) =: H ( I, ω ) , where the choice of the spaces V i (this choice is specified in Section 3.1.2 below) is such, that theinclusions (3.4) are valid, and such that (cid:83) l ∈ N V l = L ( I, ω ) holds, where L ( I, ω ) denotes a spaceof square integrable functions on I with weight function ω . The latter inclusion and convergencestatements for our choice of spaces V i are justified in Theorem 3.2 in the following section.To specify the choice of V l in (3.4), we first specify our choice of partitions of the unit interval I at discretization level L ∈ N : For this, let L ∈ N denote the discretization level, and let T be agiven initial partition of I . We assume that for any L >
0, the family {T , . . . , T L } of partitions ofthe unit interval I is such, that for each 0 < l < L the partition T l is obtained from the partition T l − by bisection of each of its subintervals. Hence, for any l ∈ { , . . . , L } , the partition T l has N l = C l subintervals, where C denotes the number of intervals of the initial partition T . Notethat with this discretization, for any l ∈ { , . . . , L } the dimension of the space V l (specified inSection 3.1.2) is dim V l = C l =: N l , ≤ l < L, dim V L = C L =: N (3.5)for a constant C >
0. Furthermore, the codimension on each level l is M l := N l +1 − N l , ≤ l < L. (3.6)3.1.2. Wavelets for L -spaces on an interval. As announced above, we now specify the choice of V l in (3.4). For the univariate approximation spaces of the H x := L ( I, ω x ) , H y := L ( I ) suchas V x := H ( I, ω x ) , V y := H ( I ) on the interval I , we recapitulate basic concepts and defini-tions of (bi-)orthogonal, compactly supported wavelets from [13, Section 2], and [69, Section 6.2].We consider two-parameter wavelet systems { ψ l,k } , l = 0 , . . . , ∞ , k ∈ { , . . . M l } of compactlysupported functions ψ l,k . Here the first index, l , denotes the “level” of refinement resp. reso-lution: wavelet functions ψ l,k with large values of the level index are well-localized in the sensethat diam supp( ψ l,k ) = O (2 − l ). The second index, k ∈ M l , measures the localization of wavelet ψ l,k within the interval I at scale l and ranges in the index set M l . We refer to [68, Section 4] IRICHLET FORMS AND FINITE ELEMENT METHODS FOR THE SABR MODEL 13 for a graphic illustration of the construction in a specific example. In order to achieve maximalflexibility in the construction of wavelet systems (which can be used to satisfy other requirements,such as minimizing their support size or to minimize the size of constants in norm equivalences,see (3.11),(3.12) cf. [69, Section 6.2]), we propose wavelet bases, which are biorthogonal in L ( I ).These consist of a primal wavelet system { ψ l,k } , l = 0 , ..., ∞ , k ∈ M l which is a Riesz basis of L ( I ) (and which enter explicitly in the space discretizations) and a corresponding dual waveletsystem { (cid:101) ψ l,k } , l = 0 , ..., ∞ , k ∈ M l , which are not used explicitly in the algorithms, see [69,Section 6.2]. The primal wavelet bases ψ l,k span the finite dimensional spaces V l = span { ψ i,j | ≤ i ≤ l, ≤ j ≤ M l } , l > , that is V l = V l − (cid:77) W l , (3.7)that is, V l = V l − (cid:76) W l inductively, where W l = span { ψ l, , . . . , ψ l,M l } . Hence,(3.8) V L = (cid:77) ≤ l
0, which only depends on the weight function. (w2) The boundary wavelets ψ k,l ( x ) and dual wavelets (cid:101) ψ k,l ( x ) are denoted by the index set ∇ Ll = { k ∈ N , γ − ≤ k ≤ l − , ∈ supp ψ k,l } , (cid:102) ∇ lL = { k ∈ N , γ − ≤ k ≤ l − , ∈ supp (cid:101) ψ k,l } , and the boundary wavelets and their duals satisfy the conditions | ψ k,l ( x ) | ≤ C ψ l/ (2 l x ) γ , | ( ψ k,l ) (cid:48) ( x ) | ≤ C ψ l/ (2 l x ) γ − , γ ∈ N , k ∈ ∇ Ll | ˜ ψ k,l ( x ) | ≤ C ψ l/ (2 l x ) ˜ γ , | ( ˜ ψ k,l ) (cid:48) ( x ) | ≤ C ψ l/ (2 l x ) ˜ γ − , ˜ γ ∈ N , k ∈ (cid:101) ∇ Ll , where γ, (cid:101) γ ∈ N are parameters such that α + γ > − and − α + ˜ γ > − is fulfilled forthe parameter α in ( w v ∈ V has a representation as a series and any v L ∈ V L as a linear combination(3.9) v L = L (cid:88) l =0 M l (cid:88) j =1 v lj ψ l,j , v L ∈ V L ; v = ∞ (cid:88) l =0 M l (cid:88) j =1 v lj ψ l,j , v ∈ V, where v lj = ( v, (cid:101) ψ l,j ) L ( I ) , cf. [13, Equation (2.9)]. Approximations of functions in V are obtainedby truncating the wavelet expansion, cf. [64, Equation (3.13)] or [44, p. 161]. Definition 3.1 (Projection Operator) . For a subspace V of a (possibly weighted) Hilbert space L ( I, ω ) over an interval I (cf. (3.4)) the projection to the univariate finite element discretizationspace P L : V → V L is defined by truncating the wavelet expansion at refinement level L (3.10) P L v := L (cid:88) l =0 M l (cid:88) j =1 v lj ψ l,j , v ∈ V. Norm equivalences.
Wavelet norm equivalences are akin to the classical Parseval relationin Fourier analysis (see [44, Section 12.1.2]) allowing to express Sobolev norms in terms of sumsof its Fourier coefficients. Norm equivalences are relevant for the construction of the mass- andstiffness matrices and for approximation estimates in the error analysis, cf. [13, Equations (2.5),and (2.6)]. In the unweighted univariate case there hold the standard norm equivalences(3.11) || u || H s ( I ) ≈ ∞ (cid:88) l =0 2 l (cid:88) k =0 ls | u lk | . For the stock price dimension V x we need norm equivalences in weighted spaces, for which therequirements ((3.1.2) and (3.1.2), cf. [13, Section 3]) are posed on the wavelets and their duals. Theorem 3.2 (Weighted norm equivalences) . Under the assumptions of Section 3.1.2 on thewavelet basis of the discretization spaces, there holds for any u ∈ L ( I, ω ) the norm equivalence ofits L ( I, ω ) -norm and of the discrete l ω -norm of its coefficients with respect to the wavelet basis. || u || L ( I,ω ) ≈ ∞ (cid:88) l =0 M l (cid:88) k =0 ω (2 − l k ) | u lk | (3.12) Proof.
See [13, Theorem 3.3], and also [13, Theorem 5.1]. (cid:3) With a view to error analysis, it is conventional (cf. [64, Section 3.1] and [44, Section 3.6.1]) to considerfunctions in V , which have additional regularity: In the unweighted univariate case one considers the classicalSobolev spaces H t ( I ), t = 0 , ,
2, where the subscript denotes Dirichlet boundary conditions.
IRICHLET FORMS AND FINITE ELEMENT METHODS FOR THE SABR MODEL 15
Bivariate setting.
For the bivariate case we set, similarly as above, without loss of generality G = (0 , × (0 , H t ( G, ω ), t = 0 , , H tx ( I, ω xt )and H ty ( I ) via tensor products, see Appendix A.2 for explicit constructions. We define the two-dimensional discretization spaces as the tensor product of the univariate discretization spaces forthe x and y coordinates(3.13) V L := V L x ⊗ V L y . Therefore, it holds similarly to (3.7) that(3.14) V L = span (cid:8) ψ l , k : 0 ≤ l x , l y ≤ L, ≤ k x ≤ M lx , ≤ k y ≤ M ly (cid:9) , where ψ l , k ( x, y ) = ψ ( l x ,l y ) , ( k x ,k y ) := ψ l x ,k x ( x ) ψ l y ,k y ( y ), for ( x, y ) ∈ (0 , × (0 , W l x = span { ψ l x , , . . . , ψ l x ,M lx } and W l y = span { ψ l y , , . . . , ψ l y ,M ly } are the correspondingcomplement spaces. Hence it follows from (3.8) with (3.13) directly that in the bivariate case(3.15) V L = (cid:18) (cid:77) ≤ l x ≤ L W l x (cid:19) ⊗ (cid:18) (cid:77) ≤ l y ≤ L W l y (cid:19) = (cid:77) ≤ l x ,l y ≤ L W l x ⊗ W l y . Every element u ∈ L ( G ) has a series representation cf. [44, Equation (13.6)](3.16) u = ∞ (cid:88) l x ,l y =0 M lx (cid:88) k x =1 M ly (cid:88) k y =1 u lk ψ l , k , where l , k = ( l x , l y ) , ( k x , k y ) , and the projection operator P L : V → V L , u (cid:55)→ P L u =: u L is defined as in (3.10) by truncating atlevel L the above series representation,(3.17) u L = L (cid:88) l x ,l y =0 M lx (cid:88) k x =1 M ly (cid:88) k y =1 u lk ψ l , k , u ∈ V, l , k = ( l x , l y ) , ( k x , k y ) . It is immediate from the underlying tensor product structure and from the norm-equivalences(3.11) of the one-dimensional case, that in the bivariate case there holds the norm equivalence(3.18) || u || H s ( G ) ≈ ∞ (cid:88) l x ,l y =0 2 lx (cid:88) k x =1 2 lx (cid:88) k x =1 (cid:0) s x l x + 2 s y l y (cid:1) | u lk | , for the Sobolev spaces H k ( G ), k = 0 , ,
2, where l , k = ( l x , l y ) , ( k x , k y ) , as in (3.9), and s = ( s x , s y )cf. [44, Equation 13.8]. For notational simplicity we stated here the unweighted version of thebivariate norm equivalence. Note however, that passing from the univariate to the bivariate caseis a direct consequence of the tensor product structure (see Section A.2) and is valid both inunweighted and weighted Sobolev spaces analogously.Now we are in a position to calculate the matrices M x , B x , and S x such as M y , B y , and S y as building blocks of the mass- and stiffness matrices M and A appearing in the semi-discreteproblem (3.3) with respect to the constructed wavelet basis { ψ l , k } : for the y -coordinate, thematrices in the appropriate weighted L ω -norm read(3.19) M yω y,M := (cid:18)(cid:82) ω y ( y ) ψ ly,ky ( y ) ψ l (cid:48) y,k (cid:48) y ( y ) ω (2 − ly k y ) ω (2 − l (cid:48) y k (cid:48) y ) dy (cid:19) ≤ l (cid:48) y ,l y ≤ L ; 0 ≤ k (cid:48) y ≤ l (cid:48) y , ≤ k y ≤ ly S yω y,S := (cid:18)(cid:82) ω y ( y ) ψ (cid:48) ly,ky ( y ) ψ (cid:48) l (cid:48) y,k (cid:48) y ( y ) ω (2 − ly k y ) ω (2 − l (cid:48) y k (cid:48) y ) dy (cid:19) ≤ l (cid:48) y ,l y ≤ L ; 0 ≤ k (cid:48) y ≤ l (cid:48) y , ≤ k y ≤ ly B yω y,B := (cid:18)(cid:82) ω y ( y ) ψ (cid:48) ly,ky ( y ) ψ l (cid:48) y,k (cid:48) y ( y ) ω (2 − ly k y ) ω (2 − l (cid:48) y k (cid:48) y ) dy (cid:19) ≤ l (cid:48) y ,l y ≤ L ; 0 ≤ k (cid:48) y ≤ l (cid:48) y , ≤ k y ≤ ly , Note that H = L ( G, x µ/ ) ⊂ L ( G ). and for the x -coordinate, the corresponding matrices are(3.20) M xω x,M := (cid:18)(cid:82) ω x ( x ) ψ lx,kx ( x ) ψ l (cid:48) x,k (cid:48) x ( x ) ω (2 − lx k x ) ω (2 − l (cid:48) x k (cid:48) x ) dx (cid:19) ≤ l (cid:48) x ,l x ≤ L ; 0 ≤ k (cid:48) x ≤ l (cid:48) x , ≤ k x ≤ lx S xω x,S := (cid:18)(cid:82) ω x ( x ) ψ (cid:48) lx,kx ( x ) ψ (cid:48) l (cid:48) x,k (cid:48) x ( x ) ω (2 − lx k x ) ω (2 − l (cid:48) x k (cid:48) x ) dx (cid:19) ≤ l (cid:48) x ,l x ≤ L ; 0 ≤ k (cid:48) x ≤ l (cid:48) x , ≤ k x ≤ lx B xω x,B := (cid:18)(cid:82) ω x ( x ) ψ (cid:48) lx,kx ( x ) ψ l (cid:48) x,k (cid:48) x ( x ) ω (2 − lx k x ) ω (2 − l (cid:48) x k (cid:48) x ) dx (cid:19) ≤ l (cid:48) x ,l x ≤ L ; 0 ≤ k (cid:48) x ≤ l (cid:48) x , ≤ k x ≤ lx , where ω x , and ω y denote weight functions in the respective dimensions. Our equations (3.19) and(3.20) closely follow the construction of [13, equation (3.12)] for the univariate weighted matrices.We take these as building blocks for the construction of our bivariate matrices, displayed in (3.22)and (3.23), with the appropriate choice of weight functions ( ω x,M = x µ , ω y,M = 1 in (3.22), andfurther ω y,M = e y , ω y,S = 1 , ω y,B = e y such as ω x,M = x µ , ω x,S = x β + µ , ω x,B = x β + µ in (3.23)).Similar constructions in the standard case can be found in [44]. In our case, the stiffness matrix A in the semi-discrete problem (3.3) is(3.21) A = (cid:0) A ( l (cid:48) , k (cid:48) ) , ( l , k ) (cid:1) ≤ l (cid:48) x ,l x ≤ L ; 0 ≤ k (cid:48) x ≤ l (cid:48) x , ≤ k x ≤ lx := ( a ( ψ l , k , ψ l (cid:48) , k (cid:48) )) ≤ l (cid:48) x ,l x ≤ L ; 0 ≤ k (cid:48) x ≤ l (cid:48) x , ≤ k x ≤ lx , where a ( · , · ) denotes the SABR-bilinear form (2.14). With respect to the weighted multiresolutionbasis { ψ l , k ψ l (cid:48) , k (cid:48) } defined in this section, the mass matrix reads(3.22) M = M xx µ ⊗ M y , and the stiffness matrix A takes the form A = (cid:0) Q xx S xx β + µ ⊗ M ye y + Q xy B xx β + µ ⊗ B ye y + Q yy M xx µ ⊗ S y (cid:1) + (cid:0) c x B xx β + µ − ⊗ M ye y + c x B xx β + µ ⊗ M ye y + c y M xx µ ⊗ B y (cid:1) , (3.23)where the coefficients are ( Q xx , Q xy , Q yy ) = ( , ρν, ν ) and ( c x , c x , c y ) = ( β + µ , ρν, ν ).3.2. Time discretization and the fully discrete scheme.
In this section we define a θ -schemefor our time discretization to introduce the fully discrete scheme of our finite element method.Furthermore, following [64] we conclude (see Proposition 3.4 below) that the stability of the θ -scheme remains valid in our setting. For this, we introduce the following dual norm for theapproximation spaces:(3.24) || f || ∗ := sup v L ∈ V L ( f, v L ) V ∗ ×V || v L || V , v L (cid:54) = 0 , f ∈ ( V L ) ∗ , furthermore, for T < ∞ and M ∈ N we consider the following uniform time-step and time mesh(3.25) k := TM , and t m = mk, m = 0 , . . . , M. The θ -scheme for the time-discretization and the fully discrete scheme are described as follows: Definition 3.3 ( θ -scheme and the fully discrete scheme) . Given the initial data u L := u (0 ,L ) = P L u , for the projector in (3.17) for m = 0 , . . . , M − u m +1 L ∈ V L such that for all v L ∈ V L :(3.26) k ( u m +1 L − u mL , v L ) V ∗ ×V + a ( θu m +1 L + (1 − θ ) u mL , v L ) = ( θg ( t m +1 ) + (1 − θ ) g ( t m ) , v L ) V ∗ ×V Hence, the fully discrete finite element scheme for the SABR model reads(3.27) ( k M + θ A ) u m +1 = k M u m − (1 − θ ) A u m + g m + θ , m = 0 , , . . . , M − , where M denotes the mass matrix (3.22), A the stiffness matrix (3.23), and u m is the coefficientvector of u mL with respect to the basis of V L . Proposition 3.4 (Stability of the θ -scheme) . For ≤ θ ≤ let the constants C and C satisfy (3.28) 0 < C < , C ≥ − C , IRICHLET FORMS AND FINITE ELEMENT METHODS FOR THE SABR MODEL 17 and for ≤ θ < denote λ A = sup v L ∈ V L || v L || H || v L || ∗ , and the constants C and C be such that (3.29) 0 < C < − σ, C ≥ − C ) σ − σ − C where σ := k (1 − θ ) λ A < . Then the sequence { u mL } Mm =0 of solutions of the θ -scheme 3.26 satisfy the stability estimate (3.30) || u ML || H + C k M − (cid:88) m =0 || u m + θk || H ≤ || u L || H + C k M − (cid:88) m =0 || g m + θ || ∗ . Proof.
The proof of this proposition is analogous to the one given in [64, Proposition 4.1]. (cid:3) Error estimates
Let V L , the finite dimensional approximation space of the solution space V , be as in Section 3.Furthermore, consider u m ( x ) := u ( t m , x ), with t m , m = 0 , . . . , M as in (3.25) and u mL as in thefully discrete scheme (3.26). In this section we estimate the error(4.1) e mL := u m ( x ) − u mL ( x ) = ( u m − P L u m ) + ( P L u m − u mL ) =: η m + ξ mL , for the the time-points t m , m = 0 , . . . , M , where P L : V → V L denotes the projection (3.10) on thefinite element space via truncated wavelet expansion. In Section 4.1 we derive estimates for theerror η m , m = 0 , . . . , M (approximation estimates). Section 4.2 is devoted to concluding from theresults of Section 4.1 corresponding error estimates for ξ mL , m = 0 , . . . , M and the convergence ofthe fully discrete scheme.4.1. Approximation estimates.
The crucial ingredient of our error analysis is the derivation ofapproximation estimates which measure the error η m = ( u m − P L u m ) between the true solution ofthe pricing equation at times t m , m = 0 , . . . , M and its projection to the discretization space in asuitably chosen norm. In the unweighted case, such approximation estimates—as in (4.2) below—in the usual H k ( I ) (resp. H k ( G )) norm, k = 0 , , V ⊂ H ⊂ V ∗ constructed in Section 2.2 and the corresponding discretization spaces in Section 3.1.4.1.1. Approximation estimates in the unweighted case.
In the unweighted univariate case it iswell-known that there exists for all u ∈ H l ( I ) with l = 0 , , u L in the correspondingdiscretization space V L , such that u L = P L u and for k = 0 , l = 0 , , l ≥ k it holds that(4.2) || u − P L u || H k ( I ) ≤ C − ( l − k ) L || u || H l ( I ) , where P L denotes the projector (3.10). The existence of such an element is provided by thenorm-equivalences (3.11). The same estimates hold in the two-dimensional case G = I × I ,see [44, Theorem 13.1.2.], in particular the approximation rate (2 − L ) ( l − k ) only depends on thediscretization level 2 − L and is independent of the dimension of the domain G . Analogously to theunivariate case, for k x = 0 , l x = 0 , , l x ≥ k x such as for k y = 0 , l y = 0 , , l y ≥ k y ,it holds that(4.3) || u − P L u || H k ( G ) ≤ (cid:40) C − ( l − k ) ∗ L || u || H l ( G ) , if k (cid:54) = 0 or l x , l y (cid:54) = 2 ,C − ( l − k ) ∗ L L / || u || H l ( G ) else , where we denote ( l − k ) ∗ := min { l x − k x , l y − k y } and where P L is the projection operator (3.17).The estimate (4.3) is a direct consequence of (4.2) and the tensor product construction (3.13).Analogous arguments obtain as above for bivariate norm-equivalences (3.18), cf. [44, Chapter 13]. Approximation estimates in the weighted case.
In the weighted setting, the order of approx-imation may depend on the norms of the weighted Sobolev spaces, in which we measure the error.In accord with usual conventions we consider functions in V with additional regularity for ourerror analysis in the weighted setting. For this purpose we consider weighted Sobolev spaces upto second order, H k ( G, ω ), k = 0 , ,
2, where the weight ω is yet to be chosen suitably to oursetting. We aim to establish approximation estimates of the following form: for any u ∈ H k ( G, ω ), k = 0 , , C > || u − P L u || H k ( G,ω ) ≤ (cid:40) C − c ω ( l − k ) ∗ L || u || H l ( G,ω ) , for k (cid:54) = 0 or l x , l y (cid:54) = 2 ,C − c ω ( l − k ) ∗ L L / || u || H l ( G,ω ) else , for a constant c ω ∈ R + , which may depend on the choice of the weights in H k ( G, ω ), k = 0 , , l − k ) ∗ = min { l x − k x , l y − k y } as in (4.3). In the following, we will first prove theone-dimensional version of the approximation property in weighted spaces. More specifically,we show statements analogous to (4.2) on the weighted Sobolev spaces H k ( I, x µ/ ), k = 0 , , x coordinate. For this, we pass for a function u ( t, x, y ) ∈ L ( J, V ) , t ∈ J, ( x, y ) ∈ G to u y ( t, x ) ∈ L ( J, V ) , t ∈ J, x ∈ I , for y ∈ I (see Appendix A.2)). Note that for the y coordinate, theunweighted setting prevails and the estimate (4.2) is valid. We then proceed to the bivariate case G = I × I by constructing tensor products of univariate multiresolution finite element spaces .Analogously as in the unweighted case (cf. equation (4.3)), the minimum of the obtained one-dimensional estimate then yields the estimate of (4.4), for the bivariate case H k ( G, ω ), k = 0 , , Definition 4.1.
Consider an interval I = (0 , R ), R > H kj ( I, x µ/ ) := { u : I → R measurable : D a u ∈ L ( I, x µ/ aβj ) , a ≤ k } , k = 0 , , , for j = 0 ,
1, with the norm || u || H kj ( I,x µ/ ) := (cid:88) a ≤ k (cid:90) I | D a u ( x ) | x µ + aβj dx, k = 0 , , . (4.6)To ease notation we shall henceforth denote H kj =0 by H k and H kj =1 by H k , k = 0 , , Remark 4.2.
Note that for the spaces H and V in Remark 2.8 and for the spaces H kj , k = 0 , , j = 0 , H = H ((0 , R ) , x µ/ ) = H ((0 , R ) , x µ/ ) and the weighted space V satisfies V = H ((0 , R ) , x µ/ ) and V ⊃ H ((0 , R ) , x µ/ ) such as the estimate(4.7) || v || V ≤ C R || v || H ((0 ,R ) ,x µ/ ) , v ∈ V for any finite R > C R > Proposition 4.3.
The projection operator in (3.10) satisfies for k = 0 , and l = 0 , , , l ≥ k theestimate (4.8) || u − P L u || H k ( I,x µ/ ) ≤ C − ( l − k ) L || u || H l ( I,x µ/ ) , u ∈ H ( I, x µ/ ) . It is immediate from Remark 4.2, that the approximation estimate (4.8) readily applies tothe CEV model. Furthermore, combining estimate for the x coordinate with the unweightedestimate (4.2) for the y coordinate and taking tensor products (see Appendix A.2.1 for an explicitconstruction of the bivariate spaces) yields that the approximation estimate (4.3) remains valid inthe (weighted) bivariate case. In contrast to this, measuring the error in the norms || u || H kj =1 ( I,x µ/ ) , k = 0 , ,
2, the approximation in the x coordinate dominates the y coordinate, see Remark 4.4. Remark 4.4.
If we consider j = 1 in Definition 4.1, we do not assume any additional integrabilityrequirements on our solution up to first order derivatives. In this case we obtain (weaker) approx-imation estimates where the order of approximation depends on the parameter β : the projection Cf. [64, Section 3.1] and [44, Section 3.6.1] and see also the unweighted case above. See Appendix A.2.1 for an explicit construction.
IRICHLET FORMS AND FINITE ELEMENT METHODS FOR THE SABR MODEL 19 operator P L : V → V L in (3.10) satisfies for k = 0 , l = 0 , , l ≥ k the estimate(4.9) || u − P L u || H k ( I,x µ/ ) ≤ C − (1 − β )( l − k ) L || u || H l ( I,x µ/ ) , u ∈ H ( I, x µ/ ) . A proof of this remark is delegated to the Appendix C.1.
Proof of Proposition 4.3.
Let u ∈ V and consider P L u ∈ V L . Then it is immediate from (3.9) and(3.10) that u − P L ( u ) = ∞ (cid:88) l = L +1 2 l (cid:88) j =1 u lj ψ l,j . It directly follows from the norm equivalence (3.11) and the scaling of the wavelet basis-elementsto unit norm in L ( I ) that the derivatives satisfy || ( u − P L u ) (cid:48) || L ( I,x µ/ ) = ˜ C ∞ (cid:88) l = L +1 l l (cid:88) k =0 ω (2 − l k ) | u lk | ≥ ˜ C L ∞ (cid:88) l =0 2 l (cid:88) k =0 ω (2 − l k ) | u lk | = ˜ C L || u − P L u || L ( I,x µ/ ) . (4.10)Now with C = C and recalling that we have set h = 2 − L , the following relation is obtained:(4.11) || u − P L u || L ( I,x µ/ ) ≤ Ch || ( u − P L u ) (cid:48) || L ( I,x µ/ ) ≤ Ch || u (cid:48) || L ( I,x µ/ ) ≤ Ch (cid:16) || u || L ( I,x µ/ ) + || u (cid:48) || L ( I,x µ/ ) (cid:17) = Ch || u || H ( I,x µ/ ) . Analogously, replacing ( u − P L u ) (cid:48) by ( u − P L u ) (cid:48)(cid:48) and ( u − P L u ) by ( u − P L u ) (cid:48) in equations (4.10)and (4.11), one obtains(4.12) || ( u − P L u ) (cid:48) || L ( I,x µ/ ) ≤ Ch || ( u − P L u ) (cid:48)(cid:48) || L ( I,x µ/ ) ≤ Ch || u (cid:48)(cid:48) || L ( I,x µ/ ) , and adding up (4.11) and (4.12) yields:(4.13) || u − P L u || H ( I,x µ/ ) = || u − P L u || L ( I,x µ/ ) + || ( u − P L u ) (cid:48) || L ( I,x µ/ ) ≤ Ch (cid:80) k =0 || u ( k ) || L ( I,x µ/ ) = Ch || u || H ( I,x µ/ ) . Finally, concatenating (4.10) and (4.11) directly yields:(4.14) || u − P L u || L ( I,x µ/ ) ≤ Ch || ( u − P L u ) (cid:48) || L ( I,x µ/ ) ≤ C (2 − L ) || ( u − P L u ) (cid:48)(cid:48) || L ( I,x µ/ ) ≤ C (2 − L ) || u (cid:48)(cid:48) || L ( I,x µ/ ) ≤ C (2 − L ) || u || H ( I,x µ/ ) . (cid:3) Discretization error and convergence of the finite element method.
In this sectionwe apply the estimates of Section 4.1 to derive estimates on the discretization error (cf. equa-tion (4.1)) and to conclude the convergence of the proposed finite element approximation of thevariational solution of the SABR pricing equations. For the estimates of the discretization errorwe follow the (unweighted) analysis of [64, Section 5]. Corresponding proofs prevail with minormodifications and are provided—accommodated to our setting and notations—in the Appendix Cfor easy reference.
Lemma 4.5.
For u ∈ C ( ¯ J ; H ( G, a )) , the errors ξ mL are the solutions of the θ -scheme:Given ξ L := P L u − u L , for m = 0 , . . . , M − find ξ m +1 L ∈ V L such that for all v L ∈ V L : k ( ξ m +1 L − ξ mL , v L ) V×V ∗ + a ( θξ m +1 L + (1 − θ ) ξ mL , v L ) =: ( r m , v L ) V×V ∗ (4.15)The proof of the above Lemma relies on the observation that the errors ξ satisfy the samePDE as the solutions u , therefore the θ schemes for ξ and for u are analogous . The followingcorollary is a direct consequence of Lemma 4.5 and of the the stability of the θ -scheme, establishedin Proposition 3.4. See Appendix C for details, where we also included for completeness a reminder of the proof of [64, Lemma5.1]—accommodated to our notation—which directly carries over to the present situation.
Corollary 4.6.
There exist constants C and C independent from the discretization level L andtime mesh k such that the solutions of (4.15) satisfy the estimate (4.16) || ξ ML || H + C k M − (cid:88) m =0 || ξ m + θL || a ≤ || ξ N || H + C k M − (cid:88) m =0 || r m || ∗ , where for any f ∈ ( V L ) ∗ , || f || ∗ := sup v L ∈V L ( f,v L ) V×V∗ || v L || a , cf. (3.24) and where r m , m = 0 , . . . , M − denote the weak residuals defined in equation (4.15) . The weak residual r m , m = 0 , . . . , M − r m , v L ) V×V ∗ := ( r m , v L ) V×V ∗ + ( r m , v L ) V×V ∗ + a ( r m , v L ) , where the components r m , r m and r m , m = 0 , . . . , M are defined as( r m , v L ) V×V ∗ := ( k ( u m +1 − u m ) − ˙ u m + θ , v L ) V×V ∗ , ( r m , v L ) V×V ∗ := ( k ( P L u m +1 − P L u m ) + k ( u m +1 − u m ) , v L ) V×V ∗ , ( r m , v L ) V×V ∗ := a ( P L u m + θ − u m + θ , v L ) . (4.17)Using this decomposition facilitates the following estimates for the residuals. Lemma 4.7 (Norm estimates for the residuals) . Consider the weak residuals r m of the θ -scheme (4.15) for m = 0 , . . . , M − . Furthermore, assume that u ∈ C ( ¯ J ; H j ( G, x µ/ )) ∩ C ( J ; H j ( G, x µ/ )) , j = 0 , , where H kj ( G, x µ/ ) , k = 0 , , , j = 0 , are the spaces in (A.9) . Then there holds the estimate (4.18) || r m || ∗ ≤ C k / (cid:16)(cid:82) t m +1 t m || ¨ u || ∗ ds (cid:17) / , θ ∈ [0 , k / (cid:16)(cid:82) t m +1 t m || ... u || ∗ ds (cid:17) / , θ = +2 − L (cid:18) Ck / (cid:16)(cid:82) t m +1 t m || ˙ u || H j ( G,x µ/ ) ds (cid:17) / + C || u m + θ || H j ( G,x µ/ ) (cid:19) , A proof of this Lemma is provided in the Appendix C. With these preparations, we are in aposition to prove the main result of this section:
Theorem 4.8 (Convergence of the finite element approximation: SABR) . Assume that u ∈ C ( ¯ J ; H j ( G, x µ/ )) ∩ C ( J ; H j ( G, x µ/ )) , j = 0 , , where H kj ( G, x µ/ ) , k = 0 , , , j = 0 , are the spaces in (A.9) , and assume further that theapproximation u (0 ,L ) ∈ V L of the initial data is quasi optimal, that is || ξ L || H = || u − u (0 ,L ) || H ≤ C − L || u || H . (4.19) Let u m ( z ) = u ( t m , z ) , z ∈ G for t m , m = 0 , . . . M be as in (3.25) let u mL denote the solution ofthe fully discrete scheme (3.26) , and let the approximation space V L be as in Section 3. Then, thefollowing error bounds hold: || u M − u ML || H + k M − (cid:88) m =0 || u m + θ − u m + θL || a ≤ C − j (1 − β )2 L max ≤ t ≤ T || u ( t ) || H j + 2 − j (1 − β )2 L (cid:90) T || ˙ u ( s ) || H j ds + C (cid:40) k (cid:82) T || ¨ u ( s ) || ∗ ds, ≤ θ ≤ k (cid:82) T || ... u ( s ) || ∗ ds, θ = for j = 0 , , where u m + θL = θu m +1 L + (1 − θ ) u mL , and u m + θ = θu m +1 + (1 − θ ) u m . Remark 4.9 (Convergence of the finite element approximation: CEV) . The same estimates holdfor the CEV model, when we replace in Theorem 4.8 the spaces H kj ( G, x µ/ ), k = 0 , , j = 0 , H kj ( I, x µ/ ), k = 0 , , j = 0 , IRICHLET FORMS AND FINITE ELEMENT METHODS FOR THE SABR MODEL 21
Proof of Theorem 4.8, and Remark 4.9.
For the proofs we follow [44, Theorem 3.6.5], and [64,Theorem 5.4.] with the appropriate adjustments. For brevity we shall prove both cases H j =0 and H j =1 together, and denote the generic convergence of order by 2 − Lc ω := 2 − L − Lj (1 − β ) as in(4.4), where c ω = 1 for H j =0 and c ω = (1 − β ) for H j =1 . By Corollary 4.6, || e ML || H = || u m ( x ) − u mL ( x ) || H = || η m + ξ mL || H ≤ (cid:16) || η m || H + k M − (cid:88) m =0 || η m + θ || a + || ξ mL || H + k M − (cid:88) m =0 || ξ m + θL || a (cid:17) . This yields || e ML || H + C k M − (cid:88) m =0 || e m + θ || a ≤ (cid:16) || η M || H + C k M − (cid:88) m =0 || η m + θ || a + || ξ mL || H + C k M − (cid:88) m =0 || ξ m + θL || a (cid:17) ≤ C (cid:16) || η M || H + k M − (cid:88) m =0 || η m + θ || a + || ξ L || H + C k M − (cid:88) m =0 || r m || ∗ (cid:17) , and by Lemma 4.7 one can further estimate the last terms to obtain C (cid:16) || η M || H + k M − (cid:88) m =0 || η m + θ || a + || ξ L || H + C k M − (cid:88) m =0 || r m || ∗ (cid:17) ≤ C (cid:40) || η M || H + k M − (cid:88) m =0 || η m + θ || a + || ξ L || H + C M − (cid:88) m =0 (2 − L ) c ω (cid:32)(cid:90) t m +1 t m || ˙ u || H j ds + || u m + θ || H j + (cid:40) k (cid:82) t m +1 t m || ¨ u || ∗ ds, θ ∈ [0 , k (cid:82) t m +1 t m || ... u || ∗ ds, θ = (cid:33) (cid:41) ≤ C || ξ L || H + C (2 − L ) c ω (cid:90) T || ˙ u || V ds + C max ≤ t ≤ T (2 − L ) c ω || u ( t ) || H j + C (cid:40) k (cid:82) T || ¨ u || ∗ ds, θ ∈ [0 , k (cid:82) T || ... u || ∗ ds, θ = , where the last step follows by || · || a ≤ || · || V ≤ C G || · || H j ( G,x µ/ ) for a C G > (cid:3) Numerical Experiments
In this section we carry out numerical experiments for option prices to investigate the robustnessof the derived finite element discretisation in a simple linear setup described in [44, Chapter 4].We find that the method performs robustly even for long maturities (which is precisely where theshort-time asymptotic formula of Hagan et al [43] is known to break down), for different regimesof the CEV parameter (both for β < . β ≥ .
5) and throughout correlations. We displaytwo examples here and remark that the convergence we find in practice (readily for this bases)outperforms the theoretically predicted ones.
The image on the left gives finite element option prices for β = 0 . ν = 1 in the uncorrelatedcase, with maturity T = 25 years and K = 1, while the image on the right shows the convergence(with respect to the mesh-width h = 2 L ) of the finite element approximation in a linear basis.Here, the image on the left gives finite element option prices for β = 0 . ν = 1 and correlation ρ = − . T = 10 years, while the image on the right shows the convergence (withrespect to the mesh-width h = 2 L ) of the finite element approximation in a linear basis.The last images display put options as an approximation of the mass at zero on the time horizon T = 10. The put options are of the form max(1 − (cid:15) X,
0) for a small (cid:15) >
0. The parameters inthis approximation are β = 0 . , ρ = 0, and ν = 1, and we chose throughout µ = − β . Appendix A. Reminder on properties of considered function spaces
IRICHLET FORMS AND FINITE ELEMENT METHODS FOR THE SABR MODEL 23
A.1.
Weighted spaces.
To make the reading self-contained, we include a reminder on someproperties of weighted Sobolev spaces used in this article and refer the reader to the monographsof [51] and [16] for full details. By a weight, we shall mean a locally integrable function ω on R such that ω ( x ) > ω gives rise to a measure (via integration ω ( E ) = (cid:82) E ω ( x ) dx ,for measurable sets E ⊂ R ). This measure is also denoted by ω . Definition A.1 (Weighted L -space) . Let ω be a weight on an open set G ⊂ R . L ( G, ω ) is theset of measurable functions u on G such that(A.1) || u || L ( G,ω ) = (cid:90) G | u ( x ) | ω ( x ) dx < ∞ Definition A.2 (Weighted Sobolev space) . Let k ∈ N and Let a = { ω a = ω a ( x ) , x ∈ G, | a | ≤ k } be a given family of weight functions on an open set G ⊂ R . We denote by W k ( G, ϕ ) the set ofall functions u ∈ L ( G, ω ) for which the weak derivatives D ( a ) u , with | a | ≤ k , belong to L ( G, ω a ).The weighted Sobolev space W k ( G, ϕ ) is a normed linear space if equipped with the norm(A.2) || u || W k ( G,ω ) = (cid:88) | a |≤ k (cid:90) G | D a u ( x ) | ω a ( x ) dx. Remark A.3. If ω ω ∈ L loc ( G ) then C ∞ ( G ) is a subset of W k ( G, ω ), and we can introduce thespace W k ( G, ω ) as the closure of C ∞ ( G ) with respect to the norm W k ( G, ω ), see also [16, 52].The class of A p weights was introduced by B. Muckenhoupt (cf. [61]). This class is relevant forthe property that it ensures C ∞ ( G ) ⊂ W k ( G, ω ). This is used in Sections 2.1 and 2.2.A weight ω is in A p if there exists a positive constant C such that for every ball B ⊂ R (A.3) (cid:18) | B | (cid:90) B ωdx (cid:19) (cid:18) | B | (cid:90) B ω − / ( p − dx (cid:19) p − ≤ C. Lemma A.4. ω ( x ) := | x | , x ∈ R is in A p if and only if − < ω < p − .Proof. See Corollary 4.4 in [72], and Corollary 2.18 in [35]. (cid:3)
Lemma A.5. If ω ∈ A p then since ω − / ( p − is locally integrable, we have L p ( G, ω ) ⊂ L loc ( G ) for every open set G ⊂ R n . Therefore, weak derivatives of functions in L p ( G, ω ) are well-defined.Furthermore, if ω ∈ A p then C ∞ ( G ) is dense in W k,p ( G, ω ) .Proof. See [73, Corollary 2.1.6 ] and [31, Theorem 1.5]. (cid:3)
A.2.
Tensor Products of Hilbert Spaces.
Let I := (0 ,
1) denote the unit interval and G := I × I . See in [44, Section 13.1] that the Hilbert spaces H k ( G ), k = 0 , , H k ( I ) via the tensor product structure: L ( G ) ∼ = (cid:0) L ( I ) ⊗ L ( I ) (cid:1) (A.4) H ( G ) ∼ = (cid:0) H ( I ) ⊗ L ( I ) (cid:1) (cid:92) (cid:0) L ( I ) ⊗ H ( I ) (cid:1) (A.5) H ( G ) ∼ = (cid:0) H ( I ) ⊗ L ( I ) (cid:1) (cid:92) (cid:0) H ( I ) ⊗ H ( I ) (cid:1) (cid:92) (cid:0) L ( I ) ⊗ H ( I ) (cid:1) (A.6)Recall from [67, Chapter II. 4] that the inner products on each of the tensor-Hilbert spaces aredefined as(A.7) (cid:104) u ⊗ u , v ⊗ v (cid:105) H ⊗ H := (cid:104) u , v (cid:105) H (cid:104) u , v (cid:105) H , for u , v ∈ H and u , v ∈ H , where H , H stand for generic Hilbert spaces (say any of thetensor products involved in (A.4),(A.5), or (A.6) above). The inner products on the intersectionspaces H k ( G ), k = 0 , , u ≡ u x ⊗ u y , and v ≡ v x ⊗ v y ∈ L ( G ) for u x , v x ∈ L ( I ) and u y , v y ∈ L ( I )we recall the following [67, Theorem II. 10. c), Chapter II. 4]: Theorem A.6.
Let ( M , µ ) and ( M , µ ) be measure spaces so that L ( M , µ ) and L ( M , µ ) are separable, then there is a unique isomorphism such that (A.8) L ( M × M , dµ ⊗ dµ ) (cid:55)−→ L ( M , dµ ; L ( M , dµ )) f ( x, y ) (cid:55)−→ ( x (cid:55)→ f ( x, · )) . A.2.1.
Explicit construction of the bivariate spaces in Section 4.1.2.
Consider our domain of in-terest G = (0 , R x ) × ( − R y , R y ), R x , R y > H kj ( G, x µ/ ) := { u : G → R measurable : ∂ | a | a u ∈ L ( I, x µ/ a x βj ) , a ≤ k } , k = 0 , , , for j = 0 ,
1, and a multiindex a with | a | = a x + a y , where a x denotes the number of derivatives indirection x and a y in direction y . The respective norms in H kj ( G, x µ/ ) for j = 0 are defined by(A.10) || u || H j =0 ( G,x µ/ ) := || x µ/ u || L || u || H j =0 ( G,x µ/ ) := || x µ/ u || L + || x µ/ ∂ y ( u ) || L + || x µ/ ∂ x ( u ) || L || u || H j =0 ( G,x µ/ ) := || x µ/ u || L + || x µ/ ∂ y ( u ) || L + || x µ/ ∂ x ( u ) || L + || x µ/ ∂ yy u || L + || x µ/ ∂ xy u || L + || x µ/ ∂ yy u || L , and for j = 1 by(A.11) || u || H j =1 ( G,x µ/ ) := || x µ/ u || L || u || H j =1 ( G,x µ/ ) := || x µ/ u || L + || x µ/ ∂ y ( u ) || L + || x β + µ/ ∂ x ( u ) || L || u || H j =1 ( G,x µ/ ) := || x µ/ u || L + || x µ/ ∂ y ( u ) || L + || x β + µ/ ∂ x ( u ) || L + || x µ/ ∂ yy u || L + || x β + µ/ ∂ xy u || L + || x β + µ/ ∂ yy u || L . Remark A.7.
Note, that similarly as in Remark 4.2, the spaces H and V in Definitions 2.2and 2.4 and the spaces H kj ( G ), k = 0 , , j = 0 , H = H j =0 ( G, x µ/ ) = H j =1 ( G, x µ/ ), and the weighted space V satisfies V = H j =1 ( G, x µ/ ) and V ⊃ H j =0 ( G, x µ/ ), furthermore, on any bounded domain G there holds the estimate(A.12) || v || V ≤ C G || v || H j =0 ( G,x µ/ ) , v ∈ V , C G > . Lemma A.8.
The spaces in (A.9) , k = 0 , , , j = 0 , can be constructed as tensor products ofthe spaces (4.5) and the usual (unweighted) Sobolev spaces H k ( G ) , k = 0 , , , j = 0 , via (A.4)(A.5) and (A.6) as follows H j ( G, x µ/ ) ∼ = (cid:16) H j ( I, x µ/ ) ⊗ H ( I ) (cid:17) H j ( G, x µ/ ) ∼ = (cid:16) H j ( I, x µ/ ) ⊗ H ( I ) (cid:17) (cid:92) (cid:16) H j ( I, x µ/ ) ⊗ H ( I ) (cid:17) H j ( G, x µ/ ) ∼ = (cid:16) H j ( I, x µ/ ) ⊗ H ( I ) (cid:17) (cid:92) (cid:16) H j ( I, x µ/ ) ⊗ H ( I ) (cid:17) (cid:92) (cid:16) H j ( I, x µ/ ) ⊗ H ( I ) (cid:17) . Appendix B. Non-symmetric Dirichlet forms
We include here a condensed reminder of some of the basic concepts on non-symmetric Dirichletforms, which are used in previous sections, in particular in Theorem 2.19. For full details see [70].
Definition B.1 (Symmetric closed form) . A pair ( E , D ( E )) is called a symmetric closed form onthe Hilbert space ( H , ( · , · ) H ), if D ( E ) is a dense linear subspace of H , and E : D ( E ) × D ( E ) → R isa non-negative definite symmetric bilinear form, which is closed on H . That is, D ( E ) is a completemetric space with respect to the norm E ( · , · ) / := ( E ( · , · ) + ( · , · ) H ) / . Definition B.2 (Coercive closed form) . A pair ( E , D ( E )) is called a coercive closed form on theHilbert space H if D ( E ) is a dense linear subspace of H , and E : D ( E ) × D ( E ) → R is a bilinearform such that the following two conditions hold: • Its symmetric part (cid:101) E ( u, v ) := ( E ( u, v ) + E ( v, u )) is a symmetric closed form on H . IRICHLET FORMS AND FINITE ELEMENT METHODS FOR THE SABR MODEL 25 • The pair ( E , D ( E )) satisfies the so-called weak sector condition : there exists a continuityconstant K > |E ( u, v ) | ≤ K E ( u, u ) / E ( v, v ) / for all u, v ∈ D ( E ) . Remark B.3.
Recall the continuity property (2.17) (also considered in Section 2.2) |E ( u, v ) | ≤ K E ( u, u ) / E ( v, v ) / for all u, v ∈ D ( E )implies the weak sector condition (B.1) above. Definition B.4 (Dirichlet form) . Consider a Hilbert space ( H , ( · , · ) H ) of the form H = L ( E, m ),where (
E, m ) is a measure space. A coercive closed form ( E , D ( E )) on L ( E, m ) is called a (non-symmetric) Dirichlet form , if for all u ∈ D ( E ) ⊂ E , one has the contraction properties (B.2) u + ∧ ∈ D ( E ) and E ( u + u + ∧ , u − u + ∧ ≥ E ( u − u + ∧ , u + u + ∧ ≥ , where for any u, v : E → R , we have set(B.3) u ∧ v := inf( u, v ) , u ∨ v := sup( u, v ) , u + := u ∨ , u − := − ( u ∧ . A coercive closed form satisfying one of the two inequalities in (B.2) is called - Dirichlet form .If ( E , D ( E )) is in addition symmetric, that is E = (cid:101) E , where (cid:101) E denotes the symmetric part of E (recall (cid:101) E ( u, v ) := ( E ( u, v ) + E ( v, u ))), then ( E , D ( E )) is called a symmetric Dirichlet form.In the latter case, the contraction property in condition (B.2) reduces to(B.4) E ( u + ∧ , u + ∧ ≤ E ( u, u ) . See [70, Section 4, Def. 4.5]
Theorem B.5.
Let ( E , D ( E )) be a coercive closed form on a Hilbert space ( H , ( · , · ) H ) with conti-nuity constant K > . Define the domain (B.5) D ( A ) := { u ∈ D ( E ) | v (cid:55)→ E ( u, v ) is continuous w.r.t. ( · , · ) / H on D ( E ) } . For any u ∈ D ( A ) , let Au denote the unique element in H such that (B.6) ( − Au, v ) = E ( u, v ) for all v ∈ D ( E ) . Then A is the generator of the unique strongly continuous contraction resolvent ( G α ) α> on H which satisfies (B.7) G α ( H ) , ⊂ D ( E ) and E ( G α f, u ) + α ( G α f, u ) H = ( f, u ) H for all f ∈ H , u ∈ D ( E ) , α > . Furthermore, since ( E , D ( E )) is a coercive closed form on ( H , ( · , · ) H ) , there exists a further uniquestrongly continuous contraction resolvent ( ˆ G α ) α> on H , which satisfies (B.8) ˆ G α ( H ) ⊂ D ( E ) and ( f, u ) H = E ( u, ˆ G α f ) + α ( u, ˆ G α f ) H for all f ∈ H , u ∈ D ( E ) , α > . In particular, ˆ G α is the adjoint of G α for all α > . That is (B.9) ( G α f, g ) H = ( f, ˆ G α g ) H for all f, g ∈ H , and similarly, for the (unique) strongly continuous contraction semigroups ( P t ) t ≥ , ( ˆ P t ) t ≥ corre-sponding to ( G α ) α> and ( ˆ G α ) α> respectively it holds that (B.10) ( P t f, g ) H = ( f, ˆ P t g ) H for all f, g ∈ H , t ≥ . Proof.
See: [70, Theorem 2.8, Corollary 2.10 and Proposition 2.16]. (cid:3)
Definition B.6 (Contraction Properties) . Let ( H , ( · , · ) H ) be a Hilbert space where H = L ( E, m ),and where (
E, m ) is a measure space. For any f, g ∈ L ( E, m ) we write f ≤ g or f < g forany m -classes f, g of functions on E , if the respective inequality holds m -a.e. for correspondingrepresentatives. See [70, Section 4]. ( i ) Let G be a bounded linear operator on L ( E, m ) with domain D ( G ) = L ( E, m ). Then G is called sub-Markovian , if for all f ∈ L ( E, m ) the condition 0 ≤ f ≤ ≤ Gf ≤ ii ) A strongly continuous contraction semigroup ( P t ) t ≥ resp. resolvent ( G α ) α> is called sub-Markovian if all P t , t ≥ αG α , α > iii ) A closed, densely defined operator A on ( L ( E, m ) , ( · , · ) H ) is called Dirichlet operator if(
Au, ( u − + ) H ≤ u ∈ D ( A ) ⊂ E .See [70, Section 4, Def. 4.1]. Theorem B.7.
Consider a Hilbert space ( H , ( · , · ) H ) of the form H = L ( E, m ) , where ( E, m ) isa measure space. Let ( G α ) α> be a strongly continuous contraction resolvent on ( L ( E, m ) , ( · , · ) H ) with corresponding generator A and semigroup ( P t ) t ≥ . Furthermore, let ( E , D ( E )) be a coerciveclosed form on L ( E, m ) with continuity constant K > and corresponding resolvent ( G α ) α> .Then the following are equivalent: ( i ) ( P t ) t ≥ is sub-Markovian. ( ii ) A is a Dirichlet operator. ( iii ) ( G α ) α> is sub-Markovian. ( iv ) for all u ∈ D ( E ) , u + ∧ ∈ D ( E ) and E ( u + u + ∧ , u − u + ∧ ≥ , that is ( E , D ( E )) is a -Dirichlet form.If in the above statements the operators ( G α ) α> (resp. ( P t ) t ≥ and A ) are replaced by theiradjoints ( ˆ G α ) α> (resp. ( ˆ P t ) t ≥ and ˆ A ), then the analogous equivalences hold, where in ( iv ) theentries of E are interchanged. Hence, if ( iii ) (resp. ( ii ) or ( i ) ) holds both for ( G α ) α> (resp. A or ( P t ) t ≥ ) and its adjoint ( ˆ G α ) α> (resp. ˆ A or ( ˆ P t ) t ≥ ), then the coercive closed form ( E , D ( E )) is a (non-symmetric) Dirichlet form.Proof. See [70, Section 4 Proposition 4.3 and Theorem 4.4]. (cid:3)
Appendix C. Further Proofs
Proof of Lemma 4.5.
The statement of the Lemma is by the decomposition (4.17) essentially aconsequence of the fact that the errors ξ satisfy the same PDE as the functions u : The variationalformulation implies(C.1) ( ˙ u m + θ , v ) V×V ∗ + a ( u m + θ , v ) = ( g m + θ , v ) V×V ∗ ∀ v ∈ V . Using the definition (4.1) of ξ , we rewrite (4.15) in the form: k ( ξ m +1 L − ξ mL , v L ) V×V ∗ + a ( θξ m +1 L + (1 − θ ) ξ mL , v L ) = k (( P L u m +1 − u m +1 L ) − ( P L u m − u mL ) , v L ) V×V ∗ + a ( P L u m + θ + u m + θL , v L ) = (cid:16) P L u m +1 − P L u m k , v L (cid:17) V×V ∗ + a ( P L u m + θ , v L ) − (cid:0) k ( u m +1 L − u mL , v L ) V×V ∗ − a ( u m + θL , v L ) (cid:1) where we used u m + θ := θu m +1 + (1 − θ ) u m and the linearity of the projector: θP L u m +1 − (1 − θ ) P L u m +1 L = P L u m + θ . Furthermore, by the θ -scheme (3.26) for u , and by (C.1) (cid:16) P L u m +1 − P L u m k , v L (cid:17) V×V ∗ + a ( P L u m + θ , v L ) − (cid:18)(cid:16) u m +1 L − u mL k , v L (cid:17) V×V ∗ − a ( u m + θL , v L ) (cid:19) = (cid:16) P L u m +1 − P L u m k , v L (cid:17) V×V ∗ + a ( P L u m + θ , v L ) − ( g m + θ , v L ) V×V ∗ = (cid:16) P L u m +1 − P L u m k , v L (cid:17) V×V ∗ + a ( P L u m + θ , v L ) − a ( u m + θ , v L ) − ( ˙ u m + θ , v L ) V×V ∗ = (cid:16) P L u m +1 − P L u m k − u m +1 − u m k , v L (cid:17) V×V ∗ + a (cid:0) P L u m + θ − u m + θ , v L (cid:1) + (cid:16) u m +1 − u m k − ˙ u m + θ , v L (cid:17) V×V ∗ = ( r , v L ) V×V ∗ + ( r , v L ) V×V ∗ + ( r , v L ) V×V ∗ . (cid:3) IRICHLET FORMS AND FINITE ELEMENT METHODS FOR THE SABR MODEL 27
Proof of Lemma 4.7.
We adapt [44, Section 3.6.2] and [64, Lemma 5.3] to the situation at handand confirm that the estimates of [64, Lemma 5.3] carry over to the weighted case. Analogouslyto the classical (unweighted) case, the statement of the Lemma follows from || r m || ∗ ≤ || r m || ∗ + || r m || ∗ + || r m || ∗ and the corresponding norm estimates for the decomposition (4.17). The estimateof the residual r is(C.2) || r || ∗ = || k ( u m +1 − u m ) − ˙ u m + θ || ∗ ≤ k (cid:16)(cid:82) t m +1 t m | s − (1 − θ ) t m +1 − θt m | || ¨ u || ∗ ds (cid:17) ≤ C θ k / (cid:16)(cid:82) t m +1 t m || ¨ u || ∗ ds (cid:17) / . In case θ = , partial integration yields the refined estimate(C.3) || r || ∗ = || k ( u m +1 − u m ) − ˙ u m + θ || ∗ ≤ k (cid:16)(cid:82) t m +1 t m | ( t m +1 − s )( t m − s ) | || ... u || ∗ ds (cid:17) ≤ Ck / (cid:16)(cid:82) t m +1 t m || ... u || ∗ ds (cid:17) / . The norm of the residual r is bounded by(C.4) || r || ∗ ≤ − L Ck / (cid:16)(cid:82) t m +1 t m || ˙ u || V ds (cid:17) / . The bound (C.4) follows from (3.24) and from the estimate(C.5) | ( r , v L ) V×V ∗ | = | ( k ( P L u m +1 − P L u m ) + k ( u m +1 − u m ) , v L ) V×V ∗ |≤ C || k ( P L u m +1 − P L u m ) + k ( u m +1 − u m ) || ∗ || v L || a ≤ Ck || ( I − P L ) (cid:82) t m +1 t m ˙ u ( s ) ds || ∗ || v L || a ≤ Ck / (cid:16)(cid:82) t m +1 t m || ˙ u − P L ˙ u || H ds (cid:17) / || v L || a (4.4) ≤ Ck / (2 − L ) c a (cid:16)(cid:82) t m +1 t m || ˙ u || H j ds (cid:17) / || v L || a , where the last step follows from the approximation property (4.8) resp. (4.9) and the estimate(A.12). The norm of the residual r allows for the upper bound(C.6) || r || ∗ ≤ C (2 − L ) c a || u m + θ || H j . The estimate (C.6) follows from (3.24) and from(C.7) | ( r , v L ) V×V ∗ | = | a ( P L u m + θ − u m + θ , v L ) |≤ || P L u m + θ − u m + θ || a || v L || a (4.4) ≤ C (2 − L ) c a || u m + θ || H j || v L || a . The second step follows from a simple polarisation argument | ( u, v ) V×V ∗ | ≤ (cid:12)(cid:12) ( a ( u + v, u + v ) − a ( u − v, u − v )) (cid:12)(cid:12) ≤ (cid:12)(cid:12) ( a ( u, u ) + a ( v, v )) (cid:12)(cid:12) ≤ | ( a ( u, u ) a ( v, v )) | for u = P L u m + θ − u m + θ and v = v L , and in the last step we used || · || a ≤ || · || V (which holds bycontinuity of the bilinear form cf. Lemma 2.14) and the approximation property (4.8) resp. (4.9). (cid:3) C.1.
Approximation estimates for CEV: Alternative.
If we do not assume any additionalintegrability requirements on our solution up to first order derivatives, we stated in Remark 4.4that we obtain such approximation estimates where the order of approximation depends on theparameter β . The proof of the estimates (4.9) in Remark 4.4 is similar to that of Proposition 4.3and is included here: Proof of Remark 4.4.
For j = 1 and µ ∈ [max {− , − β } , − β ] the inequalities in (4.10) become || ( u − P L u ) (cid:48) || L ( I,x µ/ β ) = ˜ C ∞ (cid:88) l = L +1 l l (cid:88) k =0 (2 − l k ) µ +2 β | u lk | ≥ ˜ C L (1 − β ) ∞ (cid:88) l =0 2 l (cid:88) k =0 l ( β − ( µ/ β )) ( k ) µ +2 β | u lk | ≥ ˜ C L (1 − β ) ∞ (cid:88) l =0 2 l (cid:88) k =0 (2 − l k ) µ | u lk | = ˜ C L (1 − β ) || u − P L u || L ( I,x µ/ ) . (C.8)Similarly, (C.8) after replacing ( u − P L u ) (cid:48) by ( u − P L u ) (cid:48)(cid:48) and ( u − P L u ) by ( u − P L u ) (cid:48) reads || ( u − P L u ) (cid:48)(cid:48) || L ( I,x µ/ β ) = ˜ C ∞ (cid:88) l = L +1 l l (cid:88) k =0 (2 − l k ) µ +4 β | u lk | ≥ ˜ C L (1 − β ) ∞ (cid:88) l =0 l l (cid:88) k =0 l ( β − ( µ/ β )) ( k ) µ +2 β | u lk | ≥ ˜ C L (1 − β ) ∞ (cid:88) l = L +1 l l (cid:88) k =0 (2 − l k ) µ +2 β | u lk | = ˜ C L (1 − β ) || ( u − P L u ) (cid:48) || L ( I,x µ/ β ) . (C.9)Finally, combining (C.8) and (C.9) yields the last estimate in (4.9). (cid:3) References [1] A. Alfonsi. On the discretization schemes for the CIR (and Bessel squared) processes.
Monte Carlo Methodsand Applications , (4), 1569-3961, 2005.[2] A. Alfonsi. A high-order discretization scheme for the CIR process: application to affine term structure andHeston models. Mathematics of Computation , (269), 209-237, 2010.[3] L. Andersen. Simple and efficient simulation of the Heston stochastic volatility model. Journal of ComputationalFinance , (3), 1-42, 2008.[4] J. Andreasen and B. Huge. ZABR – Expansions for the Masses. Preprint, SSRN//1980726, 2011.[5] J. Andreasen and B. Huge. Expanded forward volatility. Risk (January issue): 101-107, 2013.[6] A. Antonov and M. Spector. Advanced Analytics for the SABR model. Preprint, SSRN//2026350, 2012.[7] A. Antonov, M. Konikov and M. Spector. SABR Spreads its Wings.
Risk (August issue), 58-63, 2013.[8] A. Antonov, M. Konikov and M. Spector. The Free Boundary SABR: Natural Extension to Negative Rates.Preprint, SSRN//2557046, 2015.[9] A. Antonov, M. Konikov and M. Spector. Mixing SABR Models for Negative Rates. SSRN//2653682,2015.[10] P. Balland and Q. Tran. SABR Goes Normal.
Risk (June issue), 76-81, 2013.[11] C. Bayer, P. Friz and R. Loeffen. Semi-Closed form cubature and applications to financial diffusion models.
Quantitative Finance , (5): 769-782, 2013.[12] C. Bayer, J. Gatheral, M. Karlsmark. Fast Ninomiya-Victoir calibration of the Double-Mean-Reverting model. Quantitative Finance (11): 1813-1829, 2013.[13] S. Beuchler, R. Schneider, and C. Schwab. Multiresolution Weighted Norm Equivalences and Applications. Numer. Math. (98): 67-97, 2004.[14] N. Bouleau and L. Denis. Energy image density property and the lent particle method for Poisson measures.
Journal of Functional Analysis , (4), 1144-1174, 2009.[15] N. Bouleau and F. Hirsch. Dirichlet Forms and Analysis on Wiener Space. De Gruyter , 1991.[16] A. C. Cavalheiro. Weighted Sobolev Spaces and Degenerate Elliptic Equations.
Boletim da SociedadeParanaense de Matematica , v. 26 1-2 , 117132, 2008.[17] J.F. Chassagneux, A. Jacquier and I. Mihaylov. An explicit Euler scheme with strong rate of convergence fornon-Lipschitz SDEs. SIAM Joural on Financial Mathematics , (1), 2016.[18] B. Chen, C. W. Oosterlee and H. van der Weide. A low-bias simulation scheme for the SABR stochasticvolatility model. International Journal of Theoretical and Applied Finance , (2), 2012.[19] R. Cont and E. Voltchkova. A finite difference scheme for option pricing in exponential L´evy Models. SIAMJournal of Numerical Analysis , (4): 1596-1626, 2005.[20] R. Cont and E. Voltchkova. Integro-differential equations for option prices in exponential L´evy models. Financeand Stochastics , (3): 299-325, 2005.[21] S. Cox, M. Hutzenthaler and A. Jentzen. Local Lipschitz continuity in the initial value and strong completenessfor nonlinear stochastic differential equations. Preprint, arXiv:1309.5595, 2013. IRICHLET FORMS AND FINITE ELEMENT METHODS FOR THE SABR MODEL 29 [22] L. D¨oring, B. Horvath, J. Teichmann. Functional analytic (ir-)regularity properties of SABR-type processes.
International Journal of Theoretical and Applied Finance , (3), 2017.[23] W. Dahmen, A. Kunoth, and K. Urban. Biorthogonal spline wavelets on the interval - stability and momentconditions. Applied and Computational Harmonic Analysis , :132-196, 1999.[24] G. Deelstra and F. Delbaen. Convergence of discretized stochastic (interest rate) processes with stochastic driftterm. Applied Stochastic Models Data Analysis (1): 77-84, 1998.[25] R. Dautray and J.-L. Lions. Mathemathical Analysis and Numerical Methods for Science and Technology, Volume 6
Evolution Problems II . Springer, 2000.[26] S. De Marco, C. Hillairet, and A. Jacquier. Shapes of implied volatility with positive mass at zero.
SIAMJournal on Financial Mathematics , forthcoming.[27] S. Dereich, A. Neuenkirch, L. Szpruch. An Euler-type Method for the Strong approximation of the Cox-Ingersoll-Ross process.
Proceedings of the Royal Society A (2140): 1105-1115, 2012.[28] P. D¨orsek and J. Teichmann. A semigroup point of view on splitting schemes for stochastic (partial) differentialequations. Preprint, arXiv/1011.2651, 2010.[29] P. Doust. No-arbitrage SABR.
The Journal of Computational Finance , (3): 3-31, 2012.[30] S. N. Ethier and T. G. Kurtz. Markov Processes: Characterization and Convergence. Wiley, 1986.[31] E. Fabes, C. Kenig, R. Serapioni. The local regularity of solutions of degenerate elliptic equations. Communi-cations in Partial Differential Equations (1), 77-116, 1982.[32] M. Forde, A. Pogudin. The large-maturity smile for the SABR and CEV-Heston models. International Journalof Theoretical and Applied Finance , (8), 2013.[33] M. Forde, H. Zhang. Sharp tail estimates for the correlated SABR model. Preprint , 2014.[34] M. Fukushima, Y. Oshima and M. Takeda. Dirichlet Forms and Symmetric Markov Processes.
De Gruyter ,1994.[35] J. Garcia-Cuerva and J.L. Rubio de Francia. Weighted norm inequalities and related topics.
Elsevier, North-Holland Mathematics Studies , , Amsterdam, 1985.[36] K. Glau, Feynman-Kac-Darstellungen zur Optionspreisbewertung in L´evy-Modellen. PhD Thesis , Albert-Ludwigs-Universit¨at Freiburg, 2010.[37] A. Gulisashvili. Left-wing asymptotics of the implied volatility in the presence of atoms.
International Journalof Theoretical and Applied Finance , (2) 2015.[38] A. Gulisashvili, B. Horvath and A. Jacquier. On the probability of hitting the boundary of a Brownian motionon the SABR plane. Electronic Communications in Probability , (75): 1-13, 2016.[39] A. Gulisashvili, B. Horvath and A. Jacquier. Mass at zero in the uncorrelated SABR model and impliedvolatility asymptotics. Preprint, 2016.[40] I. Gy¨ongy and M. R´asonyi. A note on Euler approximations for SDEs with H¨older continuous diffusion coeffi-cients. Stochastic Process. Appl. (10): 2189-2200, 2011.[41] P. Hagan, D. Kumar, A. Lesniewski, and D. Woodward. Managing smile risk.
Wilmott Magazine (Septemberissue), 84-108, 2002.[42] P. Hagan, D. Kumar, A. Lesniewski, and D. Woodward. Arbitrage-free SABR.
Wilmott Magazine (Januaryissue): 60-75, 2014.[43] P. Hagan, A. Lesniewski, and D. Woodward. Probability distribution in the SABR model of stochastic volatility.
Large Deviations and Asymptotic Methods in Finance (Editors: P. Friz, J. Gatheral, A. Gulisashvili, A.Jacquier, J. Teichmann) , Springer Proceedings in Mathematics and Statistics, , 2015.[44] N. Hilber, O. Reichmann, C. Schwab, C. Winter. Computational methods for quantitative finance: Finiteelement methods for derivative pricing.
Springer Finance , 2013.[45] D. Hobson. Comparison results for stochastic volatility models via coupling.
Finance and Stochastics , :129-152, 2010.[46] B. Horvath. Robust methods for the SABR model and related processes. PhD Thesis , ETH Z¨urich, 2016.[47] M. Hutzenthaler, A. Jentzen and M. Noll, Strong convergence rates and temporal regularity for Cox-Ingersoll-Ross processes and Bessel processes with accessible boundaries. Preprint, arXiv:1403.6385, 2014.[48] O. Islah. Solving SABR in exact form and unifying it with LIBOR market model. Preprint, SSRN//1489428,2009.[49] B. Jourdain. Loss of martingality in asset price models with lognormal stochastic volatility.
Internat. Journalof Theoretical and Applied Finance , : 767-787, 2004.[50] O. Kallenberg. Foundations of Modern Probability. Springer series in statistics: Probability and its applica-tions , 2002.[51] A. Kufner. Weighted Sobolev Spaces. Teubner-Texte zur Mathematik,
Bd. 31 , Teubner Verlagsgesellschaft ,Leipzig, 1980.[52] A. Kufner and B. Opic. How to define reasonably weighted Sobolev spaces.
Commentationes MathematicaeUniversitatis Carolinae , (3): 537-554 , 1984.[53] F. Le Floc’h and G. J. Kennedy. Finite Difference Techniques for Arbitrage Free SABR. Preprint,ssrn.com/abstract=2402001, 2014.[54] J.-L. Lions and E. Magenes. Probl`emes aux limites non homog`enes et applications. Travaux et RecherchesMath´ematiquies , Dunod, Paris, 1968. [55] A.E. Lindsay and D.R. Brecher. Simulation of the CEV process and the local martingale property.
Mathematicsand Computers in Simulation , , 868-878, 2012.[56] P.-L. Lions and M. Musiela. Correlations and bounds for stochastic volatility models. Annales de l’InstitutHenri Poincar´e (C) Non Linear Analysis , (1): 1-16, 2007.[57] R. Lord. Fifty shades of SABR simulation. Presentation, Bachelier World Congress , /fiftyshadesbachelier.pdf,2014.[58] R. Lord, R. Koekkoek and D. v. Dijk. A comparison of biased simulation schemes for stochastic volatilitymodels.
Quantitative Finance (2): 177-194, 2010.[59] M. Matache, T. von Petersdorff, and C. Schwab. Fast deterministic pricing of options on L´evy driven assets. M2AN Math. Model. Numer. Anal. , (1): 37-71, 2004.[60] G. N. Milstein and J. Schoenmakers. Uniform approximation of the Cox-Ingersoll-Ross process. Preprint,arXiv:1312.0876v1, 2013.[61] B. Muckenhoupt. Weighted norm inequalities for the Hardy maximal functions. Transactions of the AmericanMathematical Society , , 207-226, 1972.[62] J. Ob(cid:32)l´oj. Fine-tune your Smile: Correction to Hagan et al. Wilmott Magazine (May issue), 2008.[63] L. Paulot. Asymptotic implied volatility at the second order with application to the SABR model.
LargeDeviations and Asymptotic Methods in Finance (Editors: P. Friz, J. Gatheral, A. Gulisashvili, A. Jacquier,J. Teichmann) , Springer Proceedings in Mathematics and Statistics,
Volume 110 , 2015.[64] T. von Petersdorff and C. Schwab. Wavelet discretizations of parabolic integrodifferential equations.
SIAM J.Numer. Anal. , (1): 159-180, 2003.[65] R. Rebonato. A simple approximation for the no-arbitrage drifts for LMM-SABR-family interest-rate models.Preprint, papers.ssrn.com/sol3/papers.cfm?astract-id=2560241, 2014.[66] R. Rebonato, K. McKay and R. White. The SABR/LIBOR Market Model: Pricing, Calibration and Hedgingfor Complex Interest-Rate Derivatives. Wiley , 2009.[67] M. Reed and B. Simon. Functional Analysis (Methods of modern mathematical physics), . AcademicPress, London, 1980.[68] O. Reichmann and C. Schwab. Wavelet solution of degenerate Kolmogoroff forward equations for exotic con-tracts in finance.
Recent Developments in Computational Finance: Foundations, Algorithms and Applications(Editors: Th. Gerstner and P. Kloeden) , Volume 14 , 2013.[69] O. Reichmann. Numerical option pricing beyond L´evy.
PhD Thesis , ETH Z¨urich, 2012.[70] M. R¨ockner and Z.-M. Ma. Introduction to the Theory of (Non-Symmetric) Dirichlet Forms.
Springer , 1991.[71] M. R¨ockner and N. Wielens. Dirichlet forms-closability and change of speed measure.
Infinite DimensionalAnalysis and Stochastic Processes ( Editor: S. Albeverio ), Research Notes in Mathematics. Pitman, BostonLondon Melbourne: 119-144, 1985.[72] A. Torchinsky. Real-Variable Methods in Harmonic Analysis.
Academic Press, San Diego, California , 1986.[73] B. O. Turesson. Nonlinear Potential Theory and Weighted Sobolev Spaces.
Springer Lecture Notes in Mathe-matics , , 2000.[74] D.Talay and L. Tubaro. Expansion of the global error for numerical schemes solving stochastic differentialequations. Stochastic Analysis and Applications , (4) 94-120, 1990.[75] P. Wilmott, S. Howison and J. Dewynne. The mathematics of financial derivatives: a student introduction.Cambridge University Press, Oxford, 1995. Department of Mathematics, Imperial College London
E-mail address : [email protected] European Investment Bank
E-mail address ::