A fourier interpolation method for numerical solution of FBSDEs: Global convergence, stability, and higher order discretizations

Oyono Ngou, Polynice; Hyndman, Cody Article A fourier interpolation method for numerical solution of FBSDEs: Global convergence, stability, and higher order discretizations Journal of Risk and Financial Management Provided in Cooperation with: MDPI – Multidisciplinary Digital Publishing Institute, Basel Suggested Citation: Oyono Ngou, Polynice; Hyndman, Cody (2022) : A fourier interpolation method for numerical solution of FBSDEs: Global convergence, stability, and higher order discretizations, Journal of Risk and Financial Management, ISSN 1911-8074, MDPI, Basel, Vol. 15, Iss. 9, pp. 1-32, https://doi.org/10.3390/jrfm15090388 This Version is available at: https://hdl.handle.net/10419/274909 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by/4.0/ Citation: Oyono Ngou, Polynice, and Cody Hyndman. 2022. A Fourier Interpolation Method for Numerical Solution of FBSDEs: Global Convergence, Stability, and Higher Order Discretizations. Journal of Risk and Financial Management 15: 388. https://doi.org/10.3390/ jrfm15090388 Academic Editor: Alexandru M. Badescu Received: 20 May 2022 Accepted: 30 August 2022 Published: 31 August 2022 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2022 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). Journal of Risk and Financial Management Article A Fourier Interpolation Method for Numerical Solution of FBSDEs: Global Convergence, Stability, and Higher Order Discretizations † Polynice Oyono Ngou and Cody Hyndman * Department of Mathematics and Statistics, Concordia University, 1455 Boulevard de Maisonneuve Ouest, Montr´ eal, QC H3G 1M8, Canada *Correspondence: [email protected] † A previous version of this paper was titled “Global convergence and stability of a convolution method for numerical solution of BSDEs” arXiv:1410.8595v1. Abstract: The convolution method for the numerical solution of forward-backward stochastic differential equations (FBSDEs) was originally formulated using Euler time discretizations and a uniform space grid. In this paper, we utilize a tree-like spatial discretization that approximates the BSDE on the tree, so that no spatial interpolation procedure is necessary. In addition to suppressing extrapolation error, leading to a globally convergent numerical solution for the FBSDE, we provide explicit convergence rates. On this alternative grid the conditional expectations involved in the time discretization of the BSDE are computed using Fourier analysis and the fast Fourier transform (FFT) algorithm. The method is then extended to higher-order time discretizations of FBSDEs. Numerical results demonstrating convergence are presented using a commodity price model, incorporating seasonality, and forward prices. Keywords: forward-backward stochastic differential equations; numerical solutions; fast Fourier transform Mathematics Subject Classification (2000): Primary: 60H10, 65C30; Secondary: 60H30 1. Introduction A variety of numerical methods for backward stochastic differential equations (BSDEs) and forward-backward stochastic differential equations (FBSDEs) have been developed recently. Different applications call for innovative techniques for the efficient resolution of these systems. In finance and economics BSDEs are used for option pricing and hedging El Karoui et al. (1997), reflected BSDEs are used for modeling American options El Karoui et al. (1997), and quadratic BSDEs play an important role in continuous-time recursive utility Duffie and Epstein (1992) and other utility maximization problems. Following the establishment of the well-posedness of BSDEs by Pardoux and Peng (1992) the first numerical procedures to emerge were partial differential equation (PDE) based methods such as the finite difference approach of Douglas et al. (1996). The PDE method is mainly devoted to coupled problems as is the spectral method of Ma et al. (2008). More recent numerical methods based on machine learning (Beck et al. 2019;Han and Long 2020;E et al. 2017,2019,2022) have been applied to the numerical solution of BSDEs which, given the connections between the theory of PDEs and their representations as BSDEs, in turn provides solutions to high dimensional PDEs. Spatial discretization methods were initiated by Chevance (1997) with a quantization approach to conditional expectations. However, it is only since Zhang (2004) and Bouchard and Touzi (2004) that a sound time discretization of (decoupled) FBSDEs is available. The quantization approach was then used in the multidimensional framework of Bally and Pages (2003); Bally et al. J. Risk Financial Manag. 2022,15, 388. https://doi.org/10.3390/jrfm15090388 https://www.mdpi.com/journal/jrfm J. Risk Financial Manag. 2022,15, 388 2 of 32 (2005) and the coupled FBSDE case of Delarue and Menozzi (2006). The theoretical basis for a multinomial approach for BSDEs was introduced by Briand et al. (2001) and Ma et al. (2002) and followed in practice by Peng and Xu (2011). Monte Carlo methods are the most prolific approach for numerical solutions of (F)BSDEs. They include the backward scheme of Zhang (2004), the Malliavin approach of Bouchard and Touzi (2004) and Crisan et al. (2010), the least-square regression approach of Gobet et al. (2005) and the iterative schemes Bender and Denk (2007) and Bender and Zhang (2008). Additional approaches to numerical solution of BSDEs include the cubature method Crisan and Manolarakis (2012,2014), Fourier-cosine expansions Huijskens et al. (2016); Ruijter and Oosterlee (2015,2016), and the convolution method Hyndman and Oyono Ngou (2017). In introducing the convolution method, Hyndman and Oyono Ngou (2017) developed a local discretization error which includes an extrapolation error. The extrapolation error component is exclusively produced by the fast Fourier transform (FFT) algorithm and the underlying trigonometric interpolation used to compute conditional expectations. To improve the performance of the convolution method it is desirable to eliminate the extrapolation error, and improve the error bound, with an alternative implementation of the FFT algorithm. In this paper, we propose an alternative space grid for the convolution method, instead of the rectangular grid in Hyndman and Oyono Ngou (2017), which eliminates the extrapolation error, leads to a globally convergent numerical solution for the (F)BSDE, and provides explicit convergence rates. We also apply the numerical method to the Runge– Kutta schemes for FBSDEs proposed by Chassagneux and Crisan (2014). The tree-like nature of the alternative grid avoids extrapolations and leads to a global error bound for the BSDE approximate solutions. Further, the implementation of the convolution method originally presented in Hyndman and Oyono Ngou (2017) is simplified by using an alternative parametric transformation to enforce the necessary periodic boundary conditions. The paper is organized as follows. Section 2reviews the explicit Euler time discretization of BSDEs, recalls the convolution method of Hyndman and Oyono Ngou (2017) for computing the necessary conditional expectations, gives a description of an alternative spatial discretization, and provides a generic implementation of the convolution method on this grid using the discrete Fourier transform. The section ends with a global error analysis. Section 3extends the Fourier interpolation method to higher order time discretizations of FBSDEs and includes the related global error analysis. Finally, Section 4presents a numerical implementation, in the context of a commodity price model, that illustrates the theoretical results. Section 5concludes. 2. The Fourier Interpolation Method In this section, we introduce an alternative grid that removes extrapolation error of Hyndman and Oyono Ngou (2017), after a quick presentation of the Euler scheme for FBSDEs. Section 2.4 presents a global error analysis of the Fourier interpolation method on this alternative grid and under the Euler scheme. 2.1. Time Discretization Let Ω,P,F,{Ft}t∈[0,T] be a complete filtered probability space generated by a d−dimensional Wiener process W. We seek a numerical solution to the FBSDE      dXt=a(t,Xt)dt +σ(t,Xt)dWt −dYt=f(t,Xt,Yt,Zt)dt −Z∗ tdWt X0=x0,YT=ξ (1) The forward drift a:[ 0, T]×Rd→Rd , the forward volatility σ:[ 0, T]×Rd→Rd×d , the driver f:[ 0, T]×Rd×R×Rd→R are deterministic functions. The initial condition is x0∈Rd and the terminal condition takes the Markovian form ξ=g(XT) where g:Rd→R . J. Risk Financial Manag. 2022,15, 388 3 of 32 The FBSDE coefficients satisfy Assumption 1so that existence and uniqueness of the FBSDE solution (X,Y,Z)is assured. Assumption 1. There exist positive constants K1 , K2K3 , and K4 such that the coefficients of the FBSDE (1) satisfy |a(t,x1)−a(t,x2)|+kσ(t,x1)−σ(t,x2)k2≤K1|x1−x2|(2) |a(t,x)|+kσ(t,x)k2≤K2(3) |f(t,x1,y,z)−f(t,x2,y,z)|≤K1|x1−x2|(4) |f(t,x,y1,z1)−f(t,x,y2,z2)|≤K1(|y1−y2|+|z1−z2|)(5) |f(t,x,y,z)|≤K3(1+|x|+|y|+|z|)(6) for any t ∈[0, T], x,x1,x2∈Rd, y,y1,y2∈R, z,z1,z2∈Rd. Moreover σ2:=σσ∗is (uniformly) invertible, continuous and bounded k(σ2(t,x))−1k2≤K4(7) for any t ∈[0, T], x ∈Rd. In addition, the terminal value is square integrable kξk2 L2:=Eh|g(XT)|2i<∞. (8) Remark 1. Assumption 1makes no explicit assumption on g except square integrability. More restrictive assumptions on g , namely that g is twice continuously differentiable, shall be specified in the main results of Sections 2and 3which provide explicit rates of convergence for the Fourier interpolation algorithms. However, similar to Crisan and Manolarakis (2012), see also Turkedjiev (2015), it should be possible to consider g non-differentiable using a mollification argument on the terminal condition. 1 However, we shall not follow this approach in this paper so as to keep the focus on the main contributions which are the overall approach, implementation of the convolution method on the tree-like grid, developing the Runge–Kutta discretization schemes, and explicit convergence rates. The time discretization of the FBSDE (1) on the time partition π={t0= 0 <t1< . . . <tn=T}consists of the explicit Euler scheme given by                  Xπ 0=x0 Xπ ti+1=Xπ ti+a(ti,Xπ ti)∆i+σ(ti,Xπ ti)∆Wi Zπ tn=0, Yπ tn=ξπ Zπ ti=1 ∆iEhYπ ti+1∆Wi|Ftii Yπ ti=EhYπ ti+1|Ftii+f(ti,Xπ ti,EhYπ ti+1|Ftii,Zπ ti)∆i (9) with ∆i=ti+1−ti and ∆Wi=Wti+1−Wti . Under an additional Lipschitz condition on the function g we have, from Zhang (2004) and Bouchard and Touzi (2004), that the quadratic discretization error E2 π:=max 0≤i<nE"sup t∈[ti,ti+1]Yt−Yπ ti 2#+ n−1 ∑ i=0 EZti+1 tiZs−Zπ ti 2ds(10) is of first-order in time, i.e., E2 π=O(|π|)(11) where |π|=max i ∆i. (12) J. Risk Financial Manag. 2022,15, 388 4 of 32 Following Hyndman and Oyono Ngou (2017) and Oyono Ngou (2014), the approximate solution ui and the approximate gradient ˙ ui at time node ti , i= 0, 1, . . . , n− 1, are given by ui(x) = ˜ ui(x) + ∆if(ti,x,˜ ui(x),σ(ti,x)˙ ui(x)) (13) σ(ti,x)˙ ui(x) = 1 ∆i EhYπ ti+1σ(ti,Xπ ti)∆Wi|Xπ ti=xi =1 ∆iZRd(y−∆ia(ti,x))ui+1(x+y)hi(y|x)dy (14) where the intermediate solution ˜ uiis given by ˜ ui(x) = EhYπ ti+1|Xπ ti=xi =ZRdui+1(y+x)hi(y|x)dy, (15) un=gand hiis a Gaussian density hi(y|x) = (2π)−d 2k∆iσ2(ti,x)k−1 2 2exp−1 2∆i y∗(σ2(ti,x))−1y, (16) and where y=y−∆ia(ti,x)with characteristic function φi(ν,x) = expi∆iν∗a(ti,x)−1 2∆iν∗σ2(ti,x)ν. (17) Let F and F−1 denote the Fourier transform operator and the inverse Fourier transform operator, respectively, F[θ](ν) = ZRde−ix∗νθ(x)dx (18) F−1[θ](x) = 1 (2π)dZRdeiν∗xθ(ν)dν. (19) The density function hiand the characteristic function φisatisfy the relations hi(y|x) = 1 (2π)dZRde−iν∗yφi(ν,x)dν(20) yhi(y|x) = −1 (2π)dZRde−iν∗yi∇νφi(ν,x)dν =∆ia(ti,x)hi(y|x) + ∆iσ∗(ti,x) (2π)dZRde−iν∗yiνφi(ν,x)dν(21) Using the relationships between the characteristic function and the density function (20) and (21) leads to the representation ˜ ui(x) = F−1[F[ui+1](ν)φi(ν,x)](x)(22) ˙ ui(x) = σ∗(ti,x)F−1[F[ui+1](ν)iνφi(ν,x)](x)(23) for Equations (14) and (15) under integrability condition on the approximate solution ui+1 . In the sequel, we restrict the analysis to the one-dimensional case with d=1. 2.2. Space Discretization The space discretization is performed on a tree-like grid using three parameters: the increment length l> 0, the even number N∈N∗ of space steps on the increment length, J. Risk Financial Manag. 2022,15, 388 5 of 32 and the initial number N0∈N of increment intervals. Hence, the space step is constant and uniform on the grid ∆x=l N. (24) At the time node ti , i= 0, 1, . . . , n , the space domain is restricted on an interval of length Nilcentred at x0and discretized uniformly with NiNspace steps where Ni=N0+i(25) giving the nodes xik =x0−Nil 2+k∆x,k=0, 1, . . . , NiN. (26) In particular, the relation xik =xi+1,k+N 2,k=0, 1, . . . , NiN. (27) holds since the restricted interval at each time node is obtained by evenly increasing the previous one with an interval of length l . If N0= 0 then the space grid at mesh time t0 is composed of the single point x00 =x0. (28) Figure 1gives examples of alternative grids. P. Oyono Ngou & C. Hyndman Fourier interpolation method for FBSDEs August 26, 2022 Figure 1: Examples of alternative grids (a) N0= 0 and N= 2. b b b b b b b b b t0t1t2 l x0,0=x0 x1,0 x1,2 x2,0 x2,4 (b) N0= 1 and N= 2 b b b b b b b b b b b b b b b t0t1t2 l x0,0 x0,1=x0 x0,2 x1,0 x1,4 x2,0 x2,6 2.3 Implementation In order to compute numerical approximations of equations (2.22) and (2.23) at time node ti,i= 0,1,...,n−1, we introduce the generic functions θi:R→R,ψ:R2→Cand θi+1 :R→Rsuch that θi(x) = F−1[F[θi+1](ν)ψ(ν, x)] (x).(2.31) We assume that the function θi+1 satisfies the periodicity boundary value equalities of Assumption 2.2. Assumption 2.2. The generic function θi+1 satisfies θi+1 (xi+1,0) = θi+1 xi+1,NNi+1 (2.32) ∂θi+1 ∂x (xi+1,0) = ∂θi+1 ∂x xi+1,NNi+1 .(2.33) Hence, the Fourier integral F[θi+1](ν) = Z∞ −∞ e−iνxθi+1(x)dx (2.34) is restricted on the interval [x0−Ni+1l 2, x0+Ni+1l 2] = [xi+1,0, xi+1,NNi+1 ]and discretized using the grid points {xi+1,k}Ni+1N k=0 with a quadrature rule with weights {wk}Ni+1N k=0 . As to the inverse Fourier integral of equation (2.31) we restrict it on the interval [−L 2,L 2]and discretize it with lower Riemann sums at the Fourier space grid point {νi+1,k}Ni+1N k=0 . Let Dand D−1denote the discrete Fourier transform and the inverse discrete Fourier transform respectively D[{x}N−1 i=0 ]k=1 N N−1 X j=0 e−ijk 2π Nxj(2.35) D−1[{x}N−1 i=0 ]k= N−1 X j=0 eijk 2π Nxj.(2.36) Then the discretization procedure leads to the approximation θi(xi+1,k)≈(−1)kD−1h{ψ(νi+1,j , xi+1,k)D[θi+1]j}Ni+1N−1 j=0 ik (2.37) where D[θi+1]j=Dh{(−1)s˜wsθi+1(xi+1,s)}Ni+1N−1 s=0 ij (2.38) 5 Figure 1. Examples of alternative grids. The convolution relations of Equations (22) and (23) call for a discretization of the Fourier space as well. At each mesh time ti , i= 1, 2, . . . , n , the Fourier space is restricted on an interval of length L centred at zero ( 0 ) and discretized with NiN space steps. The equidistant nodes are thus of the form νik =−L 2+k∆νi,k=0, 1, . . . , NiN(29) where ∆νi=L NiN. The Nyquist relation2holds whenever Lis such that Ll =2πN. (30) J. Risk Financial Manag. 2022,15, 388 6 of 32 2.3. Implementation In order to compute numerical approximations of Equations (22) and (23) at time node ti , i= 0, 1, . . . , n− 1, we introduce the generic functions θi:R→R , ψ:R2→C and θi+1:R→Rsuch that θi(x) = F−1[F[θi+1](ν)ψ(ν,x)](x). (31) We assume that the function θi+1 satisfies the periodicity boundary value equalities of Assumption 2. Assumption 2. The generic function θi+1satisfies θi+1(xi+1,0)=θi+1xi+1,NNi+1(32) ∂θi+1 ∂x(xi+1,0)=∂θi+1 ∂xxi+1,NNi+1. (33) Hence, the Fourier integral F[θi+1](ν) = Z∞ −∞e−iνxθi+1(x)dx (34) is restricted on the interval [x0−Ni+1l 2 , x0+Ni+1l 2]=[xi+1,0 , xi+1,NNi+1] and discretized using the grid points {xi+1,k}Ni+1N k=0 with a quadrature rule with weights {wk}Ni+1N k=0 . As to the inverse Fourier integral of Equation (31) we restrict it on the interval [−L 2 , L 2] and discretize it with lower Riemann sums at the Fourier space grid point {νi+1,k}Ni+1N k=0. Let D and D−1 denote the discrete Fourier transform and the inverse discrete Fourier transform, respectively D[{x}N−1 i=0]k=1 N N−1 ∑ j=0 e−ijk 2π Nxj(35) D−1[{x}N−1 i=0]k= N−1 ∑ j=0 eijk 2π Nxj. (36) Then the discretization procedure leads to the approximation θi(xi+1,k)≈(−1)kD−1hψ(νi+1,j,xi+1,k)D[θi+1]jNi+1N−1 j=0ik(37) where D[θi+1]j=Dh{(−1)s˜ wsθi+1(xi+1,s)}Ni+1N−1 s=0ij(38) and the weights {˜ w}Ni+1N−1 k=0are given by ˜ wk=wk+wNi+1Nδk,0. (39) where δstands for the Kronecker delta. The relation of Equation (27) allows us to write θi(xik)≈(−1)k+N 2D−1hψ(νi+1,j,xik)D[θi+1]jNi+1N−1 j=0ik+N 2 (40) for k=0, 1, . . . , NiN. In Equation (40), the generic function ψ depends on the space node xik . If the relation generalizes for all space nodes xik , k= 0, 1, . . . , NiN , the function values θi(xik) , k= 0, 1, . . . , NiN , can not be computed with a single direct FFT procedure. Instead, a separate FFT procedure using the values of the generic function ψ at xik is needed to compute the function value θi(xik) . Nonetheless, the vector-matrix representation of the J. Risk Financial Manag. 2022,15, 388 7 of 32 FFT procedure in Equation (40) allows the computation of all function values θi(xik) with a matrix multiplication. In the vector-matrix representation, Equation (40) is θi(xik) = (−1)k+N 2ˆ Fk+N 2 Ψ(xik)D[θi+1] where ˆ Fk+N 2 is the (k+N 2) th row of the Ni+1N dimension inverse FFT matrix ˆ F and Ψ(xik) is the Ni+1N dimension diagonal matrix built with the values {ψ(νi+1,j , xik)}Ni+1N−1 j=0 . Let Θ(i)be the NiNdimension vector of the function values θi(xik)such that Θ(i) 1+k=θi(xik)(41) for k=0, 1, . . . , NiN. The matrix representation gives Θ(i)=ˆ Ψ(i)D[θi+1](42) where ˆ Ψ(i)is the (NiN+1)×Ni+1Nmatrix such that ˆ Ψ(i) 1+k,1+j= (−1)k+N 2¯ ωj(k+N 2) iψ(νi+1,j,xik)(43) with ¯ ωi=ei2π(Ni+1N)−1,k=0, 1, . . . , NiNand j=0, 1, . . . , Ni+1N−1. The requirements of Assumption 2can easily be satisfied. Given a function η:[a , b]→ Rand η∈ C1, if we consider the transformation ηα,β(x) = η(x) + αx2+βx(44) then the parameters α and β can be chosen such that the transform function and its derivative have equal values at the boundaries of any interval. The following lemma gives a method to select the coefficients α and β for the transform of Equation (44) such that Assumption 2holds in general. Lemma 1. Suppose the real function η∈ C1[a , b] is differentiable and let ηα,β be its transformed function as defined in Equation (44). Then α= ∂η ∂x(a)−∂η ∂x(b) 2(b−a),(45) β=η(a)−η(b) (b−a)−α(b+a)(46) solve the system of linear equations defined by (ηα,β(a) = ηα,β(b) ∂ηα,β ∂x(a) = ∂ηα,β ∂x(b).(47) Proof. The second equation of the system (47) gives Equation (45) in a straightforward manner. Equation (46) is given by the first equation of the system. Hence, the numerical discretization may be applied on the transformation uα,β i+1 at time node ti but a correction must be performed so to recover the values of the intermediate solution ˜ ui and the approximate gradient ˙ ui . The next theorem gives the representation under the transform of Equation (44). J. Risk Financial Manag. 2022,15, 388 8 of 32 Theorem 1. Let uα,β i+1 be the alternative transform defined in Equation (44) of the approximate solution ui+1 . Then the intermediate solution ˜ ui and the approximate gradient ˙ ui in Equations (22) and (23) satisfy ˜ ui(x) = F−1[F[uα,β i+1](ν)φ(ν,x)](x)−α[(x+∆ia(ti,x))2+∆iσ2(ti,x)] −β(x+∆ia(ti,x)) (48) ˙ ui(x) = σ(ti,x)F−1[F[uα,β i+1](ν)iνφ(ν,x)](x)−σ(ti,x)[2α(x+∆ia(ti,x)) + β]. (49) Proof. The proof follows the steps of Theorem 3.1 of Hyndman and Oyono Ngou (2017) using the transformation introduced in Equation (44) and the relations (20) and (21) between the density function hiand the characteristic function φi. Algorithm 1details the numerical procedure on the space grid and produces numerical solutions {uik}NiN k=0 , {˜ uik}NiN k=0 and {˙ uik}NiN k=0 for the approximate solution ui , the intermediate solution ˜ ui and the approximate gradient ˙ ui , respectively, i= 0, 1, . . . , n− 1. We next consider error estimates under the alternative discretization. Algorithm 1 Fourier Interpolation Method on Alternative Grid 1. Discretize the restricted real space [x0−Nnl 2 , x0+Nnl 2] and the restricted Fourier space [−L 2 , L 2] with NnN space steps so to have the real space nodes {xnk}NnN k=0 and {νnk}NnN k=0 2. Value un(xnk) = g(xnk) 3. For any ifrom n−1 to 0 (a) Compute αand βdefining the transform of Equation (44), such that θi+1=uα,β i+1(50) and θi+1satisfies the boundary conditions of Equations (32) and (33). (b) Compute θi(xik)through Equation (40) for k=0, 1, . . . , NiNwith ψ(ν,x) = φi(ν,x)(51) and retrieve the values ˜ uik as ˜ uik =θi(xik)−α[(xik +∆ia(ti,xik))2+∆iσ2(ti,xik)] −β(xik +∆ia(ti,xik)). (52) (c) Compute θi(xik)through Equation (40) for k=0, 1, . . . , NiNwith ψ(ν,x) = iνσ(ν,x)φi(ν,x)(53) and retrieve the values ˙ uik as ˙ uik =θi(xik)−σ(ti,xik)[2α(xik +∆ia(ti,xik)) + β]. (54) (d) Compute the values uik as uik =˜ uik +∆if(ti,xik,˜ uik,˙ uik)(55) for k=0, 1, . . . , NiNthrough Equation (13). (e) Update the real space grid with Equation (27) and the Fourier space grid by discretizing the interval [−L 2 , L 2] with NiN space steps so to have the real space nodes {xik}NiN k=0and {νik}NiN k=0. J. Risk Financial Manag. 2022,15, 388 15 of 32 0 0 0 0 0 0 γ2γ20 0 γ20 1 1 −1 2γ2 1 2γ20β11−β1 3.2. Further Simplification From the q -stage Runge–Kutta scheme for BSDEs, one notices that we have at least 2 q conditional expectations to compute at each time step. These conditional expectations can be simplified and made more suitable for numerical implementation if we consider a reasonable time discretization of the forward SDE. Hence, we make the following assumption. Assumption 3 (Forward process discretization) . The following are assumed throughout this section. 1. The forward SDE is discretized with the piecewise constant process Xπ such that for t∈[ti,ti+1)we have Xπ t=Xπ tipathwise. 2. The forward SDE time discretization with global error EX,πis of order m >0, i.e., E2 X,π:=max 0≤i≤nkXti−Xπ tik2 L2=O(|π|2m). (87) 3. The forward SDE time discretization admits the conditional characteristic functions φi: Rd×Rd→C φi(ν,x) = Eeiν∗Xπ ti+1−Xπ ti|Xπ ti=x(88) and Φi,j:Rd×Rd→Cd Φi,j(ν,x) = EHϕj ti,j,γj∆ieiν∗Xπ ti+1−Xπ ti|Xπ ti=x(89) for 0≤i<n and 1<j≤q+1with ϕq+1=ϕ1. 4. There are positive constants p0, q0, s0,K0and C0>0such that max inf s∈R+ d e−s∗tφi(is,x), inf s∈R+ d e−s∗tφi(−is,x)!≤e−K0∆−s0 i|t|q0, (90) ∀t∈R+ d , hence the discrete version of the forward process has conditional exponential moments. In addition, ZRd|φi(ν,x)|dν+max 1<j≤q+1ZRdΦi,j(ν,x)dν≤C0∆−p0 i. (91) It ˆ o-Taylor expansion based schemes are an example of SDE discretization satisfying the conditions of Assumption 3. A more complete presentation of these schemes can be found in Kloeden and Platen (1992). The next theorem gives a simplification of the BSDE time discretization expressions. Theorem 4. Under Assumption 3(1), the solution of the q-stage Runge–Kutta scheme satisfies {(Yπ i,j,Zπ i,j)}q+1 j=2∈ Fti(92) J. Risk Financial Manag. 2022,15, 388 16 of 32 for 0≤i<n. Consequently, we can write Zπ i,j=EtihHϕj ti,j,γj∆iYπ ti+1+∆iβj,1 f(ti+1,Xπ ti+1,Yπ ti+1,Zπ ti+1)i (93) Yπ i,j=EtihYπ ti+1+∆iαj,1 f(ti+1,Xπ ti+1,Yπ ti+1,Zπ ti+1)i +∆i j ∑ k=2 αjk f(ti,k,Xπ ti,Yπ i,k,Zπ i,k)(94) for 0≤i<n and 1<j≤q+1where ϕq+1=ϕ1,βq+1,1 =β1and αq+1,k=αk. Proof. Clearly (Yπ i,q+1 , Zπ i,q+1) = (Yπ ti , Zπ ti)∈ Fti from Equations (75) and (76). For 1 <j≤q and 0 ≤i<n, we have Yπ i,j=E"Yπ ti+1+∆i j ∑ k=1 αjk f(ti,k,Xπ ti,k,Yπ i,k,Zπ i,k)|Xπ ti,j#(starting from Equation (78)) =E"Yπ ti+1+∆i j ∑ k=1 αjk f(ti,k,Xπ ti,k,Yπ i,k,Zπ i,k)|Xπ ti# (by Assumption 3since ti,j∈[ti,ti+1)) =Eti"Yπ ti+1+∆i j ∑ k=1 αjk f(ti,k,Xπ ti,k,Yπ i,k,Zπ i,k)# so that Yπ i,j∈ Fti. Similar arguments also show that Zπ i,j∈ Ftistarting from Equation (77). Since {(Yπ i,j , Zπ i,j)}q+1 j=2∈ Fti , we naturally get Equation (94) from Equations (78) and (76). In addition, knowing that EtihHϕj ti,j,(γi−γk)∆ii=0,1<k<j(95) leads to Equation (93) from Equations (77) and (75). As a consequence of Assumption 3, if the q− stage Runge–Kutta scheme and the forward SDE time discretization are of order m> 0 then error of the FBSDE numerical solution defined as EX,π+Eπ is of order m . We must hence choose the Runge–Kutta scheme and the SDE scheme accordingly. 3.3. Fourier Representation Following Theorem 4, the intermediate solutions {(ui,j , ˙ ui,j)}q+1 j=2 at mesh time ti , 0≤i<n, are given by ˙ ui,j(x) = EhHϕj ti,j,γj∆i˜ ui+1(Xπ ti+1,βj,1)|Xπ ti=xi(96) ui,j(x) = Eh˜ ui+1(Xπ ti+1,αj,1)|Xπ ti=xi+∆i j ∑ k=2 αjk f(ti,k,x,ui,k(x),˙ ui,k(x)) (97) for 1 <j≤q+ 1 with ϕq+1=ϕ1 , βq+1,1 =β1 and αq+1,k=αk . The approximate solution uiand approximate gradient ˙ uiat mesh time ti, 0 ≤i<n, are then ui(x) = ui,q+1(x)(98) ˙ ui(x) = ˙ ui,q+1(x)(99) with ˜ ui+1(x,α) = ui+1(x) + ∆iαf(ti+1,x,ui+1(x),˙ ui+1(x)) (100) J. Risk Financial Manag. 2022,15, 388 17 of 32 and un(x) = g(x)(101) ˙ un(x) = σ∗(T,x)∇g(x). (102) In this setting, we have that ui,j(x) = Eh˜ ui+1(Xπ ti+1,αj,1)|Xπ ti=xi+∆i j ∑ k=2 αjk f(ti,k,x,ui,k(x),˙ ui,k(x)) (103) Note that Eh˜ ui+1(Xπ ti+1,αj,1)|Xπ ti=xi=Ex ti1 (2π)dZRdeiν∗Xπ ti+1F˜ ui+1(., αj,1)(ν)dν =1 (2π)dZRdEx tiheiν∗Xπ ti+1iF˜ ui+1(., αj,1)(ν)dν(using Fubini’s theorem) =1 (2π)dZRdeiν∗xφi(ν,x)F˜ ui+1(., αj,1)(ν)dν. (104) Therefore, by (103) and (104), we have ui,j(x) = F−1F˜ ui+1(., αj,1)(ν)φi(ν,x)(x) +∆i j ∑ k=2 αjk f(ti,k,x,ui,k(x),˙ ui,k(x)) (105) whenever ˜ ui+1(., α)is Lebesgue integrable. As to the intermediate solutions ˙ ui,j, 0 ≤i<nand 1 <j≤q+1, we have ˙ ui,j(x) = EhHϕj ti,j,γj∆i˜ ui+1(Xπ ti+1,βj,1)|Xπ ti=xi =Ex tiHϕj ti,j,γj∆i 1 (2π)dZRdeiν∗Xπ ti+1F˜ ui+1(., βj,1)(ν)dν =1 (2π)dZRdEx tihHϕj ti,j,γj∆ieiν∗Xπ ti+1iF˜ ui+1(., βj,1)(ν)dν(using Fubini’s theorem) =1 (2π)dZRdeiν∗xΦi,j(ν,x)F˜ ui+1(., βj,1)(ν)dν =F−1F˜ ui+1(., αj,1)(ν)Φi,j(ν,x)(x)(106) for an integrable function ˜ ui+1(., α). Even if the expressions in Equations (105) and (106) appear too general, they are implementable with the Fourier interpolation method for d= 1 in various particular cases. One can retrieve the characteristics φi and Φi,j and also perform the corrections due to the transform of Equation (44) for many SDE time discretizations. The following lemma helps in retrieving the conditional characteristics. Lemma 2. The conditional characteristics Φi,jare Φi,j(ν,x) = Ex tiH∗ i,jiνeiν∗Xπ ti+1−Xπ ti(107) with Hi,j=1 γj∆iZti+1 ti,j DsXπ ti+1ϕj s−ti,j γj∆i!ds (108) where DsXπ ti+1is the Malliavin derivative of Xπ ti+1given Xπ ti=x. J. Risk Financial Manag. 2022,15, 388 18 of 32 Proof. The lemma is proved by applying the duality formula and the chain rule successively to Equation (89). 3.3.1. Half-Order Itˆ o-Taylor Schemes The Euler scheme constitutes the main example of half-order It ˆ o-Taylor scheme with step Xπ ti+1=Xπ ti+a(ti,Xπ ti)∆i+σ(ti,Xπ ti)∆Wi. In addition, we have that Ds∆i=0d×1 and Ds∆Wi=Id×d for s∈(ti , ti+1) where 0 and I are the zero matrix and the identity matrix, respectively. Hence, DsXπ ti+1=σ(ti,Xπ ti) so we get, from Equation (107), that Φi,j(ν,x) = σ∗(ti,x)iνEx tieiν∗Xπ ti+1−Xπ ti(since ϕj∈ B0), =σ∗(ti,x)iνφi(ν,x). (109) The conditional characteristic function is explicitly given by φi(ν,x) = exp∆iiν∗a(ti,x)−1 2ν∗σ2(ti,x)ν (110) since the increment has a Gaussian distribution. Equations (105) and (106) along with the characteristics of Equations (110) and (109) define the Fourier method under half-order It ˆ o-Taylor schemes and the method is implementable in one dimension ( d= 1) with the procedure given in Section 2.3. The following theorem generalizes the result of Theorem 1to Runge–Kutta schemes under half-order Itˆ o-Taylor schemes. Theorem 5. Let ˜ uα,β i+1( ., y) be the alternative transform defined in Equation (44) of the approximate solution ˜ ui+1( ., y) . Then the intermediate solutions ui,j and ˙ ui,j in Equations (105) and (106) satisfy ui,j(x) = F−1[F[˜ uα,β i+1(., αj,1)](ν)φi(ν,x)](x) −α[(x+∆ia(ti,x))2+∆iσ2(ti,x)] −β(x+∆ia(ti,x)) +∆i j ∑ k=2 αjk f(ti,k,x,ui,k(x),˙ ui,k(x)) (111) ˙ ui,j(x) = σ(ti,x)F−1[F[˜ uα,β i+1(., βj,1)](ν)iνφi(ν,x)](x) −σ(ti,x)[2α(x+∆ia(ti,x)) + β]. (112) under a half-order Itˆo-Taylor scheme when d =1. 3.3.2. First-Order Itˆ o-Taylor Schemes Consider the first-order scheme Xπ ti+1=Xπ ti+a(ti,Xπ ti)∆i+σ(ti,Xπ ti)∆Wi+σ2(ti,Xπ ti)Zti+1 tiZt ti dWudWt. Then knowing that DsZti+1 tiZt ti dWudWt=D(∆Wi), (113) J. Risk Financial Manag. 2022,15, 388 19 of 32 using the fundamental theorem of calculus where D(x) is the diagonal matrix composed with the elements of x , for s∈(ti , ti+1) , the Malliavin derivative of the discretized forward process is given by DsXπ ti+1=σ(ti,x) + σ2(ti,x)D(∆Wi). (114) Equation (107) leads to Φi,j(ν,x) = σ∗(ti,x)iνEx tieiν∗Xπ ti+1−Xπ ti +Ex tiD(∆Wi)σ2(ti,x)iνeiν∗Xπ ti+1−Xπ ti(since ϕj∈ B0) =σ∗(ti,x)iνφi(ν,x) + Ex tiD(σ2(ti,x)iν)∆Wieiν∗Xπ ti+1−Xπ ti = (I−D(σ2(ti,x)iν)∆i)−1σ∗(ti,x)iνφi(ν,x)(115) since Ex ti∆Wieiν∗Xπ ti+1−Xπ ti=σ∗(ti,x)∆iEx tiiνeiν∗Xπ ti+1−Xπ ti +D(σ2(ti,x)iν)∆iEx ti∆Wieiν∗Xπ ti+1−Xπ ti using the duality formula, so that Ex ti∆Wieiν∗Xπ ti+1−Xπ ti=ζi(ν,x)∆iσ∗(ti,x)Ex tiiνeiν∗Xπ ti+1−Xπ ti =ζi(ν,x)∆iσ∗(ti,x)iνφi(ν,x) with ζi(ν,x) = (I−D(σ2(ti,x)iν)∆i)−1. (116) As to the conditional characteristic φi, it can be easily derived as φi(ν,x) = det(ζi(ν,x))1 2exp1 2iν∗ζ−1 i(ν,x)1d×1+iν∗κi(x)(117) where κi(x) = a(ti,x)∆i−1 2(σ2(ti,x)∆i+1)1d×1 knowing that Xπ ti+1−Xπ ti , given Xπ ti=x , is an affine function of a multivariate non-central χ2random variable with 1 degree of freedom and non-centrality parameters 1. Equations (105) and (106) along with the expressions in Equations (115) and (117) characterize the method under first order discretizations on the forward process. The procedure introduced in Section 2.3 allows us to do the necessary computations given the characteristics φiand Φi,jand using the following theorem. Theorem 6. Let ˜ uα,β i+1( ., y) be the alternative transform defined in Equation (44) of the approximate solution ˜ ui+1( ., y) . Then the intermediate solutions ui,j and ˙ ui,j in Equations (105) and (106) satisfy J. Risk Financial Manag. 2022,15, 388 20 of 32 ui,j(x) = F−1[F[˜ uα,β i+1(., αj,1)](ν)φi(ν,x)](x) −α(x+∆ia(ti,x))2+∆iσ2(ti,x) + 1 2∆2 iσ4(ti,x) −β(x+∆ia(ti,x)) + ∆i j ∑ k=2 αjk f(ui,k(x),˙ ui,k(x)) (118) ˙ ui,j(x) = σ(ti,x)F−1[F[˜ uα,β i+1(., βj,1)](ν)iνζi(ν,x)φi(ν,x)](x) −σ(ti,x)h2αx+∆ia(ti,x) + ∆iσ2(ti,x)+βi(119) under a first-order Itˆo-Taylor scheme when d =1. Proof. By the definition of the alternative transform, we must have that ui,j(x) = F−1[F[˜ uα,β i+1(., αj,1)](ν)φi(ν,x)](x)−Ex tihα(Xπ ti+1)2+βXπ ti+1i +∆i1{j>2} j−1 ∑ k=2 αjk f(ui,k(x),˙ ui,k(x)). (120) Notice that Ex tihXπ ti+1i=x+∆ia(ti,x)(121) and Ex tih(Xπ ti+1)2i=Ex tihXπ ti+1i2+Varx ti[Xπ ti+1] = (x+∆ia(ti,x))2+Ex ti"σ(ti,x)∆Wi+σ2(ti,x)Zti+1 tiZt ti dWudWt2# = (x+∆iσ(ti,x))2+∆iσ2(ti,x) + 1 2∆2 iσ4(ti,x). (122) Equations (120)–(122) lead to the expression for ui,jin Equation (118). The definition of the alternative transform also requires ˙ ui,j(x) = σ(ti,x)F−1[F[˜ uα,β i+1(., βj,1)](ν)iνζi(ν,x)φi(ν,x)](x) −Ex tihHϕj ti,j,γj∆i(α(Xπ ti+1)2+βXπ ti+1)i =σ(ti,x)F−1[F[˜ uα,β i+1(., βj,1)](ν)iνζi(ν,x)φi(ν,x)](x) −σ(ti,x)Ex tih2α(Xπ ti+1) + βi −σ2(ti,x)Ex tih∆Wi2α(Xπ ti+1) + βi (using the duality formula) =σ(ti,x)F−1[F[˜ uα,β i+1(., βj,1)](ν)iνζi(ν,x)φi(ν,x)](x) −σ(ti,x)h2αx+∆ia(ti,x) + ∆iσ2(ti,x)+βi(123) using the duality formula once again. The implementation of higher order time discretization for FBSDEs on the alternative grid is described in the following algorithm. Algorithm 2produces the numerical intermediate solutions {ui,j,k}NiN k=0 , {˜ ui,j,k}NiN k=0 and {˙ ui,j,k}NiN k=0 at time step ti , 0 ≤i<n and stage j , 1 ≤j≤q+ 1 for the approximate solution ui , the intermediate solution ˜ ui and the approximate gradient ˙ ui, respectively, i=0, 1, . . . , n−1. J. Risk Financial Manag. 2022,15, 388 21 of 32 Algorithm 2 Fourier Interpolation Method on Alternative Grid for q -stage Runge–Kutta schemes 1. Discretize the restricted real space [x0−Nnl 2 , x0+Nnl 2] and the restricted Fourier space [−L 2 , L 2] with NnN space steps so to have the real space nodes {xnk}NnN k=0 and {νnk}NnN k=0 2. Value un(xnk) = g(xnk) 3. For any ifrom n−1 to 0 (a) For any j, 1 <j≤q+1 i. Compute αand βdefining the transform of Equation (44), such that θi+1=˜ uα,β i+1(., αj,1)(124) and θi+1satisfies the boundary conditions of Equations (32) and (33). ii. Compute θi(xik)through Equation (40) for k=0, 1, . . . , NiNwith ψ(ν,x) = φi(ν,x)(125) and retrieve the values ˜ ui,j,kwith the appropriate correction. iii. Compute αand βdefining the transform of Equation (44), such that θi+1=˜ uα,β i+1(., βj,1)(126) and θi+1satisfies the boundary conditions of Equations (32) and (33). iv. Compute θi(xik)through Equation (40) for k=0, 1, . . . , NiNwith ψ(ν,x) = Φi,j(ν,x)(127) and retrieve the values ˙ ui,j,kwith the appropriate correction. v. Compute the values ui,j,kas ui,j,k=˜ ui,j,k+∆i j ∑ s=2 αjs f(ti,s,xi,k,ui,s,k,˙ ui,s,k)(128) for k=0, 1, . . . , NiNthrough Equation (13). vi. Update the real space grid with Equation (27) and the Fourier space grid by discretizing the interval [−L 2 , L 2] with NiN space steps so to have the real space nodes {xik}NiN k=0and {νik}NiN k=0. (b) Set ui,k=ui,q+1,kand ˙ ui,k=˙ ui,q+1,k Algorithm 2, as compared to Algorithm 1, contains intermediate solutions between time discretization points, resulting from the the intermediate stages of the Runge–Kutta method of Chassagneux and Crisan (2014). These intermediate solutions are computed using the Fourier interpolation method on the alternative grid presented in Section 2. We close this section by examining the spatial discretization error. 3.4. Spatial Discretization Error Analysis We denote by {ui,j,k}NiN k=0 and {˙ ui,j,k}NiN k=0 the intermediate numerical solutions obtained at time step ti , i= 0, 1, . . . , n− 1 and stage j , 1 <j≤q+ 1, from the Fourier interpolation method on the alternative grid when using a q− stage Runge–Kutta scheme. In addition, {ui,j,k}NiN k=0 and {˙ui,j,k}NiN k=0 are the intermediate numerical solutions obtained at the inter- J. Risk Financial Manag. 2022,15, 388 22 of 32 mediate stage j , 1 <j≤q+ 1, of time step ti given the exact solutions ui+1 and ˙ ui+1 at ti+1 . We have from the notation previously used that the numerical solutions at tiwrite ui,k=ui,q+1,k(129) ˙ ui,k=˙ ui,q+1,k(130) and are computed from the intermediate solutions {˜ ui,k}NiN k=0 , 0 <i≤n where ˜ un,k=˜ un(xn,k). When the exact solutions ui+1and ˙ ui+1are known at ti+1, we also write ui,k=ui,q+1,k(131) ˙ui,k=˙ui,q+1,k. (132) The local (space) discretization error has the form Eik :=ui(xk)−ui,k+˙ ui(xk)−˙ui,k(133) for i= 0, 1, . . . , n− 1 and k= 0, 1, . . . , NiN . The next theorem gives a description of the local (space) discretization error bound. Theorem 7. Suppose that the driver f∈ C1,2([ 0, T]×R2) and the terminal condition g∈ C2(R) and Assumptions 1and 3are satisfied, then the Fourier interpolation method yields a local space discretization error of the form sup i,k Eik =O(∆x)+Oe−K∆−s0 ilq0(134) for some constant K> 0on the alternative grid and under the trapezoidal quadrature rule for any explicit q-stage Runge–Kutta scheme. Proof. We follow the steps in the proof of Theorem 2. The truncation error when computing the numerical solutions ˙ui,j,kis Exik tihHϕj ti,j,γj∆i˜ ui+1(ti+1,Xπ ti+1;βj,1)1|∆Xπ i|>l 2i <KExik tiHϕj ti,j,γj∆i41 4Exik tih1|∆Xπ i|>l 2i1 4 (using Cauchy-Schwarz inequality twice since ˜ ui+1(ti+1,Xπ ti+1; .)is sq. int.) <K∆−1 2 iExik tih1|∆Xπ i|>l 2i1 4(since Hϕj ti,j,γj∆iis of Gaussian distribution) ≤K∆−1 2 iinf s>0e−sl 2φi(−is) + inf s>0e−sl 2φi(is)1 4(by Chernoff’s inequality) <K∆−1 2 ie−K0∆−s0 ilq0(by Assumption 3) <Ke−C∆−s0 ilq0. The Fourier interpolation leads to a first-order space discretization error when computing the numerical solutions ˙ui,j,k since the driver f and the terminal condition g are twice differentiable. The same statements hold for the numerical solutions ui,2,k using similar arguments. By recursion and using the Lipschitz property of the driver f , the statements hold for ui,j,k , 1 <j≤q+ 1. Since the time step ti and the space node xik are arbitrary, the space truncation and discretization error bounds hold for any iand k. J. Risk Financial Manag. 2022,15, 388 23 of 32 Locally, the truncation error remains spectral. Nonetheless, it is of a unspecified index q0 in this general setting where the conditional characteristic function φi is itself unspecified. For higher order time discretizations, one can expect q0≤ 2 since the forward process increment Xπ ti+1−Xπ ti has a heavy tail distribution. Indeed, the Gaussian distribution of forward process increments and the quadratic exponential form of their characteristic functions were the main reason for the spectral convergence of index 2 of the truncation error in Section 2.4. The space discretization error though is unchanged with first-order due to the second-order differentiability of the BSDE coefficients. However, the Fourier interpolation produces a space discretization error with a higher order when the driver f and the terminal function g have the required smoothness. In general, if f∈ Cm+1 b and g∈ Cm+1 b , we can expect a space discretization error of order m which is the convergence order of the underlying Fourier interpolation. We now turn to the global space discretization error defined as in Equation (63). The next theorem gives its error bound. Theorem 8. Suppose the conditions of Theorem 7are satisfied. If the discretization is such that sup i(C0∆x π∆p0 i)≤1 (135) then the Fourier interpolation method is stable and yields a global discretization error El,∆x of the form El,∆x=O(∆x) + Oe−K|π|−s0lq0(136) where K >0for any explicit q-stage Runge–Kutta scheme. Proof. From the definition of the global space discretization error, we may write eik ≤En−i,k+un−i,k−un−i,k(137) ˙ eik ≤En−i,k+˙un−i,k−˙ un−i,k. (138) Assume that the boundary values of the function ˜ ui+1 and the sequence ˜ ui+1,s are matched on the alternative grid so that we do not have to treat the alternative transform. Under an explicit q−stage Runge–Kutta scheme, we have ˙ui,j,k−˙ ui,j,k= D−1hΦi,j(νi+1,m,xik)D[˜ ui+1−˜ ui+1,s]mNi+1N−1 m=0ik+N 2 ≤∑Ni+1N−1 m=0Φi,j(νi+1,m,xik) Ni+1Nsup k˜ ui+1(xik,β1,j)−˜ ui+1,k ≤∆x 2πZRΦi,j(ν,xi,k)dνsup k˜ ui+1(xik,β1,j)−˜ ui+1,k ≤C0∆x 2π∆p0 i sup k˜ ui+1(xik,β1,j)−˜ ui+1,k(using Assumption 3) ≤C0∆x 2π∆p0 i (1+∆iK)sup k en−i−1,k+C0∆x 2π∆p0 i ∆iKsup k ˙ en−i−1,k (since fis Lipschitz and β1,jis bounded) ≤C0∆x 2π∆p0 i (1+∆iK)sup k en−i−1,k+C0∆x 2π∆p0 i (1+∆iK)sup k ˙ en−i−1,k. (139) J. Risk Financial Manag. 2022,15, 388 24 of 32 Similarly, we get ui,2,k−ui,2,k≤ D−1h{φi(νi+1,m,xik)D[˜ ui+1−˜ ui+1,s]m}Ni+1N−1 m=0ik+N 2 ≤∆x 2πZRφi(ν,xi,k)dνsup k˜ ui+1(xik,α1,2)−˜ ui+1,k ≤C0∆x 2π∆p0 i sup k˜ ui+1(xik,α1,2)−˜ ui+1,k(using Assumption 3) ≤C0∆x 2π∆p0 i (1+∆iK)sup k en−i−1,k+C0∆x 2π∆p0 i (1+∆iK)sup k ˙ en−i−1,k so that we get ui,j,k−ui,j,k≤C0∆x 2π∆p0 i (1+∆iK)"sup k en−i−1,k+2 sup k ˙ en−i−1,k#(140) recursively for 1 <j≤q+ 1 using the Lipschitz property of the driver f and the boundedness of the Runge–Kutta coefficients. Equations (137) and (138) combined with Equations (140) and (139) lead to sup k ei,k+sup k ˙ ei,k≤2 sup i,k Eik +C0∆x π∆p0 i (1+∆n−iK) sup k ei−1,k+sup k ˙ ei−1,k! ≤2 sup i,k Eik +ζ(1+∆n−iK) sup k ei−1,k+sup k ˙ ei−1,k! where sup i(C0∆x π∆p0 i)≤ζ≤1. Gronwall’s Lemma then yields sup k ei,k+sup k ˙ ei,k≤2eTK sup i,k Eik (141) so that the scheme is stable. The result of Equation (136) follows by taking the supremum on the left hand side of Equation (141) over time steps and applying Theorem 7. In this general case, the global discretization error maintains the structure of the local discretization error under a stability condition. Equation (135) indicates that the space discretization has to be relatively as fine as the time discretization to ensure stability. Hence, stability can always be reached for any time discretization by refining the space discretization. However, the structure of the characteristic functions φi and Φij determines the relative refinement needed for the space discretization. 4. Numerical Results We test the convergence properties of the Fourier interpolation method on Runge– Kutta schemes with a problem of commodity derivative pricing under a model proposed by Lucia and Schwartz (2002). We shall test the method’s convergence and behaviour on smooth and unbounded FBSDE coefficients. The commodity spot price Xis defined by Xt=eS(t)+Vt(142) J. Risk Financial Manag. 2022,15, 388 31 of 32 Data Availability Statement: Not applicable. Acknowledgments: The authors thank the anonymous referees for their helpful comments which improved the paper. Conflicts of Interest: The authors declare no conflict of interest. Notes 1This approach was suggested by an anonymous referee of an earlier version of this paper. 2The minimum sampling rate to avoid aliasing. 3The real value ¯ Pcan be considered as the production cost (per unit) of the commodity. 4See Equation (144). References Bally, Vlad and Gilles Pag ` es. 2003. A quantization algorithm for solving multidimensional discrete-time optimal stopping problems. Bernoulli 9: 1003–49. [CrossRef] Bally, Vlad, Gilles Pag ` es, and Jacques Printems. 2005. A quantization tree method for pricing and hedging multidimensional American options. Mathematical Finance 15: 119–68. [CrossRef] Beck, Christian, Weinan E, and Arnulf Jentzen. 2019. Machine learning approximation algorithms for high-dimensional fully nonlinear partial differential equations and second-order backward stochastic differential equations. Journal of Nonlinear Science 29: 1563–19. [CrossRef] Bender, Christian, and Jianfeng Zhang. 2008. Time discretization and Markovian iteration for coupled FBSDEs. The Annals of Applied Probability 18: 143–77. [CrossRef] Bender, Christian, and Robert Denk. 2007. A forward scheme for backward SDEs. Stochastic Processes and Their Applications 117: 1793–812. [CrossRef] Bouchard, Bruno, and Nizar Touzi. 2004. Discrete-time approximation and Monte-Carlo simulation of backward stochastic differential equations. Stochastic Processes and Their Applications 111: 175–206. [CrossRef] Briand, Philippe, Bernard Delyon, and Jean M ´ emin. 2001. Donsker-type theorem for BSDEs. Electronic Communications in Probability 6: 1–14. [CrossRef] Chassagneux, Jean-Franc¸ois, and Dan Crisan. 2014. Runge-Kutta schemes for BSDEs. The Annals of Applied Probability 24: 679–720. Chevance, David. 1997. Numerical methods for backward stochastic differential equations. In Numerical Methods in Finance. Edited by Leonard Christopher Gordon Rogers and Denis Talay. Publications of the Newton Institute. Cambridge: Cambridge University Press, pp. 232–44. Crisan, Dan, and Konstantinos Manolarakis. 2012. Solving backward stochastic differential equations using the cubature method: application to nonlinear pricing. SIAM Journal on Financial Mathematics 3: 534–71. [CrossRef] Crisan, Dan, and Konstantinos Manolarakis. 2014. Second order discretization of backward SDEs and simulation with the cubature method. The Annals of Applied Probability 24: 652–78. [CrossRef] Crisan, Dan, Konstantinos Manolarakis, and Nizar Touzi. 2010. On the Monte Carlo simulation of BSDEs: An improvement on the Malliavin weights. Stochastic Processes and Their Applications 120: 1133–58. [CrossRef] Delarue, Fran c¸ ois, and St ´ ephane Menozzi. 2006. A forward-backward stochastic algorithm for quasi-linear PDEs. The Annals of Applied Probability 16: 140–84. [CrossRef] Douglas, Jim, Jr., Jin Ma, and Philip Protter. 1996. Numerical methods for forward-backward stochastic differential equations. The Annals of Applied Probability 6: 940–68. [CrossRef] Duffie, Darrell, and Larry G. Epstein. 1992. Stochastic differential utility. Econometrica 60: 353–94. [CrossRef] E, Weinan, Jiequn Han, and Arnulf Jentzen. 2017. Deep learning-based numerical methods for high-dimensional parabolic partial differential equations and backward stochastic differential equations. Communications in Mathematics and Statistics 5: 349–80. [CrossRef] E, Weinan, Jiequn Han, and Arnulf Jentzen. 2022. Algorithms for solving high dimensional PDEs: From nonlinear Monte Carlo to machine learning. Nonlinearity 35: 278–310. [CrossRef] E, Weinan, Martin Hutzenthaler, Arnulf Jentzen, and Thomas Kruse. 2019. On multilevel Picard numerical approximations for highdimensional nonlinear parabolic partial differential equations and high-dimensional nonlinear backward stochastic differential equations. Journal of Scientific Computing 79: 1534–71. [CrossRef] El Karoui, Nicole, ´ Etienne Pardoux, and Marie-Claire Quenez. 1997. Reflected backward SDEs and American options. In Numerical Methods in Finance. Edited by Leonard Christopher Gordon Rogers and Denis Talay. Publications of the Newton Institute. Cambridge: Cambridge University Press, pp. 215–31. El Karoui, Nicole, Shige Peng, and Marie-Claire Quenez. 1997. Backward stochastic differential equations in finance. Mathematical Finance 7: 1–71. [CrossRef] Gobet, Emmanuel, Jean-Phillipe Lemor, and Xavier Warin. 2005. A regression-based Monte Carlo method to solve backward stochastic differential equations. The Annals of Applied Probability 15: 2172–202. [CrossRef] J. Risk Financial Manag. 2022,15, 388 32 of 32 Han, Jiequn, and Jihao Long. 2020. Convergence of the deep BSDE method for coupled FBSDEs. Probability, Uncertainty and Quantitative Risk 5: 5. [CrossRef] Huijskens, Thomas P., Marjon J. Ruijter, and Cornelis W. Oosterlee. 2016. Efficient numerical Fourier methods for coupled forwardbackward SDEs. Journal of Computational and Applied Mathematics 296: 593–612. [CrossRef] Hyndman, Cody Blaine, and Polynice Oyono Ngou. 2017. A convolution method for numerical solution of backward stochastic differential equations. Methodology and Computing in Applied Probability 19: 1–29. [CrossRef] Kloeden, Peter E., and Eckhard Platen. 1992. Numerical Solution of Stochastic Differential Equations. Applications of Mathematics (New York). Berlin: Springer, vol. 23. Lucia, Julio J., and Eduardo S. Schwartz. 2002. Electricity prices and power derivatives: Evidence from the Nordic power exchange. Review of Derivatives Research 5: 5–50. [CrossRef] Ma, Jin, Jie Shen, and Yanhong Zhao. 2008. On numerical approximations of forward-backward stochastic differential equations. SIAM Journal on Numerical Analysis 46: 2636–61. [CrossRef] Ma, Jin, Philip Protter, Jaime San Mart ´ ın, and Soledad Torres. 2002. Numerical method for backward stochastic differential equations. The Annals of Applied Probability 12: 302–16. Oyono Ngou, Polynice. 2014. Fourier Methods for Numerical Solution of FBSDEs with Applications in Mathematical Finance. Ph.D. thesis, Concordia University, Montr´ eal, QC, Canada, January. Pardoux, ´ Etienne, and Shige Peng. 1992. Backward stochastic differential equations and quasilinear parabolic partial differential equations. In Stochastic Partial Differential Equations and Their Applications (Charlotte, NC, 1991). Lecture Notes in Control and Information Sciences. Berlin: Springer, vol. 176, pp. 200–17. Peng, Shige, and Mingyu Xu. 2011. Numerical algorithms for backward stochastic differential equations with 1-d Brownian motion: Convergence and simulations. ESAIM: Mathematical Modelling and Numerical Analysis 45: 335–60. [CrossRef] Ruijter, Marjon J., and Cornelis W. Oosterlee. 2015. A Fourier-cosine method for an efficient computation of solutions to BSDEs. SIAM Journal on Scientific Computing 37: A859–89. [CrossRef] Ruijter, Marjon J., and Cornelis W. Oosterlee. 2016. Numerical Fourier method and second-order Taylor scheme for backward SDEs in finance. Applied Numerical Mathematics 103: 1–26. [CrossRef] Turkedjiev, Plamen. 2015. Two algorithms for the discrete time approximation of Markovian backward stochastic differential equations under local conditions. Electronic Journal of Probability 20: 1–49. [CrossRef] Vaˇ s´ ıˇ cek, Oldˇ rich. 1977. An equilibrium characterization of the term structure. Journal of Financial Economics 5: 177–88. [CrossRef] Zhang, Jianfeng. 2004. A numerical scheme for BSDEs. The Annals of Applied Probability 14: 459–88. [CrossRef]