University of Southern California Dissertations and Theses

Numerical methods for high-dimensional path-dependent PDEs driven by stochastic Volterra integral equations

(USC Thesis Other)

Numerical methods for high-dimensional path-dependent PDEs driven by stochastic Volterra integral equations

PDF

Download a page range

Download transcript

Copy asset link

Request this asset

Transcript (if available)

Content NUMERICAL METHODS FOR HIGH-DIMENSIONAL PATH-DEPENDENT PDES DRIVEN BY STOCHASTIC VOLTERRA INTEGRAL EQUATIONS by Jie Ruan A Dissertation Presented to the FACULTY OF THE USC GRADUATE SCHOOL UNIVERSITY OF SOUTHERN CALIFORNIA In Partial Fulllment of the Requirements for the Degree DOCTOR OF PHILOSOPHY (APPLIED MATHEMATICS) August 2020 Copyright 2020 Jie Ruan Dedication To my family ii Acknowledgments I used to be arrogant enough to think Ph.D or even life is just a matter of personal hard work and dedication, which I thought I certainly have as a naive undergrad- uate. Not until I encountered those unexpected setbacks during my Ph.D, did I realize everything I was proud of would not have been possible without being lucky enough to be helped by others. So I take this opportunity to thank those people I am indebted to in this journey. I owe my deepest gratitude to my advisor, Prof. Jianfeng Zhang, for his patient guidance on my research and genuine care for my life. Not only did he shape my mathematical perspective with his wise style of approaching problems and down- to-earth working attitude, he also helped me through some most dicult periods in my life. During our conversations, he often impresses me with the insights of a mathematician and touches me with the sincerity of a generous person. For me, he is the living testament of the old Chinese saying \One day as your teacher, like a father for a lifetime". iii I also greatly appreciate Prof. Jin Ma, a respectable and cheerful mathemati- cian who treats me like family and oered me unreserved support when I expe- rienced troubles and confusions for the future. Besides, my sincere thanks shall go to Prof. Jinchi Lv for serving on my dissertation committee, Prof. Remigi- jus Mikulevicius and Stanislav Minsker for serving my qualifying committee, Prof. Peter Baxendale, Larry Goldstein, Stanislav Minsker, Igor Kukavica and Sergey Lototsky for sharing their knowledge with me through valuable courses. In addi- tion, I am especially grateful to Prof. Yijun Yao for leading me to this wonderful journey. I am also indebted to my seniors Jia Zhuo, Cong Wu, Jian Wang, Fanhui Xu, Fei Wang, Weisheng Xie, Rentao Sun, Eunjung Noh, and my classmates Chukiat Phonsom, Jiyeon Park, Zheng Dai, Bowen Gang, Jiajun Luo, Jiaowen Yang, Zimu Zhu, Pengbin Feng, Man Luo, Wenqian Wu, Linfeng Li, Ying Tan for their warm- hearted helps and the joy we shared along the way, which I haven't had the chance to say thanks. Lastly, my deepest love goes to my parents Minxian Lou and Zhangqing Ruan, who persevere throughout their lives despite their limited educational privilege so that I can see the bigger world and live a better life. iv Table of Contents Dedication ii Acknowledgments iii Abstract vi Chapter 1: Introduction 1 Chapter 2: Preliminary 11 2.1 Probability Setup, Notations and Assumptions . . . . . . . . . . . . 11 2.2 Classical Results in PPDEs and BSDEs . . . . . . . . . . . . . . . . 15 2.3 Deep Learning Based Method . . . . . . . . . . . . . . . . . . . . . 17 Chapter 3: Stochastic Volterra Integral Equations 23 3.1 Wellposedness . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24 3.2 Regularity and Flow Property . . . . . . . . . . . . . . . . . . . . . 40 Chapter 4: Branching Diusion Method 47 4.1 Branching Processes Setup . . . . . . . . . . . . . . . . . . . . . . . 48 4.2 Representation Results . . . . . . . . . . . . . . . . . . . . . . . . . 51 4.3 Blow-up Issues and Transformations . . . . . . . . . . . . . . . . . . 59 Chapter 5: Multilevel Picard Approximation 65 5.1 Convergence of Numerical Scheme . . . . . . . . . . . . . . . . . . . 66 5.2 Computation Cost and Selection of . . . . . . . . . . . . . . . . . 76 Chapter 6: Numerical Examples 84 Bibliography 96 v Abstract In this dissertation we study the stochastic Volterra equations and numerical meth- ods for the associated path-dependent PDEs (PPDEs for short). It consists of two main parts. We rst prove the wellposedness of path-dependent stochastic Volterra equations, and generalize the estimate results and ow properties of the classical SDEs to the Volterra ones. In the second part, we extend three dierent numerical methods for high-dimensional PDEs to our path-dependent scenarios: branching diusion method, multilevel Picard approximation and deep learning based method. For the branching diusion method, we explain the blow-up issues and introduce a transformation to extend the time duration, which is especially eective in quadratic cases. For the multilevel Picard approximation, we introduce an extra parameter in stopping strategy, which will help extend the time duration. In the end, we present several high-dimensional numerical examples with dierent emphasis, and give a comparison of the three methods as well as the guidelines about which method to choose in dierent scenarios, in terms of their applicability, scalability and reliability. vi Chapter 1 Introduction Developing ecient numerical methods for path-dependent PDEs (PPDEs, for short) with high dimensions is a challenging topic with potentially signicant impacts. The non-Markovian version of many well-known problems, including the pricing of path-dependent derivatives and the value processes of path-dependent stochastic control, have their solutions represented through PPDEs, which are essentially path-dependent Black-Scholes equations and Hamilton-Jacobi-Bellman equations for these two examples. In the semilinear cases for example, @ t ub@ x u 1 2 Tr[ T @ 2 xx u]f(t; x;u; T @ x u) = 0; u(T; x) =g(x): (1.1) where the time and spatial derivative are properly dened, the solution u to the PPDE is associated with the solutions to the following decoupled path-dependent forward-backward SDEs (FBSDEs, for short) in the sense that Y t;x s =u(s;X t;x ^s ). 8 > > > < > > > : X t;x s = x s + R s t b(r;X t;x ^r )dr + R s t (r;X t;x ^r )dW r Y t;x s =g(X t;x ) + R T s f(r;X t;x ^r ;Y t;x r ;Z t;x r )dr R T s Z t;x r dW r : (1.2) 1 Such connection was rst proved to hold between the standard semilinear PDEs and the Markovian BSDEs in Peng [36] and Pardoux-Peng [38], known as the nonlinear Feynman-Kac formula. It was later extended to between quasilinear PDEs and coupled FBSDEs, see for example, Ma-Protter-Yong [33], Pardoux-Tang [40] and Ma-Zhang-Zheng [34], and between fully nonlinear PDEs and second order BSDEs (2BSDEs for short), see for example, Cheridito et al. [12] and Soner-Touzi- Zhang [44]. It was in light of this fundamental equivalence that, Peng rst raised the question of whether some notion of path-dependent PDEs (PPDEs, for short) can be proposed to connect with the non-Markovian BSDEs, see his ICM2010 talk [41]. Interestingly Dupire [13], who are motivated from a more practical perspec- tive, e.g. pricing problems related to path-dependent options, rst introduced the insightful notion of horizontal and vertical derivatives, and generalized the classi- cal It^ o formula to functionals on paths, which lay the ground for PPDEs. Based on this work and later renement by Cont-Fournie [10], a complete viscosity the- ory for PPDEs was successfully developed in Ekren-Keller-Touzi-Zhang [18] and Ekren-Touzi-Zhang [21, 22, 23], which generalized the connection to that between PPDEs and non-Markovian problems. Moreover, recent progresses in volatility modeling suggest the path-dependent feature in PPDEs are intrinsic and ubiquitous in derivative pricing and hedging 2 problems. In the nancial world, the log-prices log(S t ) of a given asset are typically modeled as continuous semimartingale taking the form d log(S t ) = t dt + t dW t ; where the term t is the volatility process. As it is the most important ingredient, various models have been proposed for the volatility process but most of them still bear noticeable discrepancies from the empirical data. Recently, a revealing analysis of the volatility time series, conducted by Gatheral-Jaisson-Rosenbaum [26] using high frequency data, showed that log-volatility behaves essentially as a fractional Brownian motion (fBM) with Hurst exponent H of order around 0.1. This has motivated the proposal of rough fractional stochastic volatility (FSV) models in [26] and the rough Heston Models in [20], which turn out to be remark- ably consistent with nancial time series data in terms of regularity, shape of the volatility surface and the stylized long memory property. In these models, either log( t ) or 2 t are given as solutions to the stochastic Volterra integral equations (SVIEs) in the following form X s = x s + Z s t b(s;r;X r )dr + Z s t (s;r;X r )dW r ; (1.3) of which the fractional Brownian motion (fBM) is the most well-known example. This class of processes is a non-trivial generalization of classical SDEs, because 3 they are neither Markovian nor semi-martingales in general. Most importantly, they are intrinsically path-dependent, for example even the simplest conditional expectation for a state dependent function g Y t :=E[g(X T )jF t ] will now depends on the entire path of X up to time t. Therefore, any BSDE associated with such forward processes will become path dependent automatically, even if it is Markovian itself, i.e. the g = g(X T ) and f = f(t;X t ;Y t ;Z t ) in (1.2). Indeed, in Viens-Zhang [45], they extended Dupire's functional It^ o formula for SVIEs and showed the solutions to the associated BSDEs are also represented through a corresponding PPDEs. The main objective of this research is to develop ecient numerical methods for PPDEs driven by SVIEs, which will apply to PPDEs associated with standard SDEs as a special case. Until recently very few researches have gone into the numerical methods for PPDEs, see Guo-Zhang-Zhuo [28], Ren-Tan [43], as opposed to the intensive studies on the numerical methods for standard PDEs in high dimensions. Apart from its being a relatively new concept, another main reason is the path-dependent or non-Markovian feature of PPDEs brings extra "high dimensionality" along time, which make it more dicult than standard ones. As for SVIEs, the regular kernel cases were rst studied by Berger-Mizel [4] and the 4 wellposedness was proved for the form (1.3). Some other directions were also explored later, see for example, Pardoux-Protter [39] for SVIEs with anticipating coecients. The wellposedness for singular kernel cases could be found in Coutin- Decreusefond [9] for a special form and Wang [46] for the more general form (1.3). And a corresponding Euler schemes were developed in Zhang [47]. In chapter 3, we prove the wellposedness of more general SVIEs with path-dependent and singular coecients (1.4), which unies all the existing wellposedness results. X s = x s + Z s t b(s;r;X ^r )dr + Z s t (s;r;X ^r )dW r ; (1.4) The main diculty is bounding the path norm when the Burkholder-Davis-Gundy (BDG) inequality in classical arguments fails and it is overcome with an embedding inequality in the same spirit as the Kolmogorov- Centsov Theorem. More impor- tantly, we recover the ow property based on ideas in Viens-Zhang [45], which will be crucial for developing numerical methods for associated PPDEs. We did not realized the coecient path-dependency was also considered in Protter [42] until the thesis was nished, but our assumptions are weaker and cover the fBM type singular kernels, which are of the main interest for rough volatility models in practice. Developing ecient numerical methods for high-dimensional PDEs, even just the standard ones, has already well known to be dicult because of the curse 5 of dimensionality. Classical methods based on spatial discretization in the PDEs literature, e.g. nite dierence methods and nite elements methods, typically become impractical ford 3. Thus researchers have either turned to probabilistic methods or tried to exploit the special structure of particular type of PDEs. The rst successful attempt was the backward Euler scheme for semilinear PDEs by Zhang [48] and Bouchard-Touzi [6]. Based on the nonlinear Feynman- Kac formula, they represented the solution through the Markovian version of (1.2) and approximated the BSDE by a sequence of conditional expectations, 8 > > > < > > > : Y i =E i [Y i+1 ] +f(i;X i ;Y i ;E i [ W i+1 W i t i+1 t i Y i+1 ]) Y n =g(X n ) (1.5) from which they obtained the rate of convergence. Ecient numerical algorithms were proposed to compute the conditional expectations, see for example, Bouchard- Touzi [6], Gobet-Lemor-Warin [27], Crisan-Manolarakis [11], Bally-Pag es-Printems [5] and Bender-Denk [2]. The scheme was also extended to fully nonlinear cases in Fahim-Touzi-Warin [24] and further to PPDEs in Guo-Zhang-Zhuo [28], but only practical for very low dimension. We shall remark that this is far from being a comprehensive review here, as there has been a large stream of literature around this scheme. Essentially all of them have to deal with the computation of condi- tional expectation, which is the main obstacle for better eciency. Overall this class of methods works for dimensions around 10. 6 Henry-Labord ere [30] rst proposed the Branching diusion method by recog- nizing the forward nature of classical representation for KPP equations see McKean [35]. Unlike the backward Euler schemes, the representation exploits the polyno- mial structure and gracefully replaces the nested conditional expectation by a purely forward simulation, which is perfectly suited for high dimension scenarios. The work was later extended to non-Markovian BSDEs, see Henry-Labord ere-Tan- Touzi [32], and to allow f being polynomial in @ x u as well, see Henry-Labord ere et al [31]. Strongly motivated by their work, we extend the Branching diusion method to PPDEs driven by SVIEs in Chapter 4, using the new ow property we derived, and obtain the representation results. One of our main contributions is to observe the absolute value of the representation corresponds to a transformed PDE, which is typically wellposed only for a smaller time duration, and this explains the blow-up issues encountered in those papers. Moreover, we propose a transforma- tion to increase the time duration and prove it always work for quadratic cases. It is possible to use more general transform but still retain the polynomial structure by enlarging the original PDE to a PDE system. This might generalize the approach for higher degree polynomials but will be left for future works. Although it requires the polynomial structure to work, the Branching diusion method is remarkably ecient and unrivaled when applicable. So it is tempting to approximate gen- eral nonlinearity with polynomials and apply the Branching method, which was explored in Bouchard-Tan-Warin [7] and Bouchard et al [8]. Convergence was 7 obtained for their adapted scheme but the curse of dimensionality reappear due to interpolation. In E et al [16, 17], they proposed the multilevel Picard approximations based on the nonlinear Feynman-Kac formula. Instead of straightforwardly computing the nested expectation, they ingeniously break up the representation according to standard Picard approximation and allocate simulations proportional to the variance of each piece. This reduce the computation fromO(M n 2 ) toO(M n ) while achieving essentially same error order by exploiting the quick convergence of Picard approximation. In practice, it was shown to work for dimensions around 100 but subject to small Lipschitz constants or equivalently small time duration. In chapter 5, we extend this approximation schemes to our PPDEs scenarios and obtain a ner rate of convergence with a much conciser proof. The other contribution of us is to introduce the exponential stopping strategy, for which the range of optimal parameter is also derived, to reduce the computation cost or equivalently increase the time duration. Another promising approach is the so-called deep learning based method. It was rst proposed for the stochastic control problem, see E-Han [14], and later recognized to be applicable to PDEs through the BSDEs formulation, see E-Han- Jentzen [15] for semilinear cases and Beck-E-Jentzen [3] for fully nonlinear cases. The main idea is to view BSDEs formulation of the PDE as a forward control problem, where the Z serves as the control process and can be approximated as 8 a function of X due to the nonlinear Feynman-Kac. In chapter 2, we extend the deep learning approach to the SVIEs settings and the crucial distinction is the Z is now approximated based on the path-dependent nonlinear Feynman-Kac derived in Viens-Zhang [45]. Empirically deep neural networks work well for the function approximation part though other simpler models could also be explored for further eciency. Compared to the backward Euler scheme with linear regression in Gobet-Lemor-Warin [27], the deep learning approach has a similar spirit but achieves much better eciency, e.g. applicable to dimension 100, by optimizing in a forward global instead of backward stepwise manner, which becomes possible thanks to the advances in computation and model tting techniques. The rest of the dissertation is organized as follows. In Chapter 2, we introduce the probability setup, notations, key assumptions, fundamental results in PPDEs and BSDEs, and our extension of the deep learning based method to PPDEs. In Chapter 3, we prove the wellposedness for path-dependent SVIEs, regularity of the solutions and the ow property. Chapter 4 is dedicated to the extension of Branching diusion method to PPDEs associated with SVIEs. We introduce the Branching processes and prove the representation results. Our observation about the blow-up issues and solutions are also provided in detail. Chapter 5 deals with extension of the multilevel Picard scheme to PPDEs. We prove the rate of convergence and estimate the computation cost, on the basis of which the selection 9 of optimal parameters is discussed. In Chapter 6, we compare three methods on numerous examples and discuss their performance from dierent perspectives. 10 Chapter 2 Preliminary In this section, we provide the common probability setup, notations and assump- tions in the thesis, while leaving local ones in their corresponding chapter for better reference. Classical results in their own literature are stated without proof. In the last section, we also present our extension of deep learning based method to PPDEs and give the numerical scheme in detail. 2.1 Probability Setup, Notations and Assump- tions Let T be the xed time horizon and ( ;F;P) be a probability space equipped with ad-dimensional Brownian motion and the associated ltration (W t ;F t ) 0tT , i.e. F t =F W t . Moreover, the space is assumed to be rich enough to allow inclu- sion of any additional random variables independent of W whenever needed. The probability space is assumed implicitly throughout the thesis to save referring to it repeatedly. The following notations will be used throughout the thesis: 11 D([0;T ];R d ): c adl ag mappings from [0;T ] toR d C([0;T ];R d ): continuous mappings from [0;T ] toR d L 2 ([0;T ];R d ): F t adapted processes ' : [0;T ] ! R d such that E R T 0 j' t j 2 dt<1 := [0;T ]C([0;T ];R d ) := (t; x)2 [0;T ]D([0;T ];R d ) : xj [t;T ] 2C([t;T ];R d ) jxj := sup 0tT jx t j for x2D([0;T ];R d ) d((t; x); (t 0 ; x 0 )) :=jtt 0 j +jx x 0 j for (t; x); (t 0 ; x 0 )2D([0;T ];R d ) @ t u(t; x) := lim #0 u(t+;x)u(t;x) for (t; x)2 h@ x u(t; x);i := lim "!0 u(t;x+"1 [t;T] )u(t;x) " for 2C([0;T ];R d ) h@ 2 xx u(t; x); ( 1 ; 2 )i := lim "!0 h@xu(t;x+" 2 1 [t;T] ); 1 ih@xu(t;x); 1 i " for 1 ; 2 2 C([0;T ];R d ) 12 C 1;2 + ( ): functionals u : ! R with continuous @ t u;@ x u;@ 2 xx u w.r.t d(;) and there exists C; > 0 and a bounded modulus of continuity function such that, for any (t; x); (t; y)2 and ; 1 ; 2 2C([0;T ];R d ) j@ t u(t; x)jC(1 +jxj ); jh@ x u(t; x);ijC(1 +jxj )j1 [t;T ] j jh@ 2 xx u(t; x); ( 1 ; 2 )ijC(1 +jxj )j 1 1 [t;T ] jj 2 1 [t;T ] j jh@ 2 xx u(t; x)@ 2 xx (t; y); (;)ij (1 +jxj +jyj)j1 [t;T ] j 2 (jx yj) C 1;2 + () := uj ;u2C 1;2 + ( ) (x t y) u = x u 1 f0utg + y u 1 ft<uTg (X ) t :=X t , (X ^s ) t :=X t 1 ftsg +X s 1 ft>sg The main object of this thesis is the PPDEs driven by SVIEs, which could be equivalently written in a concise decoupled FBSDE form, 8 > > > < > > > : X t;x s = x s + R s t b(s;r;X t;x ^r )dr + R s t (s;r;X t;x ^r )dW r Y t;x s =g(X t;x ) + R T s f(r;X t;x ^r ;Y t;x r ;Z t;x r )dr R T s Z t;x r dW r : (2.1) The following assumptions about it will be made frequently throughout the thesis 13 Assumption 2.1.1. (b;) : [0;T ] [0;T ]C([0;T ];R d )! R d R dn are measurable w.r.t the-algebra induced by usual norm, and there existsK 1 ;K 2 such that for any x; x2C([0;T ];R d ), ts 0 ;s2 [0;T ], the following holds a.s. jb(s;t; 0)j +j@ s b(s;t; 0)jK 1 (s;t); j(s;t; 0)j +j@ s (s;t; 0)jK 2 (s;t) jb(s;t; x)b(s;t; x)j +j@ s b(s;t; x)@ s b(s;t; x)jK 1 (s;t)jx xj j(s;t; x)(s;t; x)j +j@ s (s;t; x)@ s (s;t; x)jK 2 (s;t)jx xj and I p p := sup 0sT Z s 0 K 1 (s;t) +K 2 2 (s;t)dt p <1 (2.2) for some > 1 and any p> 0. Assumption 2.1.2. g(x) : C([0;T ];R d ) ! R and f = f(t; x;y) : [0;T ] C([0;T ];R d )R! R are measurable w.r.t the usual -algebra, f is uniformly Lipschitz ony with Lipschitz constantL = 1, andjg(x)j+jf(t; x; 0)jC(1+jxj ) for some 1. Remark 2.1.3. The purpose of the assumption L = 1 is not only to simplify constant factors in the estimates, but also to emphasize the equivalent role that L andT play in the convergent rate. In fact, increasing(decreasing) theL by a factor c is equivalent to increasing(decreasing) the T by the same factor in our estimate. 14 2.2 Classical Results in PPDEs and BSDEs We state some standard results from PPDEs and BSDEs literature, which are used either explicitly or implicitly throughout the thesis. All BSDEs and SVIEs are assumed to start with (0; x) to reduce superscripts in the notation. Consider the following BSDE: Y t = + Z T t f(s;Y s ;Z s )ds Z T t Z s dW s (2.3) where f : [0;T ]RR d !R is progressively measurable. Assumption 2.2.1. (i) 2L 2 (F T ) and f(; 0; 0)2L 2 ([0;T ];R) (ii) f is uniformly Lipschitz continuous in (y;z), i.e. there exists C > 0 such that for any (y 1 ;z 1 ), (y 2 ;z 2 ) jf(!;t;y 1 ;z 1 )f(!;t;y 2 ;z 2 )jC(jy 1 y 2 j +jz 1 z 2 j) dP dt a:s: Theorem 2.2.2. (BSDE Wellposedness, [19]) Under Assumption 2.2.1, the BSDE (2.3) has a unique solution (Y;Z)2L 2 ([0;T ];R)L 2 ([0;T ];R d ). Consider the following SVIE and a related process: X t = x t + Z t 0 b(t;r;X ^r )dr + Z t 0 (t;r;X ^r )dW r (2.4) t s := x t + Z t^s 0 b(s;r;X ^r )dr + Z t^s 0 (s;r;X ^r )dW r (2.5) 15 Assumption 2.2.3. (i) @ t b(t;r;);@ t (t;r;) exists for t2 [r;T ], and for ' = b;;@ t b;@ t j'(t;r; x)jC[1 +jxj ] for some constants C;> 0; (ii) SVIE (2.4) has a weak solution (X;W ) such that E[jX j p ]<1 for all p 1. Theorem 2.2.4. (Functional It^ o formula for SVIEs, [45]) Let Assumptions 2.2.3 hold and u2C 1;2 + (), then du(t;X t t ) =@ t u(t;X t t )dt + 1 2 h@ 2 xx u(t;X t t ); ( t;X ; t;X )idt +h@ x u(t;X t t );b t;X idt +h@ x u(t;X t t ); t;X idW t (2.6) where ' t;x s :='(s;t; x) emphasizes the dependence on s2 [t;T ] for ' =b;. We only state the regular kernel case here for simplicity, for singular kernel case see [45] as well. With this generalized functional It^ o formula, one can easily verify the nonlinear Feynman-Kac, i.e. Y t :=u(t;X t t ), Z t :=h@ x u(t;X t t ); t;X i 16 solves the BSDE in (2.1), provided the existence of the solution the semilinear PPDE: @ t u(t; x) + 1 2 h@ 2 xx u(t; x); ( t;x ; t;x )i +h@ x u(t; x);b t;x i +f(t; x;u(t; x);h@ x u(t; x); t;x i) = 0; (t; x)2 ; (2.7) u(T; x) =g(x): therefore it is equivalent to talk about the above PPDE and the FBSDE (2.1). 2.3 Deep Learning Based Method In this section, we focus on extending the deep learning based method to the semilinear PPDE (2.7), where the essential dierence from Markovian case is h@ x u(t; x); t;x i has to be approximated as a function of the path instead of only the state. The scheme starts with the equivalent BSDE formulation: Y t =g(X ) + Z T t f(r;X ^r ;Y r ;Z r )dr Z T t Z r dW r : 17 whereX is the solution to the SVIE (2.4). In standard BSDEs, Y T is known, then the solution Y and Z are obtained by solving the equation backwardly, so the following holds if we view the dynamic forwardly Y t =Y 0 Z t 0 f(r;X ^r ;Y r ;Z r )dr + Z t 0 Z r dW r ; Y T =g(X ): Moreover, we knowZ t =h@ x u(t;X t t ); t;x i from the nonlinear Feynman-Kac. Now suppose we already know the Y 0 and u, by solving the following SDE with y = Y 0 , v(t; x) =h@ x u(t; x); t;x i, we should be able to recover Y , in particular g(X ) =Y Y 0 ;h@xu;i T . Y y;v t =y Z t 0 f(r;X ^r ;Y y;v r ;v(r;X r r ))dr + Z t 0 v(r;X r r )dW r (2.8) In this sense, we can view the PPDE solution u(t; x)2 C 1;2 + () and its Gateux derivativeh@ x u(t; x); t;x i as the solution to a stochastic control problem with terminal targeting. More formally, under suitable regularity assumptions onf and g, the uniqueness of the BSDE implies that (u(0; x);h@ x u(t; x); t;x i) = arg min y;v E[jY y;v T g(X )j 2 ]: (2.9) Then the problem of solving u, in particular u(0; x), becomes a standard statisti- cal learning problem, as suggested by the RHS of (2.9). By sampling paths from 18 dynamics ofX and updatingy andv(t; x) to minimize the empirical square error, one can expect to approximate the solution well as the sample size and discretiza- tion steps are suciently large. More specically, we rst apply a time discretization to (2.4), let N2N and t 0 ;t 1 ;:::;t N 2 [0;T ] be real numbers such that 0 = t 0 < t 1 < ::: < t N = T , approximateX at discretization pointsfX tn g 0nN using standard Euler schemes, i.e. X t n+1 := x t n+1 + n X k=0 h b t n+1 ;t k ;X [t 0 ;t k ] (t k+1 t k ) + t n+1 ;t k ;X [t 0 ;t k ] (W t k+1 W t k ) i ; 0n<N (2.10) whereX t 0 = x t 0 andX [t 0 ;tn] (s) := P n1 k=0 X t k 1 ft k s<t k+1 g +X tn 1 ft k sTg . For approx- imation of Y y;v , we can parameterize any candidate solution y and v of (2.9) in the following form yU 0 2R; v(t n ;X [t 0 ;tn] )V n (fX t k g 0kn ; n ) :R dn !R d ; 0n<N where =f n g 0nN1 2R for some2N are the full set of parameters and is the total dimension. This is true because eventually we deals with the discretized 19 dynamics. Then Y y;v can be approximated over the discretizationfY U 0 ; tn g 0nN in the standard Euler scheme as well: Y U 0 ; t n+1 :=Y U 0 ; tn f t n ;X [t 0 ;tn] ;Y U 0 ; tn ;V n (fX t k g 0kn ; n ) (t n+1 t n ) + V n (fX t k g 0kn ; n );W t n+1 W tn ; 0n<N: (2.11) Since Y Y 0 ;h@xu;i T = g(X ), one can expect E[ Y Y 0 ;h@xu;i t N g(X [t 0 ;t N ] ) 2 ] also to be close to 0 asN gets large. Therefore, it is natural to approximateY 0 andh@ x u;i by (U 0 ; ) := arg min U 0 ; E[ Y U 0 ; t N g(X [t 0 ;t N ] ) 2 ]; (2.12) assuming the discretized square error has some regularity over its minimizer. However, the recursion derivation ofX andY U 0 ; usually results in too compli- cated distribution for (2.12) to be optimized globally. Therefore, the minimization is typically carried out by local or greedy algorithms, e.g. some variants of stochas- tic gradient descent (SGD). For example, in the plain SGD with a learning rate and batch size B = 1, we starts with some initial (U 0 0 ; 0 ) and then update in the following form as we sample new pathfX m ;Y m;U 0 ; g fromX andY U 0 ; (U m 0 ; m ) = (U m1 0 ; m1 ) r U 0 ; Y m;U 0 ; t N g(X m [t 0 ;t N ] ) 2 ; 20 see [25] for a more systematic study on SGD. Such algorithms only requires the error function to be dierentiable almost everywhere, which is typically satised as long as f and g are regular enough, andfV n ( ; n )g 0n<N are chosen to be dierentiable. Empirically deep neural networks have shown to be good choices for fV n ( ; n )g 0n<N , which is why this method is called deep learning based method, but in theory one could use any function approximator. We summarize and put the algorithm with a simple batch SGD together as follows: Algorithm 1: Deep Learning Based Method for semilinear PPDEs (2.7) Choose deep neural networks structures for fV n (; n )g 0n<N ; :=f n g 0nN1 ; Choose learning rate sequencef m g 1mM and batch size B; Initialize (U 0 0 ; 0 ); for m = 1; 2;:::;M do Sample B paths offX j;m ;Y j;m;U m1 0 ; m1 g 1jB according to (2.10), (2.11); Update (U m 0 ; m ) (U m1 0 ; m1 ) m 1 B r U 0 ; P B j=1 Y j;m;U m1 0 ; m1 t N g(X j;m [t 0 ;t N ] ) 2 ; OutputU M 0 as approximation of u(0; x) We shall remark that the key dierence between standard PDEs and PPDEs is the form ofZ dictated by the nonlinear Feynman-Kac, whereZ t =h@ x u(t;X t t ); t;x i is path dependent for PPDEs and decays to only state dependent Z t =h@ x u(t;X t );i for standard PDEs. This explains why we choose the form 21 V n (fX t k g 0kn ; n ) :R dn !R d instead ofV n (X tn ; n ) :R d !R d . Another point that is interesting and worth noticing is, unlike the typical settings with most supervised learning problem, there seems no overtting problem in this framework as we can draw as many independent paths from (X;Y). This suggests one might in theory solve the discretized minimization problem (2.12) up to arbitrary accu- racy, but we should be aware there will always be discretization error involved in the end. 22 Chapter 3 Stochastic Volterra Integral Equations In this chapter, we prove the wellposedness of the SVIEs in the following general form, where the coecients are path-dependent and possibly singular X t;x s = x s + Z s t b(s;r;X t;x ^r )dr + Z s t (s;r;X t;x ^r )dW r : (3.1) the result nicely unied all existing ones with general enough conditions on the coecients. In addition to a weaker version under Assumptions 2.1.1, we also provide a more general wellposedness theorem under weaker but complicated assumptions, which covers more general singular type kernel, e.g. (ts) H 1 2 (X ) for some Lipschitz . The main diculty for the proof is we need to estimate the path norm but the classical BDG inequality cannot be applied to SVIEs, as they are not semi-martingale. To overcome it, we shall use a embedding inequality lemma, which will play the role of BGD inequality in our proof. 23 Once the wellposedness is established, we proceed to prove the regularity and ow property of the solution to SVIEs, which extends all the classical results in SDE. The ow property is crucial in the design of our numerical scheme and the regularity results will be needed in the proof of convergence. The regularity results for regular or special form b and were also studied in the literature, see [4, 9], but the ow property is quite novel based on the work [45]. 3.1 Wellposedness We start with the lemma that bounds the sup norm of a function. Lemma 3.1.1. Let (X;jj) be a Banach space and : [0;T ]!X be -H older continuous for some > 0, then for any p> 1 sup s 0 ;s2[0;T ] j(s 0 )(s)j js 0 sj C Z T 0 Z T 0 j(u)(v)j p juvj p+2 dudv 1 p sup s2[0;T ] j(s)jCT Z T 0 Z T 0 j(u)(v)j p juvj p+2 dudv 1 p +T 1 p Z T 0 j(u)j p du 1 p : 24 Proof. For any 0s<s 0 T and 0<"<s 0 s, j(s 0 )(s)j 1 " Z s 0 s 0 " j(s 0 )(u)jdu + 1 " 2 Z s 0 s 0 " Z s+" s j(u)(v)jdudv + 1 " Z s+" s j(v)(s)jdv 1 " 2 Z s 0 s 0 " Z s+" s j(u)(v)j p juvj p+2 dudv 1 p Z s 0 s 0 " Z s+" s juvj p p1 + 2 p1 dudv 1 1 p +C" sup s 0 ;s2[0;T ] j(s 0 )(s)j js 0 sj C" h Z s 0 s 0 " Z s+" s j(u)(v)j p juvj p+2 dudv 1 p + sup s 0 ;s2[0;T ] j(s 0 )(s)j js 0 sj i Take " = (s 0 s) with C 1 2 , we can derive the rst inequality sup s 0 ;s2[0;T ] j(s 0 )(s)j js 0 sj C Z T 0 Z T 0 j(u)(v)j p juvj p+2 dudv 1 p : Since for any 0sT , j(s)j 1 T Z T 0 j(s)(u)jdu + 1 T Z T 0 j(u)jdu CT sup s 0 ;s2[0;T ] j(s 0 )(s)j js 0 sj +T 1 p Z T 0 j(u)j p du 1 p 25 combining with the rst inequality, we prove the second one sup s2[0;T ] j(s)jCT sup s 0 ;s2[0;T ] j(s 0 )(s)j js 0 sj +T 1 p Z T 0 j(u)j p du 1 p CT Z T 0 Z T 0 j(u)(v)j p juvj p+2 dudv 1 p +T 1 p Z T 0 j(u)j p du 1 p The following proposition captures a main argument we will be using repeatedly in the proof of wellposedness, thus we formalize it to reduce redundancy. Proposition 3.1.2. Let Assumption 2.1.1 holds and V s be any continuous process such that EjV j p <1 for all p> 0, further assume (;)(s;t;!) : [0;T ] [0;T ] !R d R dn are progressively measurable for each xed s and j s;t jK 1 (s;t)jV ^t j; j s;t jK 2 (s;t)jV ^t j j s 0 ;t s;t j Z s 0 s K 1 (u;t)jV ^t jdu; j s 0 ;t s;t j Z s 0 s K 2 (u;t)jV ^t jdu Then the following process U s = Z s 0 s;r dr + Z s 0 s;r dW r 26 is -H older continuous a.s. and for any independent 2 [0;T ], a.s. EjU ^ j p C p;T I p p E Z 0 jV ^r j p dr where = 1 4 , p 8 1 and I p p dened as (2.2). Proof. For any p 8 1 , 0ss 0 T and = 1 4 , EjU s 0U s j p C p E h Z s 0 s s 0 ;r dr p + Z s 0 s s 0 ;r dW r p + Z s 0 s 0 ;r s;r dr p + Z s 0 s 0 ;r s 0 ;r dW r p i C p E h Z s 0 s K 1 (s 0 ;r)jV ^r jdr p + Z s 0 s K 2 2 (s 0 ;r)jV ^r j 2 dr p 2 + Z s 0 Z s 0 s K 1 (u;r)jV ^r jdudr p + Z s 0 Z s 0 s K 2 (u;r)jV ^r jdu 2 dr p 2 i C p Z s 0 0 EjV ^r j p dr n Z s 0 s K 1 (s 0 ;r)dr p js 0 sj p(1) 1 + Z s 0 s K 2 2 (s 0 ;r)dr p 2 js 0 sj p(1) 2 1 + sup 0us 0 Z u 0 K p p1 1 (u;r)dr p1 js 0 sj p + sup 0us 0 Z u 0 K 2p p2 2 (u;r)dr p2 2 js 0 sj p o C p;T I p p Z s 0 0 EjV ^r j p drjs 0 sj 3 2 p 27 by Kolmogorov- Centsov Theorem,U s has a -H older continuous modication. For the L p -norm, EjU s j p C p E h Z s 0 s;r dr p + Z s 0 s;r dW r p i C p E h Z s 0 K 1 (s;r)jV ^r jdr p + Z s 0 K 2 2 (s;r)jV ^r j 2 dr p 2 i C p;T I p p Z s 0 EjV ^r j p dr therefore by Lemma 3.1.1, EjU ^t j p C p;T I p p E[ Z t 0 jV ^r j p dr] Since is independent, EjU ^ j p =E E[jU ^t j p ] t= C p;T I p p E E[ Z t 0 jV ^r j p dr] t= =C p;T I p p E Z 0 jV ^r j p dr: Now we are ready to present the rst wellposedness theorem for SVIEs under stronger assumptions. 28 Theorem 3.1.3. Let Assumption 2.1.1 holds, then for any t2 [0;T ] and adapted continuous process on [0;T ], the following Volterra SDE X t; s = s + Z s t b(s;r;X t; ^r )dr + Z s t (s;r;X t; ^r )dW r (3.2) has unique strong continuous solution X t; on [t;T ] such that EjX t; j p C p;T (Ej j p +I p p ): for p 8 1 and I p p dened in (2.2). Moreover if is -H older continuous for any 2 (0; 1 4 ] a.s., then so is X t; . Proof. Consider the Picard iteration for the forward Volterra SDE: X t;;0 s = s ; X t;;n s = s + Z s t b(s;r;X t;;n1 ^r )dr + Z s t (s;r;X t;;n1 ^r )dW r (3.3) to see the iteration is well-dened, we rst rewrite the X t;;n s as follows: X t;;n s = s + Z s t b(s;r; 0)dr + Z s t (s;r; 0)dW r + Z s t b(s;r;X t;;n1 ^r )b(s;r; 0)dr + Z s t (s;r;X t;;n1 ^r )(s;r; 0)dW r =:U t; s + Z s t s;r dr + Z s t s;r dW r 29 where, under the Assumptions 2.1.1, we can easily check (b(s;r; 0);(s;r; 0)) and (;) satisfy the conditions in Proposition 3.1.2 respectively with V = 1 and V =X t;;n1 , thus X t;;n is -H older continuous and EjX t;;n ^s j p C p;T EjU t; j p + Z s t EjX t;;n1 ^r j p dr C p;T Ej j p +I p p + Z s t EjX t;;n1 ^r j p dr This shows X t;;n preserves the regularity and integrability of X t;;n1 , thus the iteration is well-dened. Now we analyze the dierence X t;;n+1 s X t;;n s = Z s t b(s;r;X t;;n ^r )b(s;r;X t;;n1 ^r )dr + Z s t (s;r;X t;;n ^r )(s;r;X t;;n1 ^r )dW r =: Z s t 0 s;r dr + Z s t 0 s;r dW r again we can easily check ( 0 ; 0 ) satises the conditions in Proposition 3.1.2 with V =X t;;n X t;;n1 , thus EjX t;;n+1 ^s X t;;n ^s j p C ;p;T Z s t EjX t;;n ^r X t;;n1 ^r j p dr C n ;p;T (st) n n! (1 +Ej j p ) Therefore X t;x;n converges to some X t; with continuous path a.s. and EjX t; ^s X t;;n ^s j p 1 X m=n C n ;p;T (st) n n! 1 p p (1 +Ej j p )r n ;p;T (1 +Ej j p ) 30 for somer n ;p;T < 1 andn large enough, and easy to verify by the Lipschitz condition thatX t; is indeed the solution to (3.2). Using the estimate we derived above and applying the dierence estimate for U, V 0 in Proposition 3.1.2, we can have: EjX t;;n ^s j p C p;T Ej j p +I p p + Z s t EjX t;;n1 ^r j p dr taking the limit, EjX t; ^s j p C p;T Ej j p +I p p + Z s t EjX t; ^r j p dr by Gr onwall's inequality EjX t; j p C p;T (Ej j p +I p p ): By Proposition 3.1.2, the X t; is -H older continuous with any exponent 2 (0; 1 4 ] a.s., thus if is -H older continuous, so is X t; . For uniqueness, assume there is another strong solution ~ X t; , then X t; s ~ X t; s = Z s t b(s;r;X t; ^r )b(s;r; ~ X t; ^r )dr + Z s t (s;r;X t; ^r )(s;r; ~ X t; ^r )dW r =: Z s t 00 s;r dr + Z s t 00 s;r dW r 31 by the same argument above, we can derive EjX t; ^s ~ X t; ^s j p C ;p;T Z s t EjX t; ^r ~ X t; ^r j p dr C n ;p;T (st) n n! (1 +Ej j p ) for any n> 0. Therefore X t; = ~ X t; a.s. Next we provide a wellposedness result under weaker assumptions. Assumption 3.1.4. Suppose (b;) : R + R + C 0 ([0;T ];R d )!R d R dn are measurable w.r.t the -algebra induced by usual norm, and let K ' (s;t) := sup x;y j'(s;t; x)'(s;t; y)j jx yj +'(s;t; 0); F ' (s 0 ;s;t) := sup x j'(s 0 ;t; x)'(s;t; x)j 1 +jxj G ' (s 0 ;s;t) := sup x;y j('(s 0 ;t; x)'(s 0 ;t; y)) ('(s;t; x)'(s;t; y))j jx yj the following holds sup 0sT Z s 0 K b (s;t) +K 2 (s;t)dt<1 Z s 0 F b (s 0 ;s;t) +F 2 (s 0 ;s;t) +G 2 (s 0 ;s;t)dtC T js 0 sj for some > 1; > 0. 32 Theorem 3.1.5 (Wellposedness under weaker assumptions). Let Assumption 3.1.4 holds, then for any t2 [0;T ] and adapted continuous process on [0;T ], the SVIE (3.2) has unique strong continuous solution X t; on [t;T ] such that EjX t; j p C p;T (1 +Ej j p ): for p 4 1 . If is -H older continuous for any 2 (0; minf 1 4 ; 2 g] a.s., so is X t; . Proof. Consider the Picard iteration for the forward Volterra SDE: X t;;0 s = s ; X t;;n s = s + Z s t b(s;r;X t;;n1 ^r )dr + Z s t (s;r;X t;;n1 ^r )dW r (3.4) 33 to see the iteration is well-dened, we rst show that X t;;n preserves regularity and integrability of X t;;n1 s , for tss 0 and p 4 1 : EjX t;;n s 0 X t;;n s j p C p E h s 0 s p + Z s 0 s b(s 0 ;r;X t;;n1 ^r )dr p + Z s 0 s (s 0 ;r;X t;;n1 ^r )dW r p + Z s t b(s 0 ;r;X t;;n1 ^r )b(s;r;X t;;n1 ^r )dr p + Z s t (s 0 ;r;X t;;n1 ^r )(s;r;X t;;n1 ^r )dW r p i C p E h s 0 s p + Z s 0 s K b (s 0 ;r)(1 +jX t;;n1 ^r j)dr p + Z s 0 s K 2 (s 0 ;r)(1 +jX t;;n1 ^r j) 2 dr p 2 + Z s t F b (s 0 ;s;r)(1 +jX t;;n1 ^r j)dr p + Z s t F 2 (s 0 ;s;r)(1 +jX t;;n1 ^r j) 2 dr p 2 i C p E h j s 0 s j p + 1 + Z s 0 t jX t;;n1 ^r j p dr n Z s 0 s K b (s 0 ;r)dr p js 0 sj p1 p + Z s 0 s K 2 (s 0 ;r)dr p 2 js 0 sj p 2 1 p 2 + Z s t F b (s 0 ;s;r)dr p + Z s t F 2 (s 0 ;s;r)dr p 2 oi C p 1 + Z s 0 t jX t;;n1 ^r j p dr h j s 0 s j p +js 0 sj p1 p +js 0 sj p 2 1 p 2 +js 0 sj p +js 0 sj p 2 i C p 1 + Z s 0 t jX t;;n1 ^r j p dr js 0 sj p 34 where := minf 1 4 ; 2 g. If E[jX t;;n1 j p ] <1 for all p > 0, then X t;;n has a modication that is H older-continuous with any exponent 2 (0; ) a.s., by Kolmogorov- Centsov Theorem. On the other hand, EjX t;;n s j p C p E h j s j p + Z s t b(s;r; 0)dr p + Z s t (s;r; 0)dW r p + Z s t b(s;r;X t;;n1 ^r )b(s;r; 0)dr p + Z s t (s;r;X t;;n1 ^r )(s;r; 0)dW r p i C p E h j s j p + Z s t b(s;r; 0)dr p + Z s t (s;r; 0)dW r p + Z s t K b (s;r)jX t;;n1 ^r jdr p + Z s t K 2 (s;r)jX t;;n1 ^r j 2 dr p 2 i C p;T Ej s j p + Z s t b(s;r; 0)dr p + Z s t (s;r; 0)dW r p + Z s t EjX t;;n1 ^r j p dr and by Lemma 3.1.1, we have for 2 (0; ) EjX t;;n j p C p;T E h sup s 0 ;s2[t;T ] jX t;;n s 0 X t;;n s j p js 0 sj p + 1 Tt Z T t jX t;;n u j p du i C p;T E h Z T t Z T t jX t;;n u X t;;n v j p juvj p+2 dudv + 1 Tt Z T t jX t;;n u j p du i C p;T (1 +Ejj p + Z T t EjX t;;n1 ^r j p dr): Therefore the Picard iteration is well-dened with the corresponding H older reg- ularity and any p-th moment of the sup norm. Now we proceed to prove the 35 convergence of the solution through estimates on drift and volatility term respec- tively. For drift term, we have E sup tus Z u t b(u;r;X t;;n ^r )b(u;r;X t;;n1 ^r ) dr p E sup tus Z u t K p p1 1 (u;r)dr p1 Z u t X t;;n ^r X t;;n1 ^r p dr C p;T Z s t EjX t;;n ^r X t;;n1 ^r j p dr the volatility term is more complicated, let n;n1 (s;r;X t;; ^r ) :=(s;r;X t;;n ^r )(s;r;X t;;n1 ^r ) then by Lemma 3.1.1, E h sup tus Z u t (u;r;X t;;n ^r )(u;r;X t;;n1 ^r )dW r p i E h C p (st) Z s t Z u t R u t n;n1 (u;r;X t;; ^r )dW r R v t n;n1 (v;r;X t;; ^r )dW r p juvj p+2 dvdu + 1 st Z s t Z u t n;n1 (u;r;X t;; ^r )dW r p du i 36 the second term is easier to bound E Z u t n;n1 (u;r;X t;; ^r )dW r p E Z u t K 2 2 (u;r)jX t;;n ^r X t;;n1 ^r j 2 dr p 2 E Z u t K p p=21 2 (u;r)dr p=21 Z u t X t;;n ^r X t;;n1 ^r p dr C p;T Z s t EjX t;;n ^r X t;;n1 ^r j p dr to bound the rst term E Z u t n;n1 (u;r;X t;; ^r )dW r Z v t n;n1 (v;r;X t;; ^r )dW r p C p E h Z u v n;n1 (u;r;X t;; ^r )dW r p + Z v t n;n1 (u;r;X t;; ^r ) n;n1 (v;r;X t;; ^r )dW r p i C p E h Z u v K 2 2 (u;r)jX t;;n ^r X t;;n1 ^r j 2 dr p 2 + Z v t G 2 (u;v;r)jX t;;n ^r X t;;n1 ^r j 2 dr p 2 i C p E h Z u v K 2 2 (u;r)dr p Z u v jX t;;n ^r X t;;n1 ^r j p drjuvj (1 1 2 p ) p 2 + Z v t G p p=21 (u;v;r)dr p=21 Z v t jX t;;n ^r X t;;n1 ^r j p dr i C p;T E h Z u v jX t;;n ^r X t;;n1 ^r j p drjuvj (1 1 2 p ) p 2 + Z v t G 2 dr p 2 Z v t jX t;;n ^r X t;;n1 ^r j p dr i C p;T Z s t EjX t;;n ^r X t;;n1 ^r j p drjuvj p 37 putting them together, we have EjX t;;n+1 ^s X t;;n ^s j p C p E h sup tus Z u t b(u;r;X t;;n ^r )b(u;r;X t;;n1 ^r ) dr p + sup tus Z u t n;n1 (u;r;X t;; ^r )dW r p i C p;T Z s t EjX t;;n ^r X t;;n1 ^r j p dr C n p;T (st) n n! C 0 p;T (1 +Ejj p ) Therefore X t;x;n converges to some X t; with continuous path a.s. and EjX t; ^s X t;;n ^s j p 1 X m=n C n p;T (st) n n! 1 p p (1 +Ej j p )r n p;T (1 +Ej j p ) for somer n p;T < 1 andn large enough, and easy to verify by the Lipschitz condition that X t; is indeed the solution to (3.2). Using the estimate we derived above: EjX t;;n ^s j p C p;T 1 +Ej j p + Z s t EjX t;;n1 ^r j p dr taking the limit, EjX t; ^s j p C p;T 1 +Ej j p + Z s t EjX t; ^r j p dr 38 by Gr onwall's inequality EjX t; j p C p;T (1 +Ej j p ): Since the X t; is -H older continuous with any exponent 2 (0; ] a.s. by the same argument for X t;;n , thus if is -H older continuous, so is X t; . For uniqueness, assume there is another strong solution ~ X t; , then X t; s ~ X t; s = Z s t b(s;r;X t; ^r )b(s;r; ~ X t; ^r )dr + Z s t (s;r;X t; ^r )(s;r; ~ X t; ^r )dW r by the same argument above, we can derive EjX t; ^s ~ X t; ^s j p C p;T Z s t EjX t; ^r ~ X t; ^r j p dr C n p;T (st) n n! (1 +Ej j p ) for any n> 0. Therefore X t; = ~ X t; a.s. Remark 3.1.6. The weaker wellposedness result is aimed for presenting the main idea and easier reference, and the stronger version is intended to cover more general singular kernel, such as (ts) H 1 2 (X ) withH2 (0; 1 2 ) and Lipschitz, which we think are important. For anyF stopping time, we can deneX ; similarly as the solution to (3.2) with integral part starting from , and its well-posedness follows from that for X 0; with b(s;r; x) =1 frg b(s;r; x), (s;r; x) =1 frg (s;r; x). 39 3.2 Regularity and Flow Property We rst state the regularity of the solutions in terms of the initial time and process. Proposition 3.2.1. Let Assumptions 2.1.1 hold, then for any s;t2 [0;T ] and continuous adapted ; with Ejj p <1, Ejj p <1, EjX s; ^u X t; ^u j p C ;p;T Ejj p +Ejj p ^Ejj p jstj : Proof. WLOG assume ts andEjj p Ejj p , rewrite X s; X t; as follows X s; u X t; u = u u + Z u t u;r dr + Z u t u;r dW r where u;r =1 frsg b(u;r;X s; ^r )b(u;r;X t; ^r ) 1 fr<sg b(u;r;X t; ^r ) u;r =1 frsg (u;r;X s; ^r )(u;r;X t; ^r ) 1 fr<sg (u;r;X t; ^r ) it is straightforward to verify that u;r , u;r satisfy the conditions in Proposition 3.1.2 with V r =1 frsg jX s; ^r X t; ^r j +1 fr<sg j1 +X t; ^r j 40 therefore it follows that EjX s; ^u X t; ^u j p C ;p;T Ejj p +Ejj p jstj + Z u s EjX s; ^r X t; ^r j p dr and by Gr onwall's ineqaulity EjX s; ^u X t; ^u j p C ;p;T Ejj p +Ejj p jstj : The second regularity result is for the adjoint concatenated path, which is the crucial piece for ow property and nonlinear Feynman-Kac. Proposition 3.2.2. Let assumptions in the Theorem 3.1.3 hold true, for anyst, dene (X t; s ) u :=X t; u 1 fusg +1 fu>sg u + Z s t b(u;r;X t; ^r )dr + Z s t (u;r;X t; ^r )dW r (3.5) as in (2.5), then X t; s is -H older continuous for any 2 (0; ), and for any 0 2 (0; 1), p max( 8 1 ; 2 1 0 ) the following estimates hold E sup tsT jX t; s j p C ;;T (1 +Ejj p ) E sup s 0 ;s2[t;T ] jX t; s 0 X t; s )j js 0 sj 0 p C ;;T (1 +Ejj p ): 41 Proof. First for any 0 s s 0 T , we can rewrite X t; s and X t; s 0 X t; s as (X t; s ) u = u + Z u t 1 f0rsg b(u;r;X t; ^r )dr + Z u t 1 f0rsg (u;r;X t; ^r )dW r (X t; s 0 X t; s ) u = Z u 0 1 fsrs 0 g b(u;r;X t; ^r )dr + Z u 0 1 fsrs 0 g (u;r;X t; ^r )dW r it is straightforward to verify their coecients satisfy the conditions in Proposition 3.1.2 respectively with V ^t = 1 ftsg (1 +jX t; ^t j) and V ^t = 1 fsts 0 g (1 +jX t; ^t j), thus both of them are -H older continuous under the assumptions and EjX t; s j p C ;T (1 +Ejj p ) EjX t; s 0 X t; s )j p C ;;T (1 +Ejj p )js 0 sj applying Lemma 3.1.1, we can conclude for any 0 2 (0; 1), p max( 8 1 ; 2 1 0 ) E sup tsT jX t; s j p C ;;T (1 +Ejj p ) E sup s 0 ;s2[t;T ] jX t; s 0 X t; s )j js 0 sj 0 p C ;;T (1 +Ejj p ) Finally the ow property is provided. 42 Theorem 3.2.3. Let Assumptions 2.1.1 hold and X t; be the solution to the Volterra SDE (3.2) as proved in Theorem 3.1.3, then X t; has the ow property, i.e. X t; (!) =X ;X t; (!) =X (!);X t; (!) (!) a:s: for anyF stopping time . Proof. First it is straightforward to check that X t; solves the following Volterra SDE X u = (X t; ) u + Z u b(u;r; X ^r )dr + Z u (u;r; X ^r )dW r : by denition,X ;X t; is also the solution, thus by uniquenessX t; =X ;X t; a.s. For the second part of the ow property, we rst consider the case s for some xeds. SinceC([0;T ];R d ) is separable, there exists partition A n;i and x n;i 2A n;i , such that[ 1 i=1 A n;i =C([0;T ];R d ) for alln andjx n;i yj 1 n for any y2A n;i . Let n = 1 X i=1 x n;i 1 fX t; s2A n;i g then n 2F s and lim n!1 EjX t; s n j p lim n!1 1 n 2p = 0 43 Since n are simple functions andfX t; s 2A n;i g2F s , we have X s; n (!) = 1 X i=1 X s;x n;i (!)1 fX t; s2A n;i g =X s; n (!) (!) a:s: Now on one hand, by Proposition 3.1.2 and Gr onwall's ineqaulity EjX s;X t; s X s; n j p C ;p;T EjX t; s n j p : on the other hand, since X s;(!) F s ?F s for any 2F s , EjX s;X t; s(!) X s; n (!) j p =E h E[ X s;X t; s(!) X s; n (!) p F s ] i =E h E[ X s;x X s;y p ] x=X t; s;y= n i C ;p;T EjX t; s n j p therefore lim n!1 EjX s;X t; s X s; n j 2p = 0: This implies X s;X t; s(!) (!) =X s;X t; s (!) a.s., thus X (!);X t; (!) (!) =X ;X t; (!) =X t; a:s: 44 for discrete stopping time . For general stopping time , we can approximate it by discrete one n such that 0 n 1 n , then E X (!);X t; (!) X n(!);X t; n (!) p = Nn X k=1 E h X (!);X t; (!) X s k ;X t; s k (!) p 1 fs k1 <s k g i = Nn X k=1 E h E X s k ;X (!);X t; (!) s k X s k ;X t; s k (!) p F s k 1 fs k1 <s k g i C ;p;T Nn X k=1 E h X (!);X t; (!) s k X t; s k p 1 fs k1 <s k g i C ;p;T Nn X k=1 E h X (!);X t; (!) s k X t; p + X t; X t; s k p 1 fs k1 <s k g i C ;;p;T Nn X k=1 E h E X (!);X t; (!) s k X t; p F + (1 +Ejj p ) s k p 2 1 fs k1 <s k g i C ;;p;T Nn X k=1 E h E X s;x s k x p (s;x)=(;X t; ) + (1 +Ejj p ) s k p 2 1 fs k1 <s k g i C ;;p;T (1 +Ejj p ) n p 2 45 the last three inequalities follows from Proposition 3.2.1 and (3.2.2). Therefore we can conclude X (!);X t; (!) = lim n!1 X n(!);X t; n (!) =X t; a:s: Remark 3.2.4. We would like to remark that the recovery of the ow property relies crucially on the nontrivial concatenated processX t; , and the statement of ow property are often deceptively easy compared to the actual meaning and proof. 46 Chapter 4 Branching Diusion Method Built upon the wellposedness and basic properties derived for SVIEs in the previous chapter, our focus in this part is to rigorously generalize the Branching diusion method to SVIE associated PPDEs. The main issue in adapting the numerical scheme is to be able to simulate the path-dependent X in a Markovian manner, and it is dealt with by the ow property Theorem 3.2.3 and carrying the concate- nated process X t; s through the simulation of child processes. Corresponding representation results are also established. After that, we discuss the blow-up issues for the Branching method. As alluded in the introduction, the main observation here is theL p norm of the representation corresponds to a transformed PPDEs or BSDEs with smaller time duration, which is why a completely well-posed equation can blow up rather quickly in Monte Carlo simulations. The solution we propose is to transform the original equation in a way that not only maintain the polynomial structure but also have positive polynomial coecients. This allows the Branching method to be still applicable to the new equation but will not have the moment degrading problem. Numerical examples conrm this and could be found in both current and later chapters. We show 47 the existence of such transformations for general quadratic case, but whether such transformations exist for higher order polynomials or to what extent they exist could be a interesting question for future research. 4.1 Branching Processes Setup To start with, we shall rst dene the branching process rigorously. Recall that our probability space ( ;F;P) is assumed to be rich enough for any countable inclusion of random variables independent of W , here are the extra ones needed for the setup a sequence of i.i.d non-negative random variables k with density function, a sequence of i.i.d non-negative integer-valued random variables I k with P(I k =`) =p ` ; a sequence of i.i.d Brownian motions W k on [0;T ]. where k2[ 1 n=1 N n , and (W ; k ;I k ; W k ) are all independent from each other. For any (t; x)2 [0;T ]C([0;T ];R d ), we recursively dened the branching process (T k ;I k ;X t;x;k ) as follows: 48 Start with index k = (1),K 0 :=;,K 1 :=f(1)g T (1) := (t + (1) )^T; W (1) s := W (1) s W (1) t ; X t;x;(1) s = x s + Z s^T (1) t b(s;r;X t;x;(1) ^r )dr + Z s^T (1) t (s;r;X t;x;(1) ^r )dW (1) r ; s2 [t;T ]: DeneK n+1 := k (`) 1`I k ;T k <T; k2K n , and for all k (`)2 K n+1 , T k(`) := (T k + k(`) )^T; W k(`) s :=W k s^T k + W k(`) s_T k W k(`) T k ; X t;x;k(`) s =X t;x;k s +1 fs>T k g Z s^T k(`) T k b(s;r;X t;x;k(`) ^r )dr + Z s^T k(`) T k (s;r;X t;x;k(`) ^r )dW k(`) r : To simplify the nal representation results, we also introduce notations T (1) :=t; T k(`) :=T k ; T k :=T k T k ; K n T := k T k =T; k2K n ; K 1 T :=[ n1 K n T ; K 1 :=[ n1 K n : Note the T k dened here should in fact be T t k , since it depends on t as well, but we will suppress it to make the notation cleaner in most cases and only denote them when it is important to emphasize such dependency. 49 The following two propositions are to ensure the desired branching process is well- dened. Proposition 4.1.1. If E[I k ] <1, then the number of particles inK 1 is nite a.s. Proof. See e.g. Theorem 1 of Athreya-Ney [1], or Harris [29]. Proposition 4.1.2. Under Assumptions 2.1.1, for any k2K 1 and 1`I k , (i) W k is a Brownian Motion on [t;T ] and W k(`) ^T k =W k ^T k , (ii) There exists such X t;x;k following the recursion, moreover X t;x;k(`) ^T k = X t;x;k ^T k and X t;x;k s = x s + Z s^T k t b(s;r;X t;x;k ^r )dr + Z s^T k t (s;r;X t;x;k ^r )dW k r : Proof. (i) Prove by induction, and show the nite dimensional distribution is the same as the Brownian motion by conditioning on T k and using the independence of T k and W k , the continuity of the path follows by denition. (ii) The existence follows from the well-posedness of the Volterra SDE: X t;x;k s = x s + Z s t b(s;r; X t;x;k ^r )dr + Z s t (s;r; X t;x;k ^r )dW k r : and by Proposition 3.1.2, X t;x;k(`) ^T k = X t;x;k ^T k . We can let X t;x;k s := x s + Z s^T k t b(s;r; X t;x;k ^r )dr + Z s^T k t (s;r; X t;x;k ^r )dW k r : 50 then the recursion relation and the property in the conclusions follow directly. 4.2 Representation Results With the branching process dened, we can now prove the representation results. Since the arguments are mainly probabilistic, we also choose the equivalent FBSDE formulation in the following form, X t;x s = x s + Z s t b(s;r;X t;x ^r )dr + Z s t (s;r;X t;x ^r )dW r (4.1) Y t;x s =g(X t;x ) + Z T s f(r;X t;x ^r ;Y t;x r )dr Z T s Z t;x r dW r : (4.2) where the nonlinearity f is polynomial on y, i.e. f(t; x;y) = L X `=1 c ` (t; x)y d ` : Then we introduce the branching representation and its nite approximation, which shall be used in the proof t;x := h Y k2K 1 T g(X t;x;k ) F (T k ) ih Y k2K 1 nK 1 T c I k (T k ;X t;x;k ^T k ) p I k (T k ) i ; (4.3) t;x n := h Y k2[ n j=1 K j T g(X t;x;k ) F (T k ) ih Y k2[ n j=1 K j nK j T c I k (T k ;X t;x;k ^T k ) p I k (T k ) ih Y k2K n nK n T Y t;x T k d I k i : 51 It is worth pointing out that, by denitionK n are random sets taking values depending on correspondingT k , thus both t;x and t;x n are in fact sum of dierent product combinations in the partition of , e.g. t;x 1 = g(X t;x;(1) ) F (T (1) ) 1 fT (1) =Tg +1 fT (1) <Tg c I (1) (T (1) ;X t;x;(1) ^T (1) ) p I (1) (T (1) ) (Y t;x T (1) ) d I (1) : Now we can provide the rst representation result, which says if solution to (4.1) exists, then it is represented by (4.3). Assumption 4.2.1. (i) E[I k ] = P ` `p ` <1, (ii) continuous, strictly positive on [0;T ] and F (T ) := R 1 T (t)dt> 0. Theorem 4.2.2. Let Assumption 4.2.1 holds. Assume for some (t; x)2 [0;T ] C([0;T ];R d ), the FBSDE (4.1) has a unique solution (X t;x ;Y t;x ;Z t;x ) and t;x n is uniformly integrable, then t;x 2L 1 and Y t;x t =E[ t;x ]: Proof. First apply expectation to both sides of (4.1) and use the trick of repre- senting the integral by evaluation at an independent stopping time Y t;x t =E h g(X t;x ) + Z T t L X `=1 c ` (r;X t;x ^r ) Y t;x r d ` dr i =E h g(X t;x;(1) ) F (T (1) ) 1 fT (1) =Tg +1 fT (1) <Tg c I (1) (T (1) ;X t;x;(1) ^T (1) ) p I (1) (T (1) ) (Y t;x T (1) ) d I (1) i =E[ t;x 1 ]; 52 Then apply the crucial technique of linearizing the power ofY t;x through indepen- dent samples: (Y t;x s ) d = E h g(X t;x ) + Z T s f(r;X t;x ^r ;Y t;x r )dr F s i d = E h g(X s;X t;x s ) + Z T s f(r;X s;X t;x s ^r ;Y s;X t;x s r )dr F s i d = d Y i=1 E h g(X s;X t;x s;i ) + Z T s f(r;X s;X t;x s;i ^r ;Y s;X t;x s;i r )dr F s i =E h d Y i=1 g(X s;X t;x s;i ) + Z T s f(r;X s;X t;x s;i ^r ;Y s;X t;x s;i r )dr F s i : where X s;X t;x s;i are independent solutions to the forward Volterra SDE (4.1). The second and the third equality follows from ow property in Proposition 3.2.2. 53 The last equality follows from independent property. Thus by tower property and twice of change of measure trick: Y t;x t =E h g(X t;x ) + Z T t L X `=1 c ` (r;X t;x ^r ) Y t;x r d ` dr i =E h g(X t;x ) + Z T t L X `=1 c ` (r;X t;x ^r )E h d ` Y i=1 g(X s;X t;x s;(i) ) + Z T s f(u;X s;X t;x s;(i) ^u ;Y s;X t;x s;(i) u )du F s i dr i =E h g(X t;x ) + Z T t L X `=1 c ` (r;X t;x ^r ) d ` Y i=1 g(X s;X t;x s;(i) ) + Z T s f(u;X s;X t;x s;(i) ^u ;Y s;X t;x s;(i) u )du dr i =E h g(X t;x;(1) ) F (T (1) ) 1 fT (1) =Tg +1 fT (1) <Tg c I (1) (T (1) ;X t;x;(1) ^T (1) ) p I (1) (T (1) ) d I (1) Y j=1 g(X t;x;(1;j) ) F (T (1;j) ) 1 fT (1;j) =Tg +1 fT (1;j) <Tg c I (1;j) (T (1;j) ;X t;x;(1;j) ^T (1;j) ) p I (1;j) (T (1;j) ) (Y t;x T (1;j) ) d I (1;j) i =E[ t;x 2 ] iterating the procedure, it can be prove by induction that Y t;x t =E[ t;x n ] for all n 1; since t;x n is uniformly integrable by the assumption, we can conclude Y t;x t = E[ t;x ]. 54 The next theorem concerns the other direction, which is the representation in (4.3) gives one solution to the BSDE (4.1). Theorem 4.2.3. Let Assumptions 2.1.1 and 4.2.1 holds. Assume for all (t; x)2 [0;T ]C([0;T ];R d ), t;x 2L 1 and E h Z T t v(r;X t;x r ) + L X `=1 jc ` j(r;X t;x ^r )v d ` (r;X t;x r )dr i <1 where v(t;x) :=Ej t;x j, then there exists a solution (Y t;x ;Z t;x ) to the BSDE (4.1) with Y t;x s =u(s;X t;x s ); u(t; x) :=E[ t;x ] 55 Proof. First by Theorem 3.1.3, X t;x and t;x are well-dened. By the denition of t;x t;x = g(X t;x;(1) ) F (T (1) ) 1 fT (1) =Tg +1 fT (1) <Tg c I (1) (T (1) ;X t;x;(1) ^T (1) ) p I (1) (T (1) ) d I (1) Y j=1 h Y k2K 1 T nK 1 ;k 2 =j g(X t;x;k ) F (T T (1) k ) ih Y k2K 1 n(K 1 T [K 1 );k 2 =j c I k (T k ;X t;x;k ^T k ) p I k (T T (1) k ) i = g(X t;x;(1) ) F (T (1) ) 1 fT (1) =Tg +1 fT (1) <Tg c I (1) (T (1) ;X t;x;(1) ^T (1) ) p I (1) (T (1) ) d I (1) Y j=1 h Y k2K 1 T nK 1 ;k 2 =j g(X T (1) (!);X t;x;(1) (!);k ) F (T T (1) k ) i h Y k2K 1 n(K 1 T [K 1 );k 2 =j c I k (T T (1) k ;X T (1) (!);X t;x;(1) (!);k ^T k ) p I k (T T (1) k ) i =: g(X t;x;(1) ) F (T (1) ) 1 fT (1) =Tg +1 fT (1) <Tg c I (1) (T (1) ;X t;x;(1) ^T (1) ) p I (1) (T (1) ) d I (1) Y j=1 T (1) (!);X t;x;(1) (!);j 56 where the second equality follows from the ow property and T k being indepen- dent ofW k , the last equality is a pure denition of T (1) (!);X t;x;(1) (!);j as the expres- sion in the product to simplify notation. Using the independent properpty of (I k ;T k ;X t;x;k ) E d I (1) Y j=1 T (1) (!);X t;x;(1) (!);j T (1) ;I (1) ;W (1) = d I (1) Y j=1 E T (1) (!);X t;x;(1) (!);j T (1) ;I (1) ;W (1) = d I (1) Y j=1 E[ s;y ] (s;y)=(T (1) ;X t;x;(1) ) =u d I (1) (T (1) ;X t;x;(1) ) the integrability of u(T (1) ;X t;x;(1) ) is guaranteed by the condition stated in the theorem and is bounded below on [0;T ], thus the conditional expectation is valid. Substituting it in the above derivation, we derive u(t; x) =E[ t;x ] =E h g(X t;x;(1) ) F (T (1) ) 1 fT (1) =Tg +1 fT (1) <Tg c I (1) (T (1) ;X t;x;(1) ^T (1) ) p I (1) (T (1) ) u d I (1) (T (1) ;X t;x;(1) ) i =E h g(X t;x ) + Z T t L X `=1 c ` (r;X t;x ^r )u d ` (r;X t;x r )dr i : 57 Now for anyF W stopping time , we can combine the ow property and similar tricks E[u(;X t;x )] =E[E[ (!);X t;x (!) jF ]] =E[ (!);X t;x (!) ] =E h g(X (!);X t;x (!);(1) ) F (T (1) ) 1 fT (1) =Tg +1 fT (1) <Tg c I (1) (T (1) ;X (!);X t;x (!);(1) ^T (1) ) p I (1) (T (1) ) u d I (1) (T (1) ;X (!);X t;x (!);(1) ) i =E h g(X t;x;(1) ) F (T (1) ) 1 fT (1) =Tg +1 fT (1) <Tg c I (1) (T (1) ;X t;x;(1) ^T (1) ) p I (1) (T (1) ) u d I (1) (T (1) ;X t;x;(1) ) i =E h g(X t;x ) + Z T L X `=1 c ` (r;X t;x ^r )u d ` (r;X t;x r )dr i : the third equality use the same trick as above to simplify the gigantic product into u d I (1) (T (1) ;X (!);X t;x (!);(1) ), by conditioning on (T (1) ;I (1) ;W (1) ;F W ), the forth equality is implied by the ow property and the integrability of (!);X t;x (!) can be seen from Ej (!);X t;x (!) j =E h jgj(X t;x;(1) ) F (T (1) ) 1 fT (1) =Tg +1 fT (1) <Tg jcj I (1) (T (1) ;X t;x;(1) ^T (1) ) p I (1) (T (1) ) v d I (1) (T (1) ;X t;x;(1) ) i E h jgj(X t;x ) + Z T t L X `=1 jc ` j(r;X t;x ^r )v d ` (r;X t;x r )dr i <1 58 Combining with the expression for u(t;x), we can derive u(t; x) =E h u(;X t;x ) + Z T L X `=1 c ` (r;X t;x ^r )u d ` (r;X t;x r )dr i therefore by martingale representation, we can conclude there exists a solution (Y t;x ;Z t;x ) to BSDE (4.1) with Y t;x s :=u(s;X t;x s ). 4.3 Blow-up Issues and Transformations Our attention to the blow-up issue was drawn in a talk for [31] and we start o by removing the nonessential features and considering the toy equation: @ t u(t;x) + 1 2 u(t;x)u 2 (t;x) = 0; u(T;x) = 1: (4.4) or its equivalent BSDE Y t = 1 Z T t Y 2 r dr Z T t Z r dW r note that this \PDE" admits a solution u(t;x) = 1=(1 +Tt), for any T . How- ever, our numerical simulation of (4.4) using the branching method exhibits large 59 variance once T gets very close to 1 and no longer seems convergent for T 1. If we look at the branching representation: u(0;x) = E h 1 fT (1) =Tg 1 F (T ) 1 fT (1) <Tg u 2 (T (1) ;W 0;x T (1) ) (T (1) ) i =E[ 0;x ] 0;x := h Y k2K 1 T 1 F (T k ) ih Y k2K 1 nK 1 T (1) (T k ) i it can then be realized thatj 0;x j is actually the branching representation for a close but dierent equation: @ t u(t;x) + 1 2 u(t;x) +u 2 (t;x) = 0; u(T;x) = 1: (4.5) with coecients of theu 2 being +1 instead of1, and the root cause of the blow-up is that this new equation is only well-posed up to T = 1. Surprisingly if we dene v =e 2(Tt) (1u), with u being the solution to (4.4), then v will solve the following equation: @ t v(t;x) + 1 2 v(t;x) +e 2(Tt) v 2 (t;x) +e 2(Tt) = 0; v(T;x) = 0: (4.6) the essential dierence between (4.4) and (4.6) is that the branching representation for the latter has the propertyj 0;x j = 0;x , because of the non-negativeness of the terminal condition and coecients in polynomial drifts. Our numerical simulations show that (4.6) can indeed be solved for T = 6 or even larger with the branching 60 method, once thev is solved,u can be solved by simple transformation as well for an extended time interval. The next two results formalize our observations and the transformation ideas. Proposition 4.3.1. Let Assumptions 2.1.1 and 4.2.1 hold, assume further that (t) =e t and t;x 2L q for some q 0 such that E h Z T t v(r;X t;x r ) +q 1q L X `=1 p 1q ` jc ` j(r;X t;x ^r )v d ` (r;X t;x r )dr i <1 where v(t;x) := Ej t;x j q , then there exists a solution (Y t;x ;Z t;x ) with Y t;x s = v(s;X t;x s ) to the following BSDE: Y t;x s =jgj q (X t;x ) +q 1q Z T s L X `=1 p 1q ` jc ` j(r;X t;x ^r )(Y t;x r ) d ` dr Z T s Z t;x r dW r : (4.7) Proof. By the Theorem 4.2.3, it suces to show thatj t;x j q can be written as the branching representation of BSDE (4.7), and this is true because of the special form the stopping density (t). To see it, by denition j t;x j q = h Y k2K 1 T jgj q (X t;x;k ) e qT k ih Y k2K 1 nK 1 T jcj q I k (T k ;X t;x;k ^T k ) p q I k q e qT k i = h Y k2K 1 T jgj q (X t;x;k ) e qT k ih Y k2K 1 nK 1 T q 1q p 1q I k jcj q I k (T k ;X t;x;k ^T k ) p I k qe qT k i = h Y k2K 1 T jgj q (X t;x;k ) F q (T k ) ih Y k2K 1 nK 1 T q 1q p 1q I k jcj q I k (T k ;X t;x;k ^T k ) p I k q (T k ) i 61 where q (t) =qe qt and F q (t) =e qt = R 1 t q (s)ds. Thereforej t;x j q is indeed a branching representation of BSDE (4.7) with dierent stopping density, and the conclusion follows from Theorem 4.2.3. Remark 4.3.2. Combined with Theorem 4.2.3, the above theorem explains why the Branching method typically has a small maturity time or requires small coecients and terminal conditions to work. Because theL 1 andL 2 norm of the representation are necessary for the Monte Carlo simulation to converge and have nite variance, but they in fact correspond to the branching representation of BSDEs (4.7), which tends to blow up quickly as all coecients and terminal taken absolute value. The next Proposition says the transformation could always be found for quadratic f. Proposition 4.3.3. Consider the BSDE in (4.1), assumef(t; x;y) =ay 2 +by +c, with (a;b;c) : [0;T ]C([0;T ];R d )! R RR and g : C([0;T ];R d )! R all being bounded. Then there exists another set of ( a; b; c; g) 0, and deterministic functions m(s);n(s) such that Y t;x s =m(s)+n(s) Y t;x s ; where Y t;x s = g(X t;x )+ Z T s a(Y t;x r ) 2 + bY t;x r + cdr Z T s Z t;x r dW r : 62 Proof. Let Y t;x s := e (Ts) (h(s)Y t;x s ) for some function h(s) and constant , which will be later deduced, then by It^ o formula d Y t;x s = Y t;x s ds +e (Ts) h 0 (s)dsdY t;x s = e (Ts) h 0 (s) Y t;x s dse (Ts) (a(Y t;x s ) 2 bY t;x s c)ds + Z t;x s dW s = ae (Ts) ( Y t;x s ) 2 ( + 2ah(s) +b) Y t;x s +e (Ts) (h 0 (s) +ah 2 (s) +bh(s) +c) ds e (Ts) Z t;x s dW s so it suces to nd and h(s) such that 8 > > > > > > < > > > > > > : + 2ah(s) +b 0 h 0 (s) +ah 2 (s) +bh(s) +c 0 h(T )g denote the upper bound for (a;b;c) as , if we let h(s) = (jgj + 1)e (Ts) 1; = 2(jgj + 1)e T + 63 then it is straightforward to verify all the three conditions are satised and by simply setting a; b; c = ae (Ts) ; + 2ah(s) +b;e (Ts) (h 0 (s) +ah 2 (s) +bh(s) +c) m(s);n(s) = (jgj + 1)e (Ts) 1;e (Ts) we are done with the proof. Remark 4.3.4. In light of Proposition 4.3.1 and 4.3.3, one possible approach to extending the time duration is to nd transformations, which ip the coecients and terminal to be positive, solve the new BSDE with branching representation and then transform the solution back. Since the L 1 norm of the new representation is itself and is guaranteed to be nite as long as the BSDE itself is well-posed, it successfully avoid some of blow-up issues. 64 Chapter 5 Multilevel Picard Approximation In this chapter, we extend the multilevel Picard approximation to PPDEs driven by SVIEs. Compared to the Markovian case, the main complication is the forward simulation of X, as it has to be carried out in a Markovian manner to keep the nested expectation valid. Again the ow property from Theorem 3.2.3 allows us to do so for SVIEs with the concatenated processX t; s , which amounts to tracking the history path ofX. In addition, we introduce the exponential stopping strategy to control the variance through extra parameters. A slightly ner convergence rate is obtained for our scheme. Apart from the convergence results, we also estimate the computation cost and and study the its relation to the parameter choice. Using a argument similar to the structural risk minimization in statistical learning, we derive the range of the optimal for the exponential stopping density. Numerical examples on this issue are also given in the last chapter. 65 5.1 Convergence of Numerical Scheme The rst Proposition is the classical result about the convergence rate for the Picard iteration, which shall be used in the proof of the main theorem. Proposition 5.1.1. Let assumption 2.1.1 and 2.1.2 holds andX t;x be the solution to the Volterra SDE (3.2), consider the following BSDE and its classical approxi- mation scheme Y t;x s =g(X t;x ) + Z T s f(r;X t;x ^r ;Y t;x r )dr Z T s Z t;x r dW r : (5.1) Y t;x;0 s = 0; Y t;x;n s =g(X t;x ) + Z T s f(r;X t;x ^r ;Y t;x;n1 r )dr Z T s Z t;x;n r dW r : (5.2) Then for n 2, the approximation error can be bounded as follows: jY t;x;n t Y t;x;n1 t j 2 C T (Tt) 2n2 (1 +jxj 2 ) 2 n2 (n 2)! ; jY t;x;n t Y t;x t j 2 C T (Tt) 2n+2 (1 +jxj 2 ) 2 n1 (n 1)! Proof. First note that under the assumptions, the well-posedness of the BSDE follows from the classical results. We proceed as the classical arguments: Y t;x;n s Y t;x;n1 s + Z T s Z t;x;n r Z t;x;n1 r dW r = Z T s f(r;X t;x ^r ;Y t;x;n1 r )f(r;X t;x ^r ;Y t;x;n2 r )dr 66 Using the independent increment property of Brownian motion and It^ o isometry, E h (Y t;x;n s Y t;x;n1 s ) 2 + Z T s (Z t;x;n r Z t;x;n1 r ) 2 dr i =E Z T s f(r;X t;x ^r ;Y t;x;n1 r )f(r;X t;x ^r ;Y t;x;n2 r )dr 2 By the uniform Lipschitz condition of f and H older inequality E(Y t;x;n s Y t;x;n1 s ) 2 L 2 (Ts) Z T s E(Y t;x;n1 r Y t;x;n2 r ) 2 dr Repeating the above inequality, we have the following error bound for consecutive approximation E(Y t;x;n s Y t;x;n1 s ) 2 (Ts) 2n2 L 2n2 (2n 3)!! sup suT E h g(X t;x ) + Z T u f(r;X t;x ^r ; 0)dr i 2 2(Ts) 2n2 L 2n2 (2n 3)!! E h g 2 (X t;x ) + (Ts) Z T s f 2 (r;X t;x ^r ; 0)dr i C T (Ts) 2n2 (1 +jxj 2 ) 2 n2 (n 2)! ReplaceY t;x;n1 s by the true solution Y t;x s and repeat the same arguments, we can derive the following error bound E(Y t;x;n s Y t;x s ) 2 2(Ts) 2n+1 L 2n (2n 1)!! E Z T s (Y t;x r ) 2 dr C T (Ts) 2n+2 (1 +jxj 2 ) 2 n1 (n 1)! 67 Now consider the multilevel Picard approximation for the BSDE (5.1): U 0 (t; x) = 0; U n+1 (t; x) = 1 M n+1 M n+1 X i=1 g(X t;x;0;i ) + 1 f 0;i t <Tg ( 0;i t t) F 0 ( 0;i t ;X t;x;0;i ^ 0;i t ;X t;x;0;i 0;i t ) + n X `=1 1 M n+1` M n+1` X i=1 1 f `;i t <Tg ( `;i t t) h (F ` F `1 )( `;i t ;X t;x;`;i ^ `;i t ;X t;x;`;i `;i t ) i : (5.3) where F ` (t; x; x) := f(t; x;U ` (t; x)), X t;x;`;i are independent simulations of solu- tions to the SVIE (3.2),X t;x s is the concatenated process dened in (3.5) and `;i t are i.i.d random variables such that `;i t t has probability density function (s) =e s ;s 0. The idea of the multilevel scheme above is to break the nonlinear Feynman-Kac representation into telescope sum of Picard iteration, i.e. u n+1 (t; x) =E g(X t;x ) + Z T t f(r;X t;x ;u n (r;X t;x ))dr =E g(X t;x ) + Z T t f(r;X t;x ;u 0 (r;X t;x ))dr + n X k=1 E Z T t f(r;X t;x ;u k (r;X t;x ))f(r;X t;x ;u k1 (r;X t;x ))dr and assign simulation resources roughly proportional to the variance of each term, i.e. smaller k gets more simulation while larger ones get less. In this way, the 68 scheme takes better advantages of the fast convergence of Picard iteration, which could be seen in our convergence theorem below. Theorem 5.1.2. Let assumption 2.1.1 and 2.1.2 holds, then the approximation error has the following bound for n 2 E(U n (t; x)Y t;x t ) 2 C T (1 +jxj 2 )e (Tt) 2 M n1 n1 X k=0 max( M ;MT + 1 ) k C (Tt) k k! : (5.4) where C is the C ;;T in Proposition 3.2.2. Proof. Let V n (t; x) := Var(U n (t; x)); B n (t; x) := E[U n (t; x)] Y t;x;n t 2 ; ~ E n (t; x) :=V n (t; x) +B n (t; x), then we have: V n+1 (t; x) = 1 M n+1 Var h g(X t;x ) + 1 f<Tg e (tt) F 0 ( t ;X t;x ^t ;X t;x t ) i + n X `=1 1 M n+1` Var h 1 f<Tg e (tt) F ` F `1 ( t ;X t;x ^t ;X t;x t ) i 69 Using Var[X]E[X 2 ], we bound the recurrent term as Var h 1 f<Tg e (tt) F ` F `1 ( t ;X t;x ^t ;X t;x t ) i E h 1 f<Tg e 2(tt) 2 F ` F `1 2 ( t ;X t;x ^t ;X t;x t ) i =E h Z T t e (rt) F ` F `1 2 (r;X t;x ^r ;X t;x r )dr i 3 E h Z T t e (rt) U ` (r;X t;x r )Y t;x;` r 2 + U `1 (r;X t;x r )Y t;x;`1 r 2 + Y t;x;` r Y t;x;`1 r 2 dr i Therefore we get the following recursion for V n (t; x): V n+1 (t; x) 1 M n+1 Var h g(X t;x ) + 1 f<Tg e (tt) F 0 ( t ;X t;x ^t ;X t;x t ) i + 3 n X `=1 1 M n+1` Z T t e (rt) E h Y t;x;` r Y t;x;`1 r 2 i dr + 3 n X `=1 1 +M1 f`6=ng M n+1` Z T t e (rt) E h U ` (r;X t;x r )Y t;x;` r 2 i dr 70 For B n (t; x), we have: B n+1 (t; x) = E[U n+1 (t; x)]Y t;x;n+1 t 2 = E Z T t f(r;X t;x ^r ;U n (r;X t;x r ))f(r;X t;x ^r ;Y t;x;n r )dr 2 (Tt)E Z T t U n (r;X t;x r )Y t;x;n r 2 dr = (Tt)E Z T t U n (r;X t;x r )Y r;X t;x r ;n r 2 dr = (Tt) Z T t E h ~ E n (r;X t;x r ) i dr Combining the above two: ~ E n+1 (t; x) 1 M n+1 Var h g(X t;x ) + 1 f<Tg e (tt) F 0 ( t ;X t;x ^t ;X t;x t ) i + 3 n X `=1 1 M n+1` Z T t e (rt) E h Y t;x;` r Y t;x;`1 r 2 i dr + 3(1 +M) n1 X `=1 1 M n+1` Z T t e (rt) E h ~ E ` (r;X t;x r ) i dr + 3 Z T t (Tt) + e (rt) M E h ~ E n (r;X t;x r ) i dr 71 Under the Assumptions 2.1.1, Var h g(X t;x ) + 1 f<Tg e (tt) F 0 ( t ;X t;x ^t ;X t;x t ) i = Var h g(X t;x ) + 1 f<Tg e (tt) f( t ;X t;x ^t ; 0) i 2E h g 2 (X t;x ) + 1 f<Tg e 2(tt) 2 f 2 ( t ;X t;x ^t ; 0) i = 2E h g 2 (X t;x ) + Z T t e (rt) f 2 (r;X t;x ^r ; 0)dr i C T e (Tt) (1 + 1 2 )(1 +jxj 2 ) substitute the estimate in the Theorem 5.1.1 and above to the recursion ~ E n+1 (t; x) C T e (Tt) (1 +jxj 2 ) M n+1 (1 + 1 2 ) + C T e (Tt) (1 +jxj 2 ) n X `=1 (Tt) 2`1 M n+1` 2 `1 (` 1)! + 3(1 +M) n1 X `=1 1 M n+1` Z T t e (rt) E h ~ E ` (r;X t;x r ) i dr + 3 Z T t (Tt) + e (rt) M E h ~ E n (r;X t;x r ) i dr Since ~ E 1 (t; x) = 1 M Var h g(X t;x ) + 1 f<Tg e (tt) f( t ;X t;x ^t ; 0) i C ;T (1 +jxj 2 ) M E ~ E 1 (r;X t;x r ) C ;T (1 +EjX t;x r j 2 )C C ;T (1 +jxj 2 ) 72 And by induction, we can show that ~ En(t;x) C T (1+jxj ) ~ E n (t), where ~ E n (t) have the following recursion relation ~ E n+1 (t) e (Tt) M n+1 (1 + 1 2 ) + e (Tt) n X `=1 (Tt) 2`1 M n+1` 2 `1 (` 1)! + C M n1 X `=1 1 M n+1` Z T t e (rt) ~ E ` (r)dr +C (Tt) + 1 M Z T t e (rt) ~ E n (r)dr Let E n (t) := M n e t ~ En(Tt) (1+1= 2 ) , then E n (t) have the following recursion relation E n+1 (t) 1 + n X `=1 M ` t 2`1 2 ` (` 1)! + C M n1 X `=1 Z t 0 E ` (r)dr +C Mt + 1 Z t 0 E n (r)dr 1 + n X `=1 MT 2 ` t `1 (` 1)! +C max( M ;MT + 1 ) n X `=1 Z t 0 E ` (r)dr where E 1 (t) 1. To derive the bound for E n (t), we consider the recursion with equality ~ A n+1 (t) = 1 + n X `=1 MT 2 ` t `1 (` 1)! +c n X `=1 Z t 0 ~ A ` (r)dr = ~ A n (t) +c Z t 0 ~ A n (r)dr + MT 2 n t n1 (n 1)! 73 where ~ A 1 (t) = 1 and ~ A 2 (t) 1 + (M=2 +c)T =:C 0 . Then consider the recursion dened by the second equality starting from n = 2 with A 2 (t) = C 0 , then by induction, we can show that A n (t) has the following recursion relation: A n+1 (t) = n1 X k=0 a n+1;k t k ; a n+1;k =a n;k + c k a n;k1 ; a n+1;n1 = c n 1 a n;n2 + (MT ) n 2 n (n 1)! moreover we can show by induction that A n (t)C 0 2 n2 n2 X k=0 c k t k k! to check it plug the term in the recursion A n+1 (t)C 0 2 n2 n2 X k=0 c k t k k! +C 0 2 n2 n1 X k=1 c k t k k! + MT 2 n t n1 (n 1)! C 0 2 n2 1 + 2 n2 X k=1 c k t k k! + (c n1 + ( MT 2 ) n1 ) t n1 (n 1)! C 0 2 n2 1 + 2 n2 X k=1 c k t k k! + 2c n1 t n1 (n 1)! C 0 2 n1 n1 X k=0 c k t k k! where the last to the second inequality use c = C max( M ;MT + 1 ) MT=2. And it is straight forward to verify that for n 2 E n (t) ~ A n (t)A n (t)C 0 2 n2 n2 X k=0 c k t k k! 74 Therefore we have ~ E n (t) (1 + 1= 2 )e (Tt) M n E n (Tt) C 0 (1 + 1= 2 )e (Tt) 2 n2 M n n2 X k=0 c k (Tt) k k! C T e (Tt) 2 M n1 n2 X k=0 max( M ;MT + 1 ) k C (Tt) k k! and the estimate for ~ E n (t; x) ~ E n (t; x) C T (1 +jxj 2 )e (Tt) 2 M n1 n2 X k=0 max( M ;MT + 1 ) k C (Tt) k k! (5.5) combining with the estimate in Theorem 5.1.1, we derive the nal estimate E(U n (t; x)Y t;x t ) 2 ~ E n (t; x) +E(Y t;x;n t Y t;x t ) 2 C T (1 +jxj 2 )e (Tt) 2 M n1 n2 X k=0 max( M ;MT + 1 ) k C (Tt) k k! + C T (Tt) 2n+2 (1 +jxj 2 ) 2 n1 (n 1)! C T (1 +jxj 2 )e (Tt) 2 M n1 n1 X k=0 C k max( M ;MT + 1 ) k 1 fk6=n1g + (MT ) k 1 fk=n1g (Tt) k k! C T (1 +jxj 2 )e (Tt) 2 M n1 n1 X k=0 max( M ;MT + 1 ) k C (Tt) k k! Remark 5.1.3. The structure of the error bound (5.4) shows the advantage of the multilevel scheme is to attribute the fast convergence of Picard iteration more eectively, i.e. 1 M of smaller order k gets more reduction from factorial 1 (nk)! . 75 Compared to naive nested expectation with same decreasing simulations M n for each recursion, the multilevel scheme reduces the computation fromO(M n 2 ) to O(M n ), which is further illustrated in the next section, but the error rate does not degrade as much. Similar error analysis has been done in [16] for uniform stopping time. The analysis here is much conciser and provides slightly ner upper bound without having exponential term in M as well as the estimate of error to Picard iterations in (5.5), which could also serves as a sanity check of the result by letting M!1. 5.2 Computation Cost and Selection of Since the main computation cost comes from simulation of independent random variables, we will use the expected number of random variable to be simulated as a surrogate in our discussion. In addition, the analysis will focus on Markovian cases 76 with the forward process being Brownian motion to avoid unnecessary complica- tions. So the BSDE in (5.1) and its multilevel Picard Scheme in (5.3) are reduced to the following form: Y t =g(W T ) + Z T t f(r;W r ;Y r )dr Z T t Z r dW r U n+1 (t; x) = 1 M n+1 M n+1 X i=1 g(W 0;i T ) + 1 f 0;i t <Tg ( 0;i t t) f( 0;i t ;W 0;i 0;i t ;U 0 ( 0;i t ;W 0;i 0;i t )) + n X `=1 1 M n+1` M n+1` X i=1 1 f `;i t <Tg ( `;i t t) h f( `;i t ;W `;i `;i t ;U ` ( `;i t ;W `;i `;i t )) f( `;i t ;W `;i `;i t ;U `1 ( `;i t ;W `;i `;i t )) i : Let R n (t;!) := # independent r.v. simulated for U n (t; x;!); C n (t) :=E[R n (t;!)] note that the number of random variables simulated does not depend on the x, so R n (t;!) and C n (t) are both well dened. The rst proposition provides bounds on the expected number of simulated random variables. Proposition 5.2.1. The following estimate holds for C n (t) C n (t) (2M + 1) n 2 3 + 1 3(n 2)! Z (Tt) 3 0 x n2 e x dx 77 Proof. By the denition ofU n (t;x), the following recursion holds for C n (t) in the exponential density case: C n+1 (t) 2 n X `=0 M n+1` +(1+ 1 M ) n X `=1 M n+1` Z T t e (st) C ` (s)ds; C 1 (t) = 2M to simplify it, we can consider the recursion with equality, let C 1 (t) = 2e t and: C n+1 (t) = 2e t n X `=0 M ` +(1 + 1 M ) n X `=1 Z T t C ` (s)ds =C n (t) +(1 + 1 M ) Z T t C n (s)ds + 2e t M n then C n (t)e t C n (t)=M n , by induction, it can be shown that for n 2: C n (t) =e t a n e T b n;0 +b n;1 (Tt) + +b n;n2 (Tt) n2 78 plug the expression into the recursion, we have C n+1 (t) =e t a n+1 e T b n+1;0 +b n+1;1 (Tt) + +b n+1;n1 (Tt) n1 =e t (a n + 2M n )e T b n;0 +b n;1 (Tt) + +b n;n2 (Tt) n2 +(1 + 1 M ) h e t e T a n e T b n;0 (Tt) + b n;1 2 (Tt) 2 + + b n;n2 n 1 (Tt) n1 i =e t h (2 + 1 M )a n + 2M n i e T h b n;0 +a n (1 + 1 M ) + b n;1 +b n;0 (1 + 1 M ) (Tt) + + b n;n2 + b n;n3 n 2 (1 + 1 M ) (Tt) n2 + b n;n2 n 1 (1 + 1 M )(Tt) n1 i from which we can deduce the recursion relation between coecients a n+1 = (2 + 1 M )a n + 2M n ; a 1 = 2 b n+1;0 =b 2;0 + (1 + 1 M ) n X i=2 a i ; b n+1;` =b `+2;` + (1 + 1=M) ` n X i=`+2 b i;`1 since b 2;0 = (1 + 1 M )a 1 , b `+2;` = (1+1=M) ` b `+1;`1 , they can be simplied as a n+1 = (2 + 1 M )a n + 2M n ; a 1 = 2 b n+1;0 = (1 + 1 M ) n X i=1 a i ; b n+1;` = (1 + 1=M) ` n X i=`+1 b i;`1 79 and have the general expression as follows a n = (2 + 1 M ) n 1 M n (2 + 1 M ) n ; n 1 b n;k = k k! h (2 + 1 M ) n k X `=0 (1 + 1 M ) ` (2 + 1 M ) k+1` (nk 2) +` ` (1 + 1 M ) k+1 nk1 X i=1 1 M i n 1i k i k k! (2 + 1 M ) n by induction, it can be shown that b n;k k k! (2 + 1 M ) nk1 , thus C(t) is upper bounded as follows C(t)e t (2 + 1 M ) n e T (2 + 1 M ) n1 n2 X k=0 k (Tt) k (2 + 1=M) k k! e T (2 + 1 M ) n h e (Tt) 1 3 n2 X k=0 k (Tt) k 3 k k! i =e T (2 + 1 M ) n h e (Tt) e (Tt) 3(n 2)! Z 1 (Tt) 3 x n2 e x dx i =e t (2 + 1 M ) n h 2 3 + 1 3(n 2)! Z (Tt) 3 0 x n2 e x dx i therefore we establish the upper and lower bounds for C(t) C n (t) (2M + 1) n 2 3 + 1 3(n 2)! Z (Tt) 3 0 x n2 e x dx 80 Next proposition simplies the bound in Theorem 5.1.2 in settings when T is relatively large, which we are most interested in, and the bound shall be used in our further discussion about the optimality of . Proposition 5.2.2. Assume C MT 2 n, then E(U n (0; x)Y 0;x;n 0 ) 2 C T (n 1)(1 +jxj 2 )e T (C T ) n2 2 M n1 (n 2)! M 1 f< M1 TM g + 2MT1 f M1 TM g n2 where C is the same one in Theorem 5.1.2. Proof. Based on the estimate in (5.4), it suces to show that n2 X k=0 max( M ;MT + 1 ) k C T k k! (n 1)(C T ) n2 (n 2)! M 1 f< M1 TM g +2MT1 f M1 TM g n2 : if < M1 TM , then for any k, max( M ;MT + 1 ) = M ; C MT k C M 2 T 2 (M 1)k C MT 2 n 1; if M1 TM , again for any k, max( M ;MT + 1 ) =MT + 1 ; MT + 1 TM 2 M 1 2MT; 2C MT 2 k C MT 2 n 1: 81 therefore the term inside the summation is increasing in both cases, the estimate follows by bounding all of them by the last term. To show the optimal is small when T is large, we xed the n and minimize the upper bound in Proposition 5.2.2 over the andM, dropping other unrelated terms. This is what the following Proposition says. Proposition 5.2.3. For any constant C 1, let ( ;M ) := arg min ;M:M n 2 3 + 1 3(n2)! R T 0 x n2 e x dx =C e T 2 M ; then < 2 T . Proof. First notice that can be viewed as the minimizer of the unconstrained problem = arg min '(); where '() = e T 2 2 3 + 1 3(n 2)! Z T 0 x n2 e x dx 1 n for any xed T > 0, it is straightforward to check that lim #0 '() = lim !+1 '() = +1 82 so the minimum is attainable and is well-dened. Taking the derivative to nd ' 0 () =e T T 2 2 3 2 3 + 1 3(n 2)! Z T 0 x n2 e x dx 1 n + n4 T n1 n 2 3 + Z T 0 x n2 e x dx n1 n = 2 3 + Z T 0 x n2 e x dx n1 n h e T T 2 2 3 2 3 + Z T 0 x n2 e x dx + n4 T n1 n i for ' 0 () = 0, so it is necessary to have T 2 2 3 < 0 =) < 2 T : Remark 5.2.4. Proposition 5.2.2 and 5.2.3 basically suggests for given computa- tion resources, as T gets larger, should be made smaller and then M is picked accordingly, to minimize the error bound, which we sort of take as an approxi- mation, for E(U n (0; x)Y 0;x;n 0 ) 2 . This is in spirit the same as the philosophy in structured risk minimization in VC theory. 83 Chapter 6 Numerical Examples This chapter is dedicated into numerical examples comparing three dierent meth- ods we studied and each one is designed to demonstrate some point we make in previous chapters about certain characteristics of the method. Both implementa- tion setup and details are given by example. All the results are run on a Thinkpad with Intel Core i7-5600U CPU. Example 6.0.1. Regular Kernel with Polynomial Generator Consider FBSDE (4.1) with: b i (s;t; x) = ii (s;t; x) =e 1 2 (st) 2 x i r ; ij (s;t; x) = 0 for i6=j g(x) = cos d X i=1 x i T + p d 10 ; f(t; x;y;z) = 0:5(y +y 3 ): results for x = 1;T = 0:5 from dierent methods are presented as follows branching: (t) = e t , p 1 = p 3 = 0:5 and each row is based on 100000 simulations of branching processes (representation) 84 Dimension Approx Y 0 Standard deviation Runtime (sec) 10 0.443666 0.003924 5.8 20 0.621878 0.004323 10.9 40 0.94833 0.005657 21.4 80 1.66949 0.032471 39.6 multilevel picard: M =n = 6, = 1 and each row is based on 10 simulations of U n (0; x) Dimension Approx Y 0 Standard deviation Runtime (sec) 10 0.437708 0.011431 61.3 20 0.620444 0.017917 131.3 40 0.947066 0.049139 251.8 80 1.74911 0.365404 483.8 deep learning based: t = 0:02, Z is approximated by a two layer NN with d+10 hidden units each, Adam optimizer with learning rate 510 4 is used, and each row is based on 10000 iterations with batch size = 64 (except last row is t = 0:01, 20000 iterations) 85 Dimension Approx Y 0 Standard deviation Loss Runtime (sec) 10 0.438486 0.002698 0.5839 88.9 20 0.624945 0.002298 0.4755 112.16 40 0.949345 0.003242 0.6067 166.5 80 1.70269 0.003075 0.4957 324.7 80 1.66776 0.004469 0.4385 2374.8 the rst example is intended to show that all three methods work wheng andf are relatively small. For this case, the branching method is the most ecient one in terms of runtime and accuracy, multilevel Picard and deep learning based method perform similarly, the former exhibits larger variance as dimension increases, the later shows noticeable bias but this could be solved by reducing time discretization size and increasing iterations. It is also worth noting that deep learning method approximates the solution well, even though the loss does not decrease to 0. Example 6.0.2. Regular Kernel with Non-Polynomial Generator Consider FBSDE (4.1) with: b i (s;t; x) = ii (s;t; x) =e 1 2 (st) 2 x i r ; ij (s;t; x) = 0 for i6=j g(x) = h min 1id x i T 30 i + h min 1id x i T 60 + 15; f(t; x;y;z) = minfy; 0g: results for = 0:1; x = 100;T = 1 from dierent methods are presented as follows 86 multilevel picard: M =n = 6, = 1 and each row is based on 10 simulations of U n (0; x) Dimension Approx Y 0 Standard deviation Runtime (sec) 10 -7.63928 0.086653 450.0 20 -12.1665 0.027597 841.9 40 -13.4780 0.014284 2089.3 80 -13.5666 0.00807 3500.3 deep learning based: t = 0:04, Z is approximated by a two layer NN with d+10 hidden units each, Adam optimizer with learning rate 510 4 is used, and each row is based on 50000 iterations with batch size = 64 Dimension Approx Y 0 Standard deviation Loss Runtime (sec) 10 -7.69397 0.005803 68.91 449.2 20 -12.1799 0.009940 21.73 497.1 40 -13.4757 0.004680 0.2480 766.3 80 -13.5745 0.001337 0.00003 1616.3 whenf is not polynomial, branching is no longer applicable, both multilevel Picard and deep learning method approximate the solution well as dimension increases. But by comparing with the rst example, we can note that the runtime scales roughly cubic for multilevel Picard and quadratic for deep learning in T , this is expected since the forward path simulation for Volterra SDE is already quadratic in 87 T , and increasing number of particles and model parameters contribute to the extra costs. To further illustrate the sensitivity to T of both methods, we can increase the instead, which increases the Lipschitz constant and is roughly equivalent to increasingT , to avoid long waiting time, the result of same setup ford = 10; = 0:5 is shown as follows Method Approx Y 0 Standard deviation Loss Runtime (sec) Multilevel Picard -4.75103 0.16502 NA 444.1 Deep Learning -5.16125 0.012077 94.13 406.7 the standard deviation for multilevel Picard is no longer negligible, the deep learn- ing still converges but it is hard to verify the approximation is indeed the true solution. In this example, backward non-linearity f and terminal g mimics the pricing example with counterparty credit risk in [17], but the forward process is replaced by Volterra SDE and thus the backward SDE becomes path dependent automatically. The modeling meaning is not quite clear, since the X is not guar- anteed to be positive, but we are merely interested in the performance of two algorithms. Example 6.0.3. Singular Kernel with Polynomial Generator Consider FBSDE (4.1) with: b i (s;t; x) = 0:5; ii (s;t; x) = (st) 0:2 ; ij (s;t; x) = 0 for i6=j 88 g(x) = p d 1 +jx T j 2 2 ; f(t; x;y;z) = 0:5(y +y 2 +y 3 ): results for x = 0;T = 0:5 from dierent methods are presented as follows branching: (t) = e t , p 1 = p 2 = p 3 = 1 3 and each row is based on 100000 simulations of branching processes (representation) Dimension Approx Y 0 Standard deviation Runtime (sec) 10 0.495727 0.001087 8.4 20 0.316559 0.000640 16.8 40 0.211855 0.000420 33.8 80 0.14448 0.000289 62.3 multilevel picard: M =n = 6, = 1 and each row is based on 10 simulations of U n (0; x) Dimension Approx Y 0 Standard deviation Runtime (sec) 10 0.491668 0.010663 89.1 20 0.31702 0.003839 171.4 40 0.211086 0.001701 353.5 80 0.144586 0.001007 656.95 deep learning based: t = 0:02, Z is approximated by a two layer NN with d+10 hidden units each, Adam optimizer with learning rate 510 4 is used, and each row is based on 6000 iterations with batch size = 64 89 Dimension Approx Y 0 Standard deviation Loss Runtime (sec) 10 0.480746 0.002755 0.0111 50.0 20 0.308088 0.002106 0.0029 58.7 40 0.206572 0.001261 0.0006 90.8 80 0.143113 0.000818 0.0002 132.7 since the main motivation for Volterra SDE is fractional Brownian motion (fBM), whose kernel is of order (st) H as t! s, it is natural to test the performance of these methods for such forward process. The results are very similar to the Example (6.0.1): branching is the best and the other two also approximate solu- tions very well. In fact according to our experience, the regularity of the kernel is typically less a concern for the convergence of these methods, than the size of f andg or their regularity, which can cause integrability issue or large variance. The solution for this example seems to atten as dimension increases due to the form of the terminal condition, so all the approximations tend to have less variance as dimension increases and the bias of deep learning decreases as well. Example 6.0.4. Singular Kernel in Recursive Pricing with Credit Risk Consider FBSDE (4.1) with: b i (s;t; x) = 0:5; ii (s;t; x) = (st) 0:2 ; ij (s;t; x) = 0 for i6=j g(x) =d min 1id e x i T ; f(t; x;y;z) =(1) minf h ; maxf l ; ( h l ) v h v l yv h + h ggRy: 90 results for = 2 3 ; h = 0:2; l = 0:02;v h = 47;v l = 65; x = 0;T = 0:5 from dierent methods are presented as follows multilevel picard: M =n = 6, = 1 and each row is based on 10 simulations of U n (0; x) Dimension Approx Y 0 Standard deviation Runtime (sec) 10 3.14419 0.010650 110.1 20 4.39056 0.012111 177.1 40 6.43652 0.013514 377.4 80 9.73755 0.023133 689.31 multilevel picard: M =n = 6, = 0:5 and each row is based on 10 simula- tions of U n (0; x) Dimension Approx Y 0 Standard deviation Runtime (sec) 10 3.13879 0.007104 59.2 20 4.39155 0.013761 112.9 40 6.42724 0.018202 212.0 80 9.73685 0.016279 475.5 deep learning based: t = 0:02, Z is approximated by a two layer NN with d+10 hidden units each, Adam optimizer with learning rate 510 4 is used, and each row is based on 10000 iterations with batch size = 64 91 Dimension Approx Y 0 Standard deviation Loss Runtime (sec) 10 3.08840 0.014758 0.74 83.2 20 4.31817 0.020544 4.61 100.1 40 6.27030 0.032328 6.62 163.9 80 9.43981 0.042955 12.51 320.8 this example aims to test these methods on the recursive pricing model with default risk due to [17], and here the conventional geometric Brownian motion is replaced by the exponential of Volterra SDE with similar singularity as in the kernel of fBM. As the results suggest, multilevel picard is less ecient but more reliable, in particular because we can prove its convergence rigorously; while the deep learning based method exhibits larger bias as dimension increases, which cannot be simply recovered in this case by reducing the discretization size. In addition, results for multilevel Picard under dierent choices of are also compared, and typically when coecients are nice, mildly small achieves best balance between run-time and accuracy. Example 6.0.5. Transformation for Blow-up in Quadratic Case Consider FBSDE (4.1) with: b i (s;t; x) = 0; ii (s;t; x) = (st) 0:2 x i r ; ij (s;t; x) = 0 for i6=j g(x) = jx T j 2 1 +jx T j 2 ; f(t; x;y;z) =y 0:5y 2 : 92 and one of its transformation Y s =e 2s (1Y s ), where corresponding ( f; g) are g(x) =e 2T 1 1 +jx T j 2 ; f(t; x;y;z) = 1:5e 2t + 0:5e 2t y 2 : results for x = 1;T = 1 from dierent methods are presented as follows branching via original Y and transformed Y : (t) =e t , p 1 =p 2 = 0:5 and each row is based on 100000 simulations of branching processes (representa- tion) Dimension Approx Y 0 Standard deviation Runtime (sec) 10 0.4543 0.1879 31.8 20 0.0449 0.1758 70.6 Dimension Approx Y 0 Approx Y 0 = 1 Y 0 Std Runtime (sec) 10 0.770933 0.229067 0.003848 19.7 20 0.757682 0.242318 0.003562 35.8 multilevel picard via original Y and transformed Y : M =n = 6, = 1 and each row is based on 10 simulations of U n (0; x) Dimension Approx Y 0 Standard deviation Runtime (sec) 10 0.13186 0.090706 620.4 20 0.266536 0.143286 1229.6 93 Dimension Approx Y 0 Approx Y 0 = 1 Y 0 Std Runtime (sec) 10 0.769447 0.230553 0.009995 624.0 20 0.749483 0.250517 0.009139 1261.74 deep learning based via original Y and transformed Y : t = 0:02, Z is approximated by a two layer NN with d + 10 hidden units each, Adam opti- mizer with learning rate 5 10 4 is used, and each row is based on 10000 iterations with batch size = 64 Dimension Approx Y 0 Standard deviation Loss Runtime (sec) 10 0.233362 0.000618 0.00567 207.1 20 0.251922 0.000476 0.00290 267.7 Dimension Approx Y 0 Approx Y 0 = 1 Y 0 Std Loss Runtime (sec) 10 NaN NA NA NA NA 20 NaN NA NA NA NA in Remark (4.3.4), we propose to nd alternative transformation with positive g and positive coecients in f to avoid the blow-up issue for branching method. This example is to illustrate this phenomenon, and we choose quadratic since it is always possible to nd such transformation in this case as suggested in Proposition 4.3.3. Branching does not seem to converge for the original BSDE even if the number of simulation is increased to 400000, but it converges for the transformed one well with best eciency; surprisingly multilevel picard seems to diverge or have too 94 large variance for the original BSDE as well, but converges for the transformed one; the deep learning based method serves as a complementary conrmation to the results, but it works in the opposite way converges for the original BSDE but the forward simulation blow up for the transformed one, which is why Y is NaN. The reason is the forward simulation in the transformed BSDE has large coecient in y 2 and the forward Euler scheme is equivalent to a Riccati equation with large positive leading coecients. Remark 6.0.6. To summarize, we present the following table that ranks the overall performance of the three methods, in terms of generality, eciency and reliability of results, there are some cases where the ranking is not precise, but they hold in general based on our experience in testing these methods. Overall they are the most promising numerical methods for Volterra type FBSDEs or path-dependent BSDEs more generally, but further researches are needed to reduce the complexity in time horizon, in order to make these methods truly widely applicable. Applicability Eciency Dimensionality Reliability of Results Multilevel Picard Branching (if applicable) Branching (if applicable) Branching (if applicable) Deep Learning Deep Learning Deep Learning Multilevel Picard Branching Multilevel Picard Multilevel Picard Deep Learning 95 Bibliography [1] Athreya, K. and Ney, P. (1972). Branching processes. Die Grundlehren der mathematischen Wissenschaften. 196 Springer-Verlag, New York. [2] Bender, C. and Denk, R. (2007). A Forward Scheme for Backward SDEs. Stoch. Process. Appl. 117 1793-1812. [3] Beck, C., E, W., and Jentzen, A. (2017). Machine learning approximation algorithms for high dimensional fully nonlinear partial dierential equations and second-order backward stochastic dierential equations. Journal of Non- linear Science. 1-57. [4] Berger, M. and Mizel, V. (1980). Volterra equation with Ito integrals, I, II. J Intergal Equations. 2 187-245; 319-337. [5] Bally, V., Pag es, G. and Printems, J. (2005). A Quantization Tree Method for Pricing and Hedging Multidimensional American Options. Math. Finance. 15 119-168. [6] Bouchard, B. and Touzi, N. (2004). Discrete-time Approximation and Monte- Carlo Simulation of Backward Stochastic Dierential Equations. Stochastic Process. Appl. 111 175-206. [7] Bouchard, B., Tan, X. and Warin, X. (2017). Numerical approximation of general Lipschitz BSDEs with branching processes. ESAIM: Proceedings and Surveys. 65 309-329. [8] Bouchard, B., Tan, X., Warin, X. and Zou, Y. (2017). Numerical approxima- tion of BSDEs using local polynomial drivers and branching processes. Monte Carlo Methods and Applications. 23 (4). 96 [9] Coutin, L. and Decreusefond, L. (2001). Stochastic Volterra Equations with Singular Kernels. Stochastic Analysis and Mathematical Physics. 50 39-50. [10] Cont, R. and Fourni e, D. (2013). Functional It^ o calculus and stochastic inte- gral representation of martingales. Ann. Probab. 41 109-133. [11] Crisan, D. and Manolarakis, K. (2012). Solving Backward Stochastic Dif- ferential Equations Using The Cubature Method: Application To Nonlinear Pricing. SIAM J. Financial Math. 1 534-571. [12] Cheridito, P., Soner, H., Touzi, N. and Victoir, N. (2007). Second-order back- ward stochastic dierential equations and fully nonlinear parabolic PDEs. Comm. Pure Appl. Math. 60 1081-1110. [13] Dupire, B. (2009). Functional It^ o calculus. Preprint. Available at papers.ssrn.com. [14] E, W. and Han, J. (2016). Deep learning Approximation for Stochastic Control Problems. Available at arXiv:1611.07422. [15] E, W., Han, J., and Jentzen, A. (2017). Deep learning-based numerical meth- ods for highdimensional parabolic partial dierential equations and backward stochastic dierential equations. Communications in Mathematics and Statis- tics 5 349-380. [16] E, W., Hutzenthaler, M., Jentzen, A., and Kruse, T. (2016). Multilevel Picard iterations for solving smooth semilinear parabolic heat equations. Available at arXiv:1607.03295. [17] E, W., Hutzenthaler, M., Jentzen, A., and Kruse, T. (2019). On mul- tilevel Picard numerical approximations for high-dimensional nonlinear parabolic partial dierential equations and high-dimensional nonlinear back- ward stochastic dierential equations. J. Sci. Comput. 79 (3) 1534-1571 [18] Ekren, I.; Keller, C., Touzi, N. and Zhang, J. (2014). On viscosity solutions of path dependent PDEs. Ann. Probab. 42 (1), 204-236. [19] EL Karoui, N., Peng, S. and Quenez, M. C. (1997). Backward Stochastic Dierential Equations in Finance. Math. Finance. 7, 1-71. [20] El, E. and Rosenbaum, M. (2016). The characteristic function of rough Heston models. Available at arXiv:1609.02108. [21] Ekren, I., Keller, C., Touzi, N. and Zhang, J. (2014). Optimal Stopping under Nonlinear Expectation. Stochastic Processes and Their Applications, 124, 3277- 3311. 97 [22] Ekren, I., Keller, C., Touzi, N. and Zhang, J. (2016). Viscosity solutions of fully nonlinear parabolic path dependent PDEs: part I. Ann. Probab. 44 (2), 1212-1253. [23] Ekren, I., Keller, C., Touzi, N. and Zhang, J. (2016). Viscosity solutions of fully nonlinear parabolic path dependent PDEs: part II. Ann. Probab. 44 (4), 2507-2553. [24] Fahim, A., Touzi, N. and Warin, X. (2011). A Probabilistic Numerical Method for Fully Nonlinear Parabolic PDEs. Ann. Appl. Probab. 21 1322-1364. [25] Goodfellow, I., Bengio, Y. and Courville, A. (2016). Deep Learning. MIT Press. [26] Gatheral, J., Jaisson, T. and Rosenbaum, M. (2014). Volatility is rough. Avail- able at arXiv:1410.3394. [27] Gobet, E., Lemor, J. and Waxin, X. (2005). A Regression-based Monte Carlo Method to Solve Backward Stochastic Dierential Equations. Ann. Appl. Probab. 15 2172-2202. [28] Guo, W., Zhang, J. and Zhuo, J. (2015). A monotone scheme for high dimen- sionaly fully nonlinear PDEs. Ann. Appl. Probab. 25 1540-1580. [29] Harris, T. (1963). The Theory of Branching Processes. Die Grundlehren der Mathematischen Wissenschaften. 119 Springer, Berlin. [30] Henry-Labord ere, P. (2012). Cutting CVA's Complexity. Risk magazine. [31] Henry-Labord ere, P., Oudjane, N., Tan, X., Touzi, N., and Warin, X. (2016). Branching diusion representation of semilinear PDEs and Monte Carlo approximation. Available at arXiv:1603.01727. [32] Henry-Labord ere, P., Tan, X. and Touzi, N. (2014). A numerical algorithm for a class of BSDEs via the branching process. Stochastic Processes and their Applications. 124 (2) 1112-1140. [33] Ma, J., Protter, P. and Yong, J. (1994). Solving forward-backward stochastic dierential equations explicitly-a four step scheme. Probab. Theory Related Fields 98 339-359. [34] Ma, J., Zhang, J. and Zheng, Z. (2008). Weak solutions for forward-backward SDEs-a martingale problem approach. Ann. Probab. 36 2092-2125. [35] McKean, H. (1975). Application of Brownian motion to the equation of Kolmogorov-Petrovskii-Piskunov. Comm. Pure Appl. Math. 28 323-331. 98 [36] Peng, S. (1990). Probabilistic Interpretation for systems of quasilinear parabolic partial dierential equations. Stochastics and Stochastics Reports. 37 61-74. [37] Pardoux, E. and Peng, S. (1990). Adapted solution of a backward stochastic dierential equation. Systems Control Lett. 14 55-61. [38] Pardoux, E. and Peng, S. (1992). Backward stochastic dierential equations and quasilinear parabolic partial dierential equations. In Stochastic Partial Dierential Equations and Their Applications (Charlotte, NC, 1991). Lecture Notes in Control and Information Sciences 176 200-217. Springer, Berlin. [39] Pardoux, E. and Protter, P. (1990). Stochastic Volterra Equations with Antic- ipating Coecients. Ann. Probab. 18 1635-1655. [40] Pardoux, E. and Tang, S. (1999). Forward-backward stochastic dierential equations and quasilinear parabolic PDEs. Probab. Theory Related Fields 114 123-150. [41] Peng, S. (2010). Backward stochastic dierential equation, nonlinear expec- tation and their applications. In Proceedings of the International Congress of Mathematicians. Volume I 393-432. Hindustan Book Agency, New Delhi. [42] Protter, P. (1985). Volterra Equations Driven by Semimartingales. Ann. Probab. 13 519-530. [43] Ren, Z., Tan, X. (2017). On the convergence of monotone schemes for path- dependent PDEs. Stochastic Process. Appl. 127 (6) 1738-1762. [44] Soner, H., Touzi, N. and Zhang, J. (2012). Wellposedness of second order backward SDEs. Probab. Theory Related Fields 153 149-190. [45] Viens, F. and Zhang, J. (2019). A martingale approach for fractional Brownian motions and related path dependent PDEs. Ann. Applied Probab. 29, 3489- 3540. [46] Wang, Z. (2007). Existence and uniqueness of solutions to stochastic Volterra equations with singular kernels and non-Lipschitz coecients. Statistics & Probability Letters. 78 (9) 1062-1071. [47] Zhang, X. (2008). Euler schemes and large deviations for stochastic Volterra equations with singular kernels. Journal of Dierential Equations. 244 (9) 2226-2250. [48] Zhang, J. (2004). A Numerical Scheme for Backward Stochastic Dierential Equations. Ann. Appl. Probab. 14 459-488. 99

Abstract (if available)

Abstract In this dissertation we study the stochastic Volterra equations and numerical methods for the associated path-dependent PDEs (PPDEs for short). It consists of two main parts. We first prove the wellposedness of path-dependent stochastic Volterra equations, and generalize the estimate results and flow properties of the classical SDEs to the Volterra ones. In the second part, we extend three different numerical methods for high-dimensional PDEs to our path-dependent scenarios: branching diffusion method, multilevel Picard approximation and deep learning based method. For the branching diffusion method, we explain the blow-up issues and introduce a transformation to extend the time duration, which is especially effective in quadratic cases. For the multilevel Picard approximation, we introduce an extra parameter in stopping strategy, which will help extend the time duration. In the end, we present several high-dimensional numerical examples with different emphasis, and give a comparison of the three methods as well as the guidelines about which method to choose in different scenarios, in terms of their applicability, scalability and reliability.

Linked assets

University of Southern California Dissertations and Theses

Conceptually similar

PDF

Probabilistic numerical methods for fully nonlinear PDEs and related topics

PDF

Pathwise stochastic analysis and related topics

PDF

Tamed and truncated numerical methods for stochastic differential equations

PDF

Topics on set-valued backward stochastic differential equations

PDF

Monte Carlo methods of forward backward stochastic differential equations in high dimensions

PDF

Conditional mean-fields stochastic differential equation and their application

PDF

Forward-backward stochastic differential equations with discontinuous coefficient and regime switching term structure model

PDF

Zero-sum stochastic differential games in weak formulation and related norms for semi-martingales

PDF

On non-zero-sum stochastic game problems with stopping times

PDF

On spectral approximations of stochastic partial differential equations driven by Poisson noise

$Stochastic differential equations driven by fractional Brownian motion and Poisson jumps$

PDF

Stochastic differential equations driven by fractional Brownian motion and Poisson jumps

PDF

Some topics on continuous time principal-agent problem

PDF

Path dependent partial differential equations and related topics

PDF

Statistical inference for second order ordinary differential equation driven by additive Gaussian white noise

PDF

Set values for mean field games and set valued PDEs

PDF

Reproducible large-scale inference in high-dimensional nonlinear models

PDF

Dynamic approaches for some time inconsistent problems

PDF

Large-scale inference in multiple Gaussian graphical models

PDF

Reinforcement learning for the optimal dividend problem

PDF

Generalized linear discriminant analysis for high-dimensional genomic data with external information

Asset Metadata

Creator Ruan, Jie (author)

Core Title Numerical methods for high-dimensional path-dependent PDEs driven by stochastic Volterra integral equations

School College of Letters, Arts and Sciences

Degree Doctor of Philosophy

Degree Program Applied Mathematics

Publication Date 07/08/2020

Defense Date 05/28/2020

Publisher University of Southern California (original), University of Southern California. Libraries (digital)

Tag high-dimensional path-dependent PDEs,numerical methods,OAI-PMH Harvest,stochastic Volterra equations

Language English

Contributor Electronically uploaded by the author (provenance)

Advisor Zhang, Jianfeng (committee chair), Lv, Jinchi (committee member), Ma, Jin (committee member)

Creator Email jieruan@usc.edu

Permanent Link (DOI) https://doi.org/10.25549/usctheses-c89-324809

Unique identifier UC11665365

Identifier etd-RuanJie-8638.pdf (filename),usctheses-c89-324809 (legacy record id)

Legacy Identifier etd-RuanJie-8638.pdf

Dmrecord 324809

Document Type Dissertation

Rights Ruan, Jie

Type texts

Source University of Southern California (contributing entity), University of Southern California Dissertations and Theses (collection)

Access Conditions The author retains rights to his/her dissertation, thesis or other graduate work according to U.S. copyright law. Electronic access is being provided by the USC Libraries in agreement with the a...

Repository Name University of Southern California Digital Library

Repository Location USC Digital Library, University of Southern California, University Park Campus MC 2810, 3434 South Grand Avenue, 2nd Floor, Los Angeles, California 90089-2810, USA