Analytical Approximation for Non-linear FBSDEs with Perturbation Scheme ∗ 2 1 Masaaki Fujii, Akihiko Takahashi † ‡ 0 2 First version: June 1, 2011 n This version: January 20, 2012 a J 0 2 ] P C Abstract . n i In this work, we have presented a simple analytical approximation scheme for f - genericnon-linearFBSDEs. Bytreatingthe interestedsystemasthelineardecoupled q FBSDE perturbed withnon-lineargeneratorandfeedbackterms, wehaveshownthat [ it is possible to carry out a recursive approximation to an arbitrarily higher order, 3 where the required calculations in each order are equivalent to those for standard v European contingent claims. We have also applied the perturbative method to the 3 PDE framework following the so-called Four Step Scheme. The method is found to 2 render the original non-linear PDE into a series of standard parabolic linear PDEs. 1 0 Duetotheequivalenceofthetwoapproaches,itisalsopossibletoderiveapproximate . analytic solution for the non-linear PDE by applying the asymptotic expansion to 6 0 the corresponding probabilistic model. Two simple examples are provided to demon- 1 strate how the perturbationworksandshowits accuracyrelativeto knownnumerical 1 techniques. The method presented in this paper may be useful for various important v: problems which have eluded analytical treatment so far. i X r a Keywords: BSDE,FBSDE,FourStepScheme,AsymptoticExpansion,MalliavinDeriva- tive, Non-linear PDE, CVA ∗This research is supported by CARF (Center for Advanced Research in Finance) and the global COE program “The research and training center for new development in mathematics.” All the contents expressed in this research are solely those of the authors and do not represent any views or opinions of any institutions. The authors are not responsible or liable in any manner for any losses and/or damages caused by theuse of any contents in this research. †Graduate School of Economics, TheUniversity of Tokyo ‡Graduate School of Economics, TheUniversity of Tokyo 1 1 Introduction In this paper, we propose a simple analytical approximation for backward stochastic dif- ferential equations (BSDEs). These equations were introduced by Bismut (1973) [1] for the linear case and later by Pardoux and Peng (1990) [10] for the general case, and have earned strong academic interests since then. They are particularly relevant for the pricing of contingent claims in constrained or incomplete markets, and for the study of recursive utilities as presented by Duffie and Epstein (1992) [2]. For a recent comprehensive study with financial applications, one may consult Yong and Zhou (1999) [15], Ma and Yong (2000) [9] and references therein. The importance of BSDEs, or more specifically non-linear FBSDEs which have non- linear generators coupled with some state processes satisfying the forward SDEs, has risen greatly in recent years also among practitioners. The collapse of major financial institu- tionsfollowed bythedrasticreformofregulations makethemwellawareoftheimportance ofcounterparty riskmanagement, creditvalueadjustments(CVA)inparticular. Evenina very simple setup, if there exists asymmetry in the credit risk between the two parties, the relevantdynamicsofportfoliovaluefollowsanon-linearFBSDEasclearly shownbyDuffie and Huang(1996) [3]. We have recently foundthat the asymmetrictreatment of collateral between the two parties also leads to a non-linear FBSDE [5]. Furthermore, in May 2010, regulators were forced to realize the importance of mutual interactions and feedback loops in the trading activities among financial firms, shocked by the astonishing flash crash of the Dow Jones index by almost 1,000 points. Once we take the feedback effects from the behavior of major players into account, we naturally end up with complicated coupled FBSDEs. Unfortunately, however, an explicit solution for a FBSDE is only known for a simple linearexample. Inthelastdecade, severaltechniqueshavebeenintroducedbyresearchers, but they tend to be quite complicated for practical applications. They either require one to solve non-linear PDEs, which are very difficult in general, or resort to quite time- consumingsimulation. Although regression based Monte Carlosimulation hasbeenrather popular among practitioners for the pricing of callable products, the appropriate choice of regressors and attaining numerical stability becomes a more subtle issue for a general FBSDE. In fact, in clear contrast to the pricing of callable products, one cannot tell if the price goes up or down when one improves the regressors, which makes it particularly difficult to select the appropriate basis functions. In this paper, we present a simple analytical approximation scheme for the non-linear FBSDEs coupled with generic Markovian state processes. We have perturbatively ex- panded the non-linear terms around the linearized FBSDE, where the expansion can be made recursively to an arbitrary higher order. In each order of approximation, the re- quired calculations are equivalent to those for standard European contingent claims. In order to carry out the perturbation scheme, we need to express the backward components explicitly in terms of the forward components for each order of approximation. For that purpose, we propose to use the asymptotic expansion of volatility for the forward com- ponents, which is now widely adopted to price various European contingent claims and compute optimal portfolios (See, for examples [7, 11, 12, 13, 14] and references therein for the recent developments and review.). In the case when the underlying processes have known distributions, of course, we can directly proceed to a higher order approximation 2 without resorting to an asymptotic expansion. We have also studied a perturbation scheme in the PDE framework, or in the so- called Four Step Scheme [8], for the generic fully-coupled non-linear FBSDEs. We have shown that our perturbation method renders the original non-linear PDE into the series of classical linear parabolic PDEs, which are easy to handlewith standard techniques. We thenprovidedthecorrespondingprobabilisticframeworkbyusingtheequivalencebetween the two approaches. We have shown that, also in this case, the required calculations in a given order are equivalent to those for the classical European contingent claims. As a by- product, by applying the asymptotic expansion method to the correspondingprobabilistic model, it was actually found possible to derive an analytic expression for the solution of the non-linear PDE up to a given order of perturbation. Therefore, our method can be interpreted as a practical implementation of the Four Step Scheme in the perturbative approach. The organization of of the paper is as follows: In Section 2, we will explain our new approximation scheme with perturbative expansion for generic decoupled FBSDEs. Then, in Section 3, we shall apply it to the two concrete examples to demonstrate how it works and test its numerical performance. One of them allows a direct numerical treatment by a simple PDE and hence it is easy to compare the two methods. In the second example, we will consider a slightly more complicated model. We compare our approximation result to the detailed numerical study recently carried out by Gobet et al. (2005) [6] using a regression-based Monte Carlo simulation. In Section 4, we explain how to use standard asymptotic expansion procedures to express the backward components explicitly when the forward components do not have known distributions. In Section 5, we will give an extension of our method to the fully coupled non-linear FBSDEs under the PDE frame- work, and then formulate the equivalent probabilistic approach in Section 5. Appendix contains slightly different scheme for coupled non-linear FBSDEs which may be useful for the actual application. 2 Approximation Scheme 2.1 Setup Let us briefly describe the basic setup. The probability space is taken as (Ω, ,P) and F T (0, ) denotes somefixedtime horizon. W = (W1, ,Wr) , 0 t T isRr-valued t t t ∗ ∈ ∞ ··· ≤ ≤ Brownian motion defined on (Ω, ,P), and ( ) stands for P-augmented natural t 0 t T F F { ≤ ≤ } filtration generated by the Brownian motion. We consider the following forward-backward stochastic differential equation (FBSDE) dV = f(X ,V ,Z )dt+Z dW (2.1) t t t t t t − · V = Φ(X ) (2.2) T T where V takes the value in R, and X Rd is assumed to follow a generic Markovian t ∈ forward SDE dX = γ (X )dt+γ(X ) dW . (2.3) t 0 t t t · Here,weabsorbedanexplicitdependenceontimetoX byallowingsomeofitscomponents can be a time itself. Φ(X ) denotes the terminal payoff where Φ(x) is a deterministic T 3 function of x. The following approximation procedures can be applied in the same way also in the presenceof coupon payments. Z and γ take values in Rr andRd r respectively, × and ”” in front of the dW represents the summation for the components of r-dimensional t · Brownian motion. Throughout this paper, we are going to assume that the appropriate regularity conditions are satisfied for the necessary treatments. 2.2 Perturbative Expansion for Non-linear Generator In order to solve the pair of (V ,Z ) in terms of X , we extract the linear term from t t t the generator f and treat the residual non-linear term as the perturbation to the linear FBSDE. We introduce the perturbation parameter ǫ, and then write the equation as (ǫ) (ǫ) (ǫ) (ǫ) (ǫ) dV = c(X )V dt ǫg(X ,V ,Z )dt+Z dW (2.4) t t t − t t t t · t (ǫ) V = Φ(X ) , (2.5) T T where ǫ = 1 corresponds to the original model by 1 f(X ,V ,Z ) = c(X )V +g(X ,V ,Z ) . (2.6) t t t t t t t t − Usually, c(X ) corresponds to the risk-free interest rate at time t, but it is not a necessary t condition. One should choose the linear term in such a way that the residual non-linear term becomes as small as possible to achieve better convergence. A possible linear term θ(X)Z in the driver f can be absorbed by the measure change and hence the simple reinterpretation of the drift term of the forward components γ results in the form (2.4). 0 See also the discussion in Appendix. Now, we are going to expand the solution of BSDE (2.4) and (2.5) in terms of ǫ: that (ǫ) (ǫ) is, suppose V and Z are expanded as t t V(ǫ) = V(0) +ǫV(1)+ǫ2V(2)+ (2.7) t t t t ··· Z(ǫ) = Z(0)+ǫZ(1)+ǫ2Z(2)+ . (2.8) t t t t ··· Onceweobtain thesolution uptothecertain order,say k forexample, thenbyputting ǫ = 1, k k V˜ = V(i), Z˜ = Z(i) (2.9) t t t t i=0 i=0 X X is expected to provide a reasonable approximation for the original model as long as the (i) residual term is small enough to allow the perturbative treatment. As we will see, V t (i) and Z , the corrections to each order can be calculated recursively using the results of t the lower order approximations. 1Or, one can consider ǫ=1 as simply a parameter convenient to count the approximation order. The actual quantity that should besmall for theapproximation is the residual part g. 4 2.3 Recursive Approximation for Perturbed linear FBSDE 2.3.1 Zero-th Order For the zero-th order of ǫ, one can easily see the following equation should be satisfied: (0) (0) (0) dV = c(X )V dt+Z dW (2.10) t t t t · t (0) V = Φ(X ) . (2.11) T T It can be integrated as Vt(0) = E e−RtTc(Xs)dsΦ(XT) Ft (2.12) h (cid:12) i which is equivalent to the pricing of a standard European contingent claim. (cid:12) (cid:12) Since we have T e−R0Tc(Xs)dsΦ(XT)= V0(0)+ e−R0uc(Xs)dsZu(0)·dWu (2.13) Z0 it can be shown that, by applying Malliavin derivative , t D T Dt e−R0Tc(Xs)dsΦ(XT) = Dt e−R0uc(Xs)dsZu(0) ·dWu+e−R0tc(Xs)dsZt(0) . (2.14) (cid:16) (cid:17) Zt (cid:16) (cid:17) Thus, by taking conditional expectation E[ ], we obtain t ·|F Zt(0) = E Dt e−RtTc(Xs)dsΦ(XT) Ft . (2.15) h (cid:16) (cid:17)(cid:12) i (cid:12) 2.3.2 First Order (cid:12) Now, let us consider the process V(ǫ) V(0). One can see that its dynamics is governed − by (ǫ) (0) (ǫ) (0) (ǫ) (ǫ) (ǫ) (0) d V V = c(X ) V V ǫg(X ,V ,Z )dt+ Z Z dW t − t t t − t − t t t t − t · t (ǫ) (0) (cid:0) V V (cid:1) = 0 . (cid:0) (cid:1) (cid:0) (cid:1) (2.16) T − T Now, by extracting the ǫ-first order terms, we can once again recover the linear FBSDE (1) (1) (0) (0) (1) dV = c(X )V dt g(X ,V ,Z )dt+Z dW (2.17) t t t − t t t t · t (1) V = 0 , (2.18) T which leads to T Vt(1) = E e−Rtuc(Xs)dsg(Xu,Vu(0),Zu(0))du Ft (2.19) (cid:20)Zt (cid:12) (cid:21) (cid:12) straightforwardly. By the same arguments in the zero-th order example, we can express (cid:12) (cid:12) the volatility term as T Zt(1) = E Dt e−Rtuc(Xs)dsg(Xu,Vu(0),Zu(0))du Ft . (2.20) (cid:20) (cid:18)Zt (cid:19)(cid:12) (cid:21) (cid:12) From these results, we can see that the required calculation is nothi(cid:12)ng more difficult than (cid:12) the zero-th order case as long as we have explicit expression for V(0) and Z(0). 5 2.3.3 Second and Higher Order Corrections We can proceed the same way to the second order correction. By extracting the ǫ-second (ǫ) (0) (1) order terms from V (V +ǫV ), one can show that t − t t ∂ (2) (2) (0) (0) (1) (0) (0) (1) (2) dV = c(X )V dt g(X ,V ,Z )V + g(X ,V ,Z ) Z dt+Z dW t t t − ∂v t t t t ∇z t t t · t t · t (cid:18) (cid:19) (2) V = 0 (2.21) T (2) is arelevantFBSDE, whichis onceagain linear in V . Asbefore, itleads tothefollowing t expression straightforwardly: Vt(2) = E T e−Rtuc(Xs)ds ∂∂vg(Xu,Vu(0),Zu(0))Vu(1) +∇zg(Xu,Vu(0),Zu(0))·Zu(1) du Ft (cid:20)Zt (cid:18) (cid:19) (cid:12) (cid:21) (cid:12) (2.22) (cid:12) (cid:12) Zt(2) = E Dt T e−Rtuc(Xs)ds ∂∂vg(Xu,Vu(0),Zu(0))Vu(1) +∇zg(Xu,Vu(0),Zu(0))·Zu(1) du Ft . (cid:20) (cid:18)Zt (cid:16) (cid:17) (cid:19)(cid:12)(cid:12)(2.2(cid:21)3) (cid:12) (cid:12) In the above calculation, we have assumed the driver function is differentiable. If this is not the case, we need to approximate it using some smooth function or apply integration- by-partstechniqueforgeneralized Wienerfunctionals(e.g. acompositefunctionalofDirac delta fucntion and a smooth Wiener functional). In exactly the same way, one can derive an arbitrarily higher order correction. Due to the ǫ in front of the non-linear term g, the system remains to be linear in the every order of approximation. However, in order to carry out explicit evaluation, we need to give Malliavin derivative explicitly in terms of the forward components. We will discuss this issue in the next. 2.4 Evaluation of Malliavin Derivative Firstly, let us introduce a d d matrix process Y , for u [t,T], as the solution for the t,u × ∈ following forward SDE: d d(Y )i = ∂ γi(X )(Y )kdu+∂ γi(X )(Y )k dW (2.24) t,u j k 0 u t,u j k u t,u j · u Xk=1(cid:16) (cid:17) (Y )i = δi , (2.25) t,t j j where ∂ denotes the differential with respect to the k-th component of X, and δi denotes k j Kronecker delta. Now, for Malliavin derivative, we want to express, for u [t,T], ∈ u E t e−Rt c(Xs)dsG(Xu) t (2.26) D F h (cid:16) (cid:17)(cid:12) i (cid:12) (cid:12) 6 in terms of X , where G is a some deterministic function of X, in general. Thank to the t known chain rule of Malliavin derivative, we have d u u t e−Rt c(Xs)dsG(Xu) = e−Rt c(Xs)ds∂iG(Xu)( tXui) D D (cid:16) (cid:17) Xi=1n u u e−Rt c(Xs)dsG(Xu) ∂ic(Xs)( tXsi)ds .(2.27) − D (cid:18)Zt (cid:19)o Thus, it is enough for our purpose to evaluate ( X ). Since we have t u D d u u ( Xi) = ∂ γi(X )( Xk)ds+ ∂ γi(X )( Xk) dW +γi(X ) (2.28) Dt u k 0 s Dt s k s Dt s · s t Xk=1(cid:16)Zt Zt (cid:17) it can be shown that X follows the next SDE: t u D d d( Xi) = ∂ γi(X )( Xk)du+∂ γi(X )( Xk) dW (2.29) Dt u k 0 u Dt u k u Dt u · u Xk=1(cid:16) (cid:17) ( Xi) = γi(X ) . (2.30) t t t D Thus, comparing to Eqs. (2.24) and (2.25), we can conclude that ( X ) = Y γ(X ) . (2.31) t u t,u t D Asaresult, combiningtheSDEfor Y andtheMarkovian propertyof X, onecanconfirm t,u that the conditional expectation u E t e−Rt c(Xs)dsG(Xu) t (2.32) D F h (cid:16) (cid:17)(cid:12) i is actually given by a some function of X . Therefore, in(cid:12) principle, both of the backward t (cid:12) components can be expressed in terms of X in each approximation order. t In fact, this is an easy task when the underlying process has a known distribution. In the next section, we present two such models, and demonstrate how our approximation scheme works. We will also compare our approximate solution to the direct numerical resultsobtainedfrom,suchasPDEandMonteCarlosimulation. However, inmoregeneric situations, we do not know the distribution of X. We will explain how to handle the problem in this case using the asymptotic expansion method for the forward components in Sec.4. 3 Simple Examples 3.1 A forward agreement with bilateral default risk As the first example, we consider a toy model for a forward agreement on a stock with bilateral default risk of the contracting parties, the investor (party-1) and its counterparty (party-2). The terminal payoff of the contract from the view point of the party-1 is Φ(S )= S K (3.1) T T − 7 where T is the maturity of the contract, and K is a constant. We assume the underlying stock follows a simple geometric Brownian motion: dS = rS dt+σS dW (3.2) t t t t wheretherisk-freeinterestrater andthevolatility σ areassumedtobepositiveconstants. The default intensity of party-i h is specified as i h = λ, h = λ+h (3.3) 1 2 where λ and h are also positive constants. In this setup, the pre-default value of the contract at time t, V , follows 2 t dV = rV dt h max( V ,0)dt+h max(V ,0)dt+Z dW t t 1 t 2 t t t − − = (r+λ)V dt+hmax(V ,0)dt+Z dW (3.4) t t t t V = Φ(S ) . (3.5) T T Now, following the previous arguments, let us introduce the expansion parameter ǫ, and consider the following FBSDE: (ǫ) (ǫ) (ǫ) (ǫ) dV = µV dt ǫg(V )dt+Z dW (3.6) t t − t t t (ǫ) V = Φ(S ) (3.7) T T dS = S (rdt+σdW ) , (3.8) t t t where we have defined µ = r+λ and g(v) = hv1 . v 0 − { ≥ } 3.1.1 Zero-th order In the zero-th order, we have (0) (0) (0) dV = µV dt+Z dW (3.9) t t t t (0) V = Φ(S ) . (3.10) T T Hence we simply obtain V(0) = E e µ(T t)Φ(S ) t − − T Ft h (cid:12) i = e µ(T t) S er(T t(cid:12)) K (3.11) − − t −(cid:12)− (cid:16) (cid:17) and Z(0) = e λ(T t)σS . (3.12) t − − t 2See, for example, [3, 5]. 8 3.1.2 First order In the first order, we have (1) (1) (0) (1) dV = µV dt g(V )dt+Z dW (3.13) t t − t t t (1) V = 0 . (3.14) T Thus, we obtain T V(1) = E e µ(T u)g(V(0))du (3.15) t − − u Ft (cid:20)Zt (cid:12) (cid:21) T (cid:12) = e−µ(T−t)h E max(Su(cid:12)(cid:12)er(T−u) K,0) t du (3.16) − − F Zt h (cid:12) i T (cid:12) = e µ(T t)h C(u;t,S )du , (cid:12) (3.17) − − t − Zt where C(u;t,S ) = S er(T t)N d (u;t,S ) KN d (u;t,S ) (3.18) t t − 1 t 2 t − d (u;t,S ) = 1 l(cid:0)n Ster(T−t(cid:1)) 1σ(cid:0)2(u t) ,(cid:1) (3.19) 1(2) t σ√u t K ± 2 − ! − (cid:16) (cid:17) and N denotes the cumulative distribution function for the standard normal distribution. We can also derive T Z(1) = e λ(T t)hσS N(d (u;t,S ))du . (3.20) t − − − t 1 t Zt 3.1.3 Second order Finally, let us consider the second order value adjustment. In this case, the relevant dynamics is given by ∂ (2) (2) (0) (1) (2) dV = µV dt g(V )V dt+Z dW (3.21) t t − ∂v t t t t (2) V = 0 . (3.22) T As a result, we have T ∂ V(2) = E e µ(u t) g(V(0))V(1) du (3.23) t − − ∂v u u Ft (cid:20)Zt (cid:18) (cid:19) (cid:12) (cid:21) T T (cid:12) = e−µ(T−t)h2 E 1 Suer(T−u) K (cid:12)(cid:12)0 C(s;u,Su) Ft dsdu (3.24) Zt Zu h { − ≥ } (cid:12) i (cid:12) which can be evaluated as (cid:12) T T V(2) = e µ(T t)h2 ∞ φ(z)C s;u,S (z,S ) dzdsdu , (3.25) t − − u t Zt Zu Z−d2(u;t,St) (cid:0) (cid:1) 9 where we have defined Su(z,St) = Ste(r−12σ2)(u−t)+σ√u−tz (3.26) and 1 φ(z) = e−12z2 . (3.27) √2π 3.1.4 Numerical comparison to PDE Forthissimplemodel,wecandirectlyevaluatethecontractvalueV bynumericallysolving t the PDE: ∂ ∂ 1 ∂2 V(t,S)+ rS V(t,S)+ σ2S2 V(t,S) µ+h1 V(t,S) = 0 ∂v ∂s 2 ∂s2 − {V(t,S)≥0} (cid:18) (cid:19) h i with the boundary conditions V(T,S) = S K − V(t,M) = e (µ+h)(T t) Mer(T t) K , M K − − − − ≫ V(t,m) = e−µ(T−t) me(cid:0)r(T−t) K , (cid:1) m K . − ≪ (cid:0) (cid:1) Figure 1: Numerical Comparison to PDE In Fig. 1, we have plot the numerical results of the forward contract with bilateral default risk with various maturities with the direct solution from the PDE. We have used r = 0.02, λ =0.01, h = 0.03, (3.28) σ = 0.2, S = 100 , (3.29) 0 where the strike K is chosen to make V(0) = 0 for each maturity. We have plot V(1) 0 for the first order, and V(1) +V(2) for the second order. Note that we have put ǫ = 1 to comparetheoriginalmodel. Onecanobservehowthehigherordercorrectionimprovesthe 10