A Polynomial Scheme of Asymptotic Expansion for Backward SDEs and Option pricing
aa r X i v : . [ q -f i n . C P ] D ec A Polynomial Scheme of Asymptotic Expansionfor Backward SDEs and Option pricing ∗ Masaaki Fujii † First version: May 2, 2014This version: December 22, 2014
Abstract
A new asymptotic expansion scheme for backward SDEs (BSDEs) is proposed. Theperturbation parameter “ ǫ ” is introduced just to scale the forward stochastic variableswithin a BSDE. In contrast to the standard small-diffusion asymptotic expansion method,the dynamics of variables given by the forward SDEs is treated exactly. Although itrequires a special form of the quadratic covariation terms of the continuous part, it allowsrather generic drift as well as jump components to exist. The resultant approximationis given by a polynomial function in terms of the unperturbed forward variables whosecoefficients are uniquely specified by the solution of the recursive system of linear ODEs.Applications to a jump-extended Heston and λ -SABR models for European contingentclaims, as well as the utility-optimization problem in the presence of a terminal liabilityare discussed. Keywords :
Stochastic Control, Asymptotic Expansion, BSDE, random measure, Heston,SABR, Utility optimization ∗ To appear in
Quantitative Finance . All the contents expressed in this research are solely those of theauthors and do not represent any views or opinions of any institutions. The authors are not responsible orliable in any manner for any losses and/or damages caused by the use of any contents in this research. † Graduate School of Economics, The University of Tokyo. e-mail: [email protected] Introduction
The backward stochastic differential equation (BSDE) was introduced by Bismut (1973) [2]under a linear setup and was later extended by Pardoux & Peng (1990) [29] to general non-linear situations. Although the research activity had been contained in a relatively smallmathematical community, it has been rapidly gaining traction with financial researchers andpractitioners, in particular, since the last financial crisis. This is because that one almostinevitably encounters BSDEs when he/she tries to handle various non-linear effects arisingfrom credit risk, collateralization, funding and regulatory costs, and various other sourcesof incompleteness arising from the new market realities. See, for example, Duffie & Huang(1996) [14], Fujii & Takahashi (2013) [16], Cr´epey (2013) [9] and a summary of recent practicaltopics in the financial industry [3, 5]. Various interesting financial applications, such as forinsurance, utility indifference pricing and optimal contract theory can be found in books[12, 7, 11, 10]. One can consult with [15, 25] as a text for general mathematical treatmentsof BSDEs.There now exists vast literature on their numerical treatments, ranging from the famousfour-step scheme proposed by Ma et al. (1994) [26], its discretized implementation by Douglaset al. (1996) [13], various Monte-Carlo techniques making use of the least-square regressionmethod (See, for example, Bouchard & Touzi (2004) [6], Bender & Denk (2007) [1], Go-bet et al. (2005) [21] and Gobet & Lemor (2010) [22]), a branching diffusion method byHenry-Labordere (2012) [23], and a particle method by Fujii & Takahashi (2012a) [17]. Un-fortunately though, many of them require a good amount of experience and deep expertiseto achieve stable results, such as an appropriate choice of basis functions, the order of regres-sions, and of course, a good programming technique.It is obvious that a simple analytical approximation method is deeply wanted. In Fujii& Takahashi (2012b) [18], we have developed a driver perturbation method combined witha standard asymptotic expansion technique for the forward SDEs. Its error estimate wasrecently provided by Takahashi & Yamada (2013) [33]. It is systematic and straightforward,but one still needs to endure long tough calculations especially for higher order corrections,which is stemming from the needs of evaluation of conditional expectations at each orderof expansion. An interesting exceptional case arises if a so-called quadratic-growth BSDE(qgBSDE) is associated with linear Gaussian forward SDEs, and at the same time, if itsterminal value is given by, at most, quadratic form of the Gaussian variables. In this case,the value function is given by a quadratic function of the Gaussian variables whose coefficientsare completely determined by the ordinary differential equations (ODEs) involving ones withRiccati form. See, for example, Schroder & Skiadas (1999) [31] as an early research. Recently,this property was applied to the mean-variance (quadratic) hedging problem by Fujii &Takahashi (2013, 2014) [19, 20], making use of the beautiful BSDE expression derived byMania & Tevzadze (2003) [27]. Notice that the Riccati equation may possibly diverge in afinite time-interval in a general setup. In such a case, one needs to shorten the maturity ofthe corresponding problem.In this paper, we propose a new scheme which approximates a solution of a BSDE by apolynomial function of the underlying variables. In a Markovian setup, it is well-known thatthe solution of a BSDE is given by a Markovian function of the underlying variables [26].Therefore, it is intuitively clear that the solution should be well approximated by a poly-2omial function for short maturities within which the size of the underlying variables (afterappropriate rescaling and shift of their means) remains relatively small. Despite the appar-ent similarity to the usual asymptotic expansion, the new scheme yields a recursive systemof linear
ODEs which can be obtained by simply matching the coefficients of the assumedpolynomial solution to those of the BSDE’s driver. Although we have to assume that theforward processes have a special form of quadratic covariation of the continuous part, theycan have rather general drifts and random jump components. In that sense, the method canbe interpreted as a generalizations of exact but exceptional quadratic-solution example to an approximate polynomial-solution technique with wider applications.The organization of the paper is as follows: In Section 2, the main idea of the polynomialexpansion scheme is explained. In order to show its usefulness and accuracy, we apply themethod to several well-known problems in the remaining part of the paper. In Section 3and 4, the proposed scheme is applied to European contingent claims using a jump-extendedHeston model and the λ -SABR model. We provide the closed expression for the recursivesystem of linear ODEs which specifies the approximate solution at an arbitrary order. InSection 5, the optimization problem for the exponential utility in the presence of a terminalliability is analyzed. The closed-form system of the linear ODEs is derived for this setup, too.Each model is associated with several illustrative numerical examples. Finally in Appendix,some details omitted in the main text are provided. Let us consider the following system of forward and backward SDEs in a probability space(Ω , F , P ): V t = H ( X T ) − Z Tt ¯ f (cid:16) s, X s , V s , ¯ Z s , Z K U ( s, z ) Q ( s, dz ) (cid:17) ds − Z Tt ¯ Z s dW s − Z Tt Z K U ( s, z ) e N ( ds, dz ) (2.1) X t = x + Z t b ( s, X s ) ds + Z t σ ( s, X s ) dW s + Z t Z K z N ( ds, dz ) (2.2)where x ∈ R is a constant, W one-dimensional standard Brownian motion and N a randommeasure whose deterministic jump distribution is given by Q ( t, · ) with some (compact) space K for its support. e N ( dt, dz ) is the corresponding P -compensated random measure e N ( dt, dz ) := N ( dt, dz ) − λ ( t, X t ) Q ( t, dz ) dt (2.3)where λ ( t, X t ) denotes the jump intensity. We assume that H : R → R , ¯ f : [0 , T ] × R → R , b : [0 , T ] × R → R and λ : [0 , T ] × R → R + are all smooth functions. In addition, we assumethat the quadratic covariation of X from its continuous part is given by, at most, quadratic3orm of X itself d h X c i t = (cid:16) σ ( t ) X t + σ ( t ) X t + σ ( t ) (cid:17) dt (2.4)where the superscript “ c ” denotes the continuous part of X . Here, ( σ i ( t )) i =1 , , is the setof deterministic functions in such a way that it guarantees the right-hand side of (2.4) isnon-negative for every possible value taken by X . We assume that the forward-backwardSDE system of (2.1) and (2.2) has a well-posed solution.Since the stochasticity of ( V t ) t ≥ is provided solely by ( X t ) t ≥ , we can rewrite the BSDEas V t = H ( X T ) − Z Tt f (cid:16) s, X s , V s , Z s , Z K U ( s, z ) Q ( s, dz ) (cid:17) ds − Z Tt Z s dX cs − Z Tt Z K U ( s, z ) N ( ds, dz ) (2.5)with appropriate redefinition of f ( · ) and Z . Here, dX ct := b ( t, X t ) dt + σ ( t, X t ) dW t denotesthe continuous part of the X ’s change. Thus, in the following, we consider the equivalentsystem given by (2 .
5) and (2 . In order to obtain a polynomial expansion, we introduce ǫ so that we can count the order ofthe underlying X . Let consider the following perturbed BSDE: V ǫt = H ( ǫX T ) − Z Tt f (cid:16) s, ǫX s , V ǫs , Z ǫs , Z K U ǫ ( s, z ) Q ( s, dz ) (cid:17) ds − Z Tt Z ǫs dX cs − Z Tt Z K U ǫ ( s, z ) N ( ds, dz ) . (2.6)Here, the superscripts ǫ in V, Z, U emphasizes that these variables are now dependent on theparameter ǫ . An important difference from the usual small diffusion asymptotic expansionmethod proposed by (Yoshida (1992a) [35], Takahashi (1999) [32], Kunitomo & Takahashi(2003) [24] for the pricing of contingent claims, Yoshida (1992b) [36] for statistical applications) based on Watanabe (1987) [34] theory is that the underlying process ( X t ) t ≥ itself is notperturbed and only its size is scaled by ǫ within the BSDE.We assume that that the expansion V ǫt = ∞ X n =0 ǫ n V [ n ] t (2.7) V [ n ] t = n X m =0 X mt m ! v [ n ] m ( t ) (2.8) For example, Z and ¯ Z is connected by the relation ¯ Z s = Z s σ ( s, X s ).
4s well defined, where ( v [ n ] m ( t ) , ≤ m ≤ n ) are all deterministic bounded functions in a giventime interval t ∈ [0 , T ]. In particular, the Itˆo-formula should be applicable to (2.8) so thatone obtains the well-defined forward SDE for ( V [ n ] t ) ≤ t ≤ T . It leads to the correspondingexpansions of the control variables Z ǫt = ∞ X n =0 ǫ n Z [ n ] t (2.9) U ǫ ( t, z ) = ∞ X n =0 ǫ n U [ n ] ( t, z ) (2.10)whose expressions can be easily derived once (2.8) is obtained.If the maturity is short enough so that size of X remains small, truncating the expansionin (2.7) at a certain order n and putting ( ǫ = 1) are expected to give an approximation of theoriginal problem. Note that, since ǫ is introduced as a combination ( ǫX ), discussing the sizeof ǫ separately from X is not useful. Let us denote the truncated n -th order approximationas (using the superscript ( n ) instead of [ n ]) V ( n ) t = n X j =0 V [ j ] t (2.11) Z ( n ) t = n X j =0 Z [ j ] t , U ( n ) ( t, z ) = n X j =0 U [ j ] ( t, z ) . (2.12)Note that all of them are given by the polynomial functions of X , at most, of the order of n .One can check the accuracy of approximation by comparing e V ( n ) ( T ) = V ( n )0 + Z T f (cid:16) s, X s , V ( n ) s , Z ( n ) s , Z K U ( n ) ( s, z ) Q ( s, dz ) (cid:17) ds + Z T Z ( n ) s dX cs + Z T Z K U ( n ) ( s, z ) N ( ds, dz ) (2.13)and H ( X T ) in a “ path-wise ” fashion. In numerical examples given in later sections, weshall observe that the polynomial approximation gives a good path-wise approximation, atleast for the paths along which ( | X t ( ω ) | ) t ≥ does not significantly grow to a big value. Inpractical applications, the capability of checking ( H ( X T ) − e V ( n ) ( T )) directly should be agreat help for setting aside an appropriate amount of risk-reserve for the hedging program tobe implemented. By the very nature of polynomial approximation, one can imagine that ahigher order expansion may yield an unstable result in a very volatile market, or for a problemwith long maturity. The above comparison gives useful information for an appropriate orderof expansion for a given situation.As we shall see below, all the functions ( v [ n ] m ( t )) m,n except v [0]0 ( t ) are specified by linear ODEs. In later sections which deal with specific models, we provide a closed form recursivesystem of linear ODEs which fixes the coefficients of the polynomials up to an arbitrary order.However, in this section, let us adopt a slightly tedious step-by-step explanation, which wehope to give a clearer image to the readers. 5 .2.2 Zero-th order
It is obvious that V [0] t = H (0) − Z Tt f ( s ) ds (2.14)where f ( s ) := f ( s, , V [0] s , , v [0]0 ( t ) = f ( t, , v [0]0 ( t ) , , , v [0]0 ( T ) = H (0) . (2.15)This is the only non-linear ODE we encounter. We assume that the finite solution exist to therelevant time interval t ∈ [0 , T ]. This should be the case for most of the natural applications,since the 0-th order problem corresponds to the market where all the underlyings are constant. Thanks to the smoothness assumption, one has V [1] t = ∂ x H (0) X T − Z Tt Z [1] s dX cs − Z Tt Z K U [1] ( s, z ) N ( ds, dz ) − Z Tt ( ∂ x f ( s ) X s + ∂ v f ( s ) V [1] s + ∂ z f ( s ) Z [1] s + ∂ u f ( s ) Z K U [1] ( s, z ) Q ( s, dz ) ) ds . (2.16)On the other hand, let us suppose the above solution is given by V [1] t = v [1]1 ( t ) X t + v [1]0 ( t ) . (2.17)Then, it yields the dynamics dV [1] t = n ˙ v [1]1 ( t ) X t + ˙ v [1]0 ( t ) o dt + v [1]1 ( t ) dX ct + v [1]1 ( t ) Z K z N ( dt, dz ) . (2.18)By comparing (2.16) and (2.18), we should have Z [1] t = v [1]1 ( t ) (2.19) U [1] ( t, z ) = v [1]1 ( t ) z . (2.20)Substituting the expressions of V [1] , Z [1] and U [1] into (2.16) and matching its driver tothe drift part of (2.18), one obtains˙ v [1]1 ( t ) = ∂ x f ( t ) + ∂ v f ( t ) v [1]1 ( t ) (2.21)˙ v [1]0 ( t ) = ∂ v f ( t ) v [1]0 ( t ) + (cid:16) ∂ z f ( t ) + ∂ u f ( t ) q ( t, (cid:17) v [1]1 ( t ) (2.22)6ith terminal conditions v [1]1 ( T ) = ∂ x H (0) and v [1]0 ( T ) = 0. Here, we have used the notation q ( t, n ) = Z K z n Q ( t, dz ) (2.23)to denote the n-th jump moment. In the reminder of the paper, we assume the existence of themoments relevant for the approximation scheme. It is clear that the assumed solution (2.17)and the corresponding control variables with the coefficients satisfying the above ODEs is infact one solution of the BSDE (2.16) as long as the forward SDE (2.18) is well-defined. Dueto the linearity of the ODEs, the solution should be unique among the assumed polynomialforms. For the 2nd and higher order corrections, the assumption on the quadratic covariation termplays an important role. As before, let us suppose that the solution takes the polynomialform: V [2] t = v [2]2 ( t ) X t
2! + v [2]1 ( t ) X t + v [2]0 ( t ) . (2.24)Then, a simple application of Itˆo-formula yields dV [2] t = (cid:18) ˙ v [2]2 ( t ) X t
2! + ˙ v [2]1 ( t ) X t + ˙ v [2]0 ( t ) (cid:19) dt + 12 v [2]2 ( t ) d h X c i t + (cid:16) v [2]2 ( t ) X t + v [2]1 ( t ) (cid:17) dX ct + Z K (cid:18) v [2]2 ( t ) ( X t − + z ) − X t − v [2]1 ( t ) z (cid:19) N ( dt, dz ) . (2.25)It should be clear that the assumption made in (2 .
4) is necessary to guarantee that the highest polynomial order assumed in V [ n ] t remains n under the dynamics of ( X t ) t ≥ . Theexpression in (2.25) now implies Z [2] t = v [2]2 ( t ) X t + v [2]1 ( t ) (2.26) U [2] ( t, z ) = v [2]2 ( t ) (cid:16) X t − z + z (cid:17) + v [2]1 ( t ) z . (2.27)7n the other hand, the 2nd order part of (2.6) leads to a BSDE V [2] t = X T ∂ x H (0) − Z Tt Z [2] s dX cs − Z Tt Z K U [2] ( s, z ) N ( ds, dz ) − Z Tt ( ∂ v f ( s ) V [2] s + ∂ z f ( s ) Z [2] s + ∂ u f ( s ) (cid:16)Z K U [2] ( s, z ) Q ( s, dz ) (cid:17) + 12 ∂ x f ( s ) X s + 12 ∂ v f ( s )[ V [1] s ] + 12 ∂ z f ( s )[ Z [1] s ] + 12 ∂ u f ( s ) (cid:16)Z K U [1] ( s, z ) Q ( s, dz ) (cid:17) + X s (cid:18) ∂ x,v f ( s ) V [1] s + ∂ x,z f ( s ) Z [1] s + ∂ x,u f ( s ) (cid:16)Z K U [1] ( s, z ) Q ( s, dz ) (cid:17)(cid:19) + V [1] s (cid:18) ∂ v,z f ( s ) Z [1] s + ∂ v,u f ( s ) (cid:16)Z K U [1] ( s, z ) Q ( s, dz ) (cid:17)(cid:19) + ∂ z,u f ( s ) Z [1] s Z K U [1] ( s, z ) Q ( s, dz ) ) (2.28)Although there appear many cross terms, the same procedures as in the first-order case ofmatching the coefficients of X in the driver of (2.28) to those in (2.25) give us a set of linearODEs.After a simple calculation, one can confirm that the relevant ODEs are given by˙ v [2]2 ( t ) = (cid:16) ∂ v f ( t ) − σ ( t ) (cid:17) v [2]2 ( t ) + ∂ v f ( t ) (cid:2) v [1]1 ( t ) (cid:3) + 2 ∂ x,y f ( t ) v [1]1 ( t ) + ∂ x f ( t ) (2.29)˙ v [2]1 ( t ) = ∂ v f ( t ) v [2]1 ( t ) + (cid:16) ∂ z f ( t ) + ∂ u f ( t ) q ( t, − σ ( t )2 (cid:17) v [2]2 ( t )+ (cid:16) ∂ v,z f ( t ) + ∂ v,u f ( t ) q ( t, (cid:17) [ v [1]1 ( t )] + ∂ v f ( t ) v [1]1 ( t ) v [1]0 ( t )+ (cid:16) ∂ x,z f ( t ) + ∂ x,u f ( t ) q ( t, (cid:17) v [1]1 ( t ) + ∂ x,v f ( t ) v [1]0 ( t ) (2.30)˙ v [2]0 ( t ) = ∂ v f ( t ) v [2]0 ( t ) + 12 (cid:16) ∂ u f ( t ) q ( t, − σ ( t ) (cid:17) v [2]2 ( t )+ (cid:16) ∂ z f ( t ) + ∂ u f ( t ) q ( t, (cid:17) v [2]1 ( t ) + 12 ∂ v f ( t )[ v [1]0 ( t )] + 12 (cid:16) ∂ z f ( t ) + ∂ u f ( t ) q ( t, + 2 ∂ z,u f ( t ) q ( t, (cid:17) [ v [1]1 ( t )] + (cid:16) ∂ v,z f ( t ) + ∂ v,u f ( t ) q ( t, (cid:17) v [1]1 ( t ) v [1]0 ( t ) (2.31)with terminal conditions v [2]2 ( T ) = ∂ x H (0) , v [2]1 ( T ) = v [2]0 ( T ) = 0 . (2.32)Given the solution for the 1st-order expansion ( v [1]1 ( t ) , v [1]0 ( t )), the above ODEs can be solved one-by-one following the order of v [2]2 → v [2]1 → v [2]0 . If the forward dynamics (2.25) of the8ypothesized solution is well-defined, then it is clear that it actually gives one solution forthe BSDE of the 2nd order (2.28). Due to the linearity of the ODE, the solution should beunique among the assumed forms.It is clear that one can repeat the procedures up to an arbitrary order. At any order n ( ≥ n -th order. Remark:
In the above example, the distribution ( Q ( t, · )) t ≥ is not necessary be deterministic. q ( t, n )can be a polynomial function of X at most of the order of n . However, it is important to notethat this point is dependent on how the jump component is introduced in the model: If X has a proportional jump, then q ( t, n ) specifying the n -th moment of the proportional jumpfactor should be independent of X . For example, one can consider the conditions to keep(2 .
28) as a 2nd-order polynomial of X . Unfortunately, we have not yet obtained a good understanding of the mathematical prop-erties of the proposed expansion. Despite the similarity to Takahashi [32] and Kunitomo &Takahashi [24] in the way that the parameter ǫ is introduced and its similar application toBSDEs in Fujii & Takahashi [18], it is not yet clear if we can simply borrow the argumentsin Takahashi & Yamada [33] for justification to the current polynomial scheme. A rigorousproof is left for an important topic for the future research.However, we would like to emphasize that the above limitation is not a significant draw-back for practical applications. The great advantage to have an explicit form of an approxi-mate solution is allowing one to test its accuracy directly for a given setup (See the discussionin Section 2.2.1.). The test like this is necessary for any methods since the convenient as-sumptions needed for the mathematical justification will be violated in realistic situationsanyway. In contrast to the proposed scheme (and the one in [18]), one can see that carryingout this check is not a simple task for purely simulation-based techniques.In addition, it is interesting to notice that we have not used any special properties ofWiener integral after expressing the BSDE in term of dX as in (2.5). If one can loosen theconditions necessary for the quadratic covariation terms, one may possibly obtain an unifiedway of approximation for rather general semimartingales.9 Pricing European Optionswith a Jump-extended Heston Model
We assume that the asset price S and its stochastic variance-factor Y have the followingdynamics under a probability space (Ω , F , Q ): S t = S + Z t S s (cid:18) σ ( s ) p ¯ Y s dW s + Z K ( e z − e N ( ds, dz ) (cid:19) (3.1)¯ Y t = 1 + Z t (cid:16) α ( s ) p ¯ Y s dB s + κ ( s )(1 − ¯ Y s ) ds (cid:17) (3.2)where Q is supposed to be a certain equivalent martingale measure chosen by market partic-ipants. W and B denote one-dimensional Q -Brownian motions with d h W, B i t = ρ ( t ) dt . e N denotes Q -compensated random measure specified by e N ( dt, dz ) = N ( dt, dz ) − ¯ λ ( t, ¯ Y t ) Q ( t, dz ) dt (3.3)with the jump intensity ¯ λ and its deterministic distribution function Q ( t, · ). σ ( · ) , α ( · ) , ρ ( · )and κ ( · ) are appropriate deterministic functions. We allow ¯ λ : [0 , T ] × R + → R + to be asmooth generic function of ¯ Y , and hence (3.1) and (3.2) do not consist of the analyticallysolvable affine system.In order to make the expansion around the origin a good approximation, we perform achange of variables X t := ln (cid:18) S t S (cid:19) Y t := ¯ Y t − . (3.4)Then, they follow the dynamics X t = Z t (cid:18) σ ( s ) p Y s + 1 dW s − h σ ( s ) Y s + 1) + λ ( s, Y s ) β ( s ) i ds (cid:19) + Z t Z K z N ( ds, dz ) (3.5) Y t = Z t (cid:16) α ( s ) p Y s + 1 dB s − κ ( s ) Y s ds (cid:17) (3.6)where β ( · ) is a deterministic function defined by β ( t ) = Z K ( e z − Q ( t, dz ) (3.7)and λ ( t, Y t ) := ¯ λ ( t, Y t + 1). 10et us consider the valuation problem for a European option in a BSDE form: V t = H ( X T ) − Z Tt ¯ Z s dW s − Z Tt ¯Γ s dB s − Z Tt Z K U ( s, z ) e N ( ds, dz ) (3.8)where H ( X T ) denotes the terminal payoff at the maturity T in terms of X , and V t denotes itspresent value at time t ( < T ). Simple redefinition of the control variables ( ¯ Z, ¯Γ), one obtains V t = H ( X T ) − Z Tt Z s dX cs − Z Tt Γ s dY s − Z Tt Z K U ( s, z ) N ( ds, dz ) − Z Tt (cid:26) Z s h σ ( s ) Y s + 1) + λ ( s, Y s ) β ( s ) i + κ ( s )Γ s Y s − λ ( s, Y s ) Z K U ( s, z ) Q ( s, dz ) (cid:27) ds . One can now apply the proposed polynomial expansion scheme to the BSDE if H ( · ) is asmooth function. Although we can directly approximate the option payoff by a polynomialfunction, we shall take an alternative road that does not involve such approximation. We aregoing to consider H ( x ) = x m for m = 1 , , , · · · . Then the corresponding value function V t gives the moments of X T . We finally use the Edgeworth expansion to get an estimate of theprobability density function of X T (and hence S T ) to calculate the standard Call and Putoptions. We consider the system of a perturbed BSDE V ǫt = H ( ǫX T ) − Z Tt Z ǫs dX cs − Z Tt Γ ǫs dY s − Z Tt Z K U ǫ ( s, z ) N ( ds, dz ) − Z Tt (cid:26) Z ǫs h σ ( s ) ǫY s + 1) + λ ( s, ǫY s ) β ( s ) i + ǫκ ( s )Γ ǫs Y s − λ ( s, ǫY s ) Z K U ǫ ( s, z ) Q ( s, dz ) (cid:27) ds (3.9)and the forward SDEs (3.5) and (3.6). We expand the solution in term of ǫ as V ǫt = ∞ X n =0 ǫ n V [ n ] t (3.10) Z ǫt = ∞ X n =0 ǫ n Z [ n ] t , Γ ǫt = ∞ X n =0 ǫ n Γ [ n ] t , U ǫ ( t, z ) = ∞ X n =0 ǫ n U [ n ] ( t, z ) . (3.11)In the next lemma, we give the solution of the above expansion in terms of a recursivesystem of linear ODEs. We denote the number of choices selecting m out of n ( ≥ m ) by C ( n,m ) = n !( n − m )! m ! . We also use the convention for the summation symbol that P ji ≡ j < i ). Lemma 1
If it exists, the polynomial solution for the expansion in (3.10) and (3.11) is niquely given by V [ n ] t = n X m =0 m X k =0 X m − kt Y kt ( m − k )! k ! v [ n ] m − k,k ( t ) (3.12) Z [ n ] t = n X m =1 m − X k =0 X m − k − t Y kt ( m − k − k ! v [ n ] m − k,k ( t ) (3.13)Γ [ n ] t = n X m =1 m X k =1 X m − kt Y k − t ( m − k )!( k − v [ n ] m − k,k ( t ) (3.14) U [ n ] ( t, z ) = n − X m =0 m X k =0 X m − kt Y kt ( m − k )! k ! (cid:16) n X l = m +1 z l − m ( l − m )! v [ n ] l − k,k ( t ) (cid:17) (3.15) with the set of deterministic functions v [ n ] m − k,k ( t ) of (0 ≤ k ≤ m ≤ n ) satisfying the followingrecursive system of linear ODEs ˙ v [ n ] m − k,k ( t ) = − I ( m ≤ n − , ≤ k ) k (cid:18) σ ( t ) v [ n ] m − k +2 ,k − ( t ) + ρ ( t ) σ ( t ) α ( t ) v [ n ] m − k +1 ,k ( t ) + α ( t ) v [ n ] m − k,k +1 ( t ) (cid:19) − I ( m ≤ n − (cid:18) σ ( t )2 v [ n ] m − k +2 ,k ( t ) + ρ ( t ) σ ( t ) α ( t ) v [ n ] m − k +1 ,k +1 ( t ) + α ( t ) v [ n ] m − k,k +2 ( t ) (cid:19) + I ( m ≤ n − , ≤ k ) k (cid:18) σ ( t ) v [ n − m − k +1 ,k − ( t ) + κ ( t ) v [ n − m − k,k ( t ) (cid:19) + I ( m ≤ n − σ ( t ) v [ n ] m − k +1 ,k ( t )+ I ( m ≤ n − k X l =0 C ( k,l ) ∂ ly λ ( t, β ( t ) v [ n − l ] m − k +1 ,k − l ( t ) − n − m X j =1 v [ n − l ] j + m − k,k − l ( t ) q ( t, j ) j ! (3.16) having the terminal conditions v [ n ] n, ( T ) = ∂ nx H (0) with all the other components zero.Proof: Let us suppose that the polynomial-form solution given in (3.12) exists. Then, theapplication of Itˆo-formula and simple rearrangements of summation yield dV [ n ] t = n X m =0 m X k =0 X m − kt Y kt ( m − k )! k ! ( ˙ v [ n ] m − k,k ( t )+ I ( m ≤ n − , ≤ k ) k (cid:18) σ t v [ n ] m − k +2 ,k − ( t ) + ρ t σ t α t v [ n ] m − k +1 ,k ( t ) + α t v [ n ] m − k,k +1 ( t ) (cid:19) + I ( m ≤ n − (cid:18) σ t v [ n ] m − k +2 ,k ( t ) + ρ t σ t α t v [ n ] m − k +1 ,k +1 ( t ) + α t v [ n ] m − k,k +2 ( t ) (cid:19) ) dt + n X m =1 m − X k =0 v [ n ] m − k,k ( t ) X m − k − t Y kt ( m − k − k ! dX ct + n X m =1 m X k =1 v [ n ] m − k,k ( t ) X m − kt Y k − t ( m − k )!( k − dY t + n − X m =0 m X k =0 X m − kt Y kt ( m − k )! k ! (cid:16) n X l = m +1 v [ n ] l − k,k ( t ) Z K z l − m ( l − m )! N ( dt, dz ) (cid:17) (3.17)12hich then implies (3.13), (3.14) and (3.15).On the other hand, extracting the n -th order part from the BSDE (3.9), one obtains V [ n ] t = X nT n ! ∂ nx H (0) − Z Tt Z [ n ] s dX cs − Z Tt Γ [ n ] s dY s − Z Tt Z K U [ n ] ( s, z ) N ( ds, dz ) − Z Tt ( σ s (cid:16) Y s Z [ n − s + Z [ n ] s (cid:17) + β ( s ) n − X l =0 ∂ ly λ ( s, l ! Y ls Z [ n − l ] s + κ s Y s Γ [ n − s − n − X l =0 ∂ ly λ ( s, l ! Y ls Z K U [ n − l ] ( s, z ) Q ( s, dz ) ) ds (3.18)Substituting the control variables Z, Γ and U with assumed form in (3.13), (3.14) and (3.15),and reordering the summation, one can confirm that (3.18) becomes V [ n ] t = X nT n ! ∂ nx H (0) − Z Tt Z [ n ] s dX cs − Z Tt Γ [ n ] s dY s − Z Tt Z K U [ n ] ( s, z ) N ( ds, dz ) − n X m =0 m X k =0 Z Tt X m − ks Y ks ( m − k )! k ! I ( m ≤ n − σ s (cid:16) I (1 ≤ k ) kv [ n − m − k +1 ,k − ( s ) + v [ n ] m − k +1 ,k ( s ) (cid:17) + k X l =0 C ( k,l ) β ( s ) ∂ ly λ ( s, v [ n − l ] m − k +1 ,k − l ( s ) + I (1 ≤ k ) kκ ( s ) v [ n − m − k,k ( s ) − k X l =0 C ( k,l ) ∂ ly λ ( s, (cid:16) n − m X j =1 v [ n − l ] j + m − k,k − l ( s ) q ( s, j ) j ! (cid:17) ds . (3.19)Then, matching the coefficients in the drift term of (3.17) to those of (3.19) yields the systemof the linear ODEs (3.16). The terminal conditions should be clear from the expression (3.19).As long as the forward SDE (3.17) is well-defined when using the solution of the ODEs(3.16), it actually gives one possible solution for the n -th order BSDE (3.18). Due to the lin-earity of the ODEs, the uniqueness of the solution within the assumed form should be clear. (cid:4) Note that the above system of ODEs can be easily solved one-by-one by evaluating in thefollowing order: n : 0 −→ n max (3.20) m : n −→ k : 0 −→ m . (3.22) Suppose that we have obtained the good estimate of moments of γ m = E [ X mT ] for m = 1 , , · · · from the truncated approximation of the BSDE (3.9) with H ( x ) = x m . The n -th order13umulant χ n is given, in terms of these moments, by χ n = n ! X { k m } ( − r − ( r − n X m =1 k m ! (cid:16) γ m m ! (cid:17) k m (3.23)where the summation P { k m } is taken for all the n-uplets of non-negative integers { k , · · · , k n } satisfying the Diophantine equation k + 2 k + · · · + nk n = n . (3.24) r is defined by r := k + k + · · · + k n .Then, the Edgeworth expansion of the X T ’s density using up to the n -th order cumulantis given by p n ( x ) = φ ( x ; µ, Σ ) n − X s =1 X { k m } s +2 r H s +2 r (cid:16) x − µ Σ (cid:17) s Y m =1 k m ! (cid:16) χ m +2 ( m + 2)! (cid:17) k m (3.25)where µ := χ , Σ := √ χ and φ ( x ; µ, Σ ) = 1 √ π Σ exp (cid:18) − (cid:16) x − µ Σ (cid:17) (cid:19) . (3.26)Here, the summation P { k m } is taken for all the s-uplets of non-negative integers satisfying k + 2 k + · · · + sk s = s (3.27)and r := k + k + · · · + k s (3.28)in (3.25). H n () denotes the Hermite polynomial defined by H n ( x ) := ( − n e x d n dx n e − x . (3.29)See, for example, Blinnikov and Moessner (1998) [4] for a simple derivation of the formulasand informative numerical examples of the density approximation from the moments.Then an approximated price of a Call option on S T with strike K based on the n -th order( n ≥
2) Edgeworth expansion is given by C Kn = Z ∞−∞ ( S e x − K ) + p n ( x ) dx = Z ∞ d ( S e Σ y + µ − K ) φ ( y ) n − X s =1 X { k m } s +2 r H s +2 r ( y ) s Y m =1 k m ! (cid:16) χ m +2 ( m + 2)! (cid:17) k m dy (3.30) We mean that the expansion using the cumulants ( χ i ) , i = 1 , , · · · , n . d := ln( K/S ) − µ Σ , φ ( y ) = 1 √ π e − y . (3.31)All the necessary integrations in (3.30) can be performed analytically thanks to the followingproperties of the Hermite polynomials: Z ∞ d φ ( y ) H n ( y ) dy = φ ( d ) H n − ( d ) (3.32) Z ∞ d e Σ y φ ( y ) H n ( y ) dy = e Σ d φ ( d ) H n − ( d ) + Σ Z ∞ d e Σ y φ ( y ) H n − ( y ) dy. (3.33)Put options can be evaluated similarly. For numerical examples, we choose a set of constant parameters and a Gaussian jump densitygiven as Q ( t, dz ) = 1 √ πσ J exp (cid:18) − (cid:16) z − µ J σ J (cid:17) (cid:19) dz . (3.34)In each of Figure 1 to 3, the approximation of the moments γ m = E [ X mT ] ( m = 1 , · · · , n = 20 based on the result of Lemma 1 is given in the left-hand panel. Each line is connecting one of the { γ m } estimated by the polynomial expansionup to the order specified by the horizontal axis. Note that the approximation of γ m becomesnon-zero only for n ≥ m . One can see that the lower-order moments converge rather quickly.The right-hand panel gives the comparison of the implied volatilities approximated by theEdgeworth expansion using the corresponding order of cumulants and the result from theMonte-Carlo simulation with 500,000 paths. We have used Put options for lower strikes bydirectly applying the corresponding formula without relying on the Put-Call parity. Thehorizontal axis denotes the size of the strikes scaled by S , i.e. K/S . In Figure 4, the highermoments ( γ , γ , γ ) are shown separately and the results for the implied volatilities aregiven in Figure 5 .As one can see from Figures 3 and 4, higher moments grow rapidly for longer maturitiesand also the rate of convergence slows down. As for higher moments there is no guaranteethat the Edgeworth expansion converges even if the moments are accurately estimated . Inaddition, by the very nature of polynomial expansion, when | γ m | ≫ It does not have a compact support but the scheme still seems to work well in this example. The estimation based on χ is omitted since it seems to give a totally useless result. Although it is similar, Gram-Charlier series generally gives much worse approximation.
Estimation of moments and implied volatilities. T = 1 , σ = 0 . , α = 0 . , ρ = − . , κ = 0 . , µ J = − . , σ J = 0 .
03 and λ ( t, Y t ) = 8( Y t + 1) . Figure 2:
Estimation of moments and implied volatilities. T = 1 , σ = 0 . , α = 0 . , ρ = 0 , κ = 0 . , µ J = − . , σ J = 0 .
03 and λ ( t, Y t ) = 8( Y t + 1) . Figure 3:
Estimation of moments and implied volatilities. T = 3 , σ = 0 . , α = 0 . , ρ = − . , κ = 0 . , µ J =0 . , σ J = 0 .
035 and λ ( t, Y t ) = 5 Y t + 10 Y t + 8. Estimation of moments. T = 5 , σ = 0 . , α = 0 . , ρ = − . , κ = 0 . , µ J = 0 . , σ J = 0 .
035 and λ ( t, Y t ) = 5 Y t + 10 Y t + 8. Figure 5:
Estimation of implied volatilities. T = 5 , σ = 0 . , α = 0 . , ρ = − . , κ = 0 . , µ J = 0 . , σ J =0 .
035 and λ ( t, Y t ) = 5 Y t + 10 Y t + 8. Before closing this section, let us study the path-wise nature of the current approximationscheme for the terminal condition H ( X T ) = X mT . For each order of moment m and expansion n , one can calculate the path-wise truncated approximation error [ X mT − e V ( n ) T ], where e V ( n ) is given by (2.13) appropriately specified for the current model. In Figure 6, we have shownthe scattered plot of this quantity for ( m = 1 and 5) with various orders of expansion n using the same setup as in Figure 3, i.e. { T = 3 , σ = 0 . , α = 0 . , ρ = − . , κ = 0 . , µ J =0 . , σ J = 0 .
035 and λ ( t, Y t ) = 5 Y t + 10 Y t + 8 } . In Table 1, the mean and standard deviationof [ X mT − e V ( n ) T ] are given for m = { , , · · · , } in the same setup. For ease of comparison, E [ X mT ] estimated by simulation is also given in the lower table for each moment. Note thatthe non-trivial approximation exists only for n ≥ m . Improvement of approximation stopseffectively at just a few higher order expansion n ≥ m , which means that the contributions ofpolynomial expansion for the target of X mT is dominated by m -th and just a couple of higherorder polynomials. This is rather natural and also consistent with the left panel of Figure 3showing the convergence of approximation series for each moment.One can observe that our scheme can provide accurate path-wise approximation of X m X T becomes more important for higher moments. For the above example, thesituation does not change meaningfully even if we use the pure diffusion model by putting λ = 0. We have observed a minor improvement of convergence only by a factor of few.Since we have used the standard Euler scheme, the corresponding simulation error may becontributing to the above result to some extent.Figure 6: Scattered plot of [ X mT − e V ( n ) T ] ( m = 1 in the left and m = 5 in the right) with various expansionorders n . The setup is the same as in Figure 3. m = 1 n = 0 n = 1 n = 2 n = 3 n = 5 n = 7 n = 10mean -0.054 − . × − − . × − . × − . × − . × − . × − stdev 0.34 0.028 4 . × − . × − . × − . × − . × − m = 2 n = 0 n = 2 n = 3 n = 4 n = 5 n = 7 n = 10mean 0.12 0.014 7 . × − − . × − . × − − . × − − . × − stdev 0.22 0.044 0.017 4 . × − . × − . × − . × − m = 3 n = 0 n = 3 n = 4 n = 5 n = 6 n = 7 n = 10mean -0.041 -0.017 − . × − . × − − . × − . × − . × − stdev 0.27 0.084 0.043 5 . × − . × − . × − . × − m = 4 n = 0 n = 4 n = 5 n = 6 n = 7 n = 8 n = 10mean 0.060 0.025 0.014 8 . × − . × − − . × − − . × − stdev 0.39 0.18 0.11 0.019 0.012 0.011 0.010 m = 5 n = 0 n = 5 n = 6 n = 7 n = 8 n = 9 n = 10mean -0.050 -0.037 -0.021 − . × − − . × − . × − − . × − stdev 0.65 0.40 0.28 0.084 0.042 0.025 0.024 m = 1 m = 2 m = 3 m = 4 m = 5 E [ X mT ] − . × − . × − − . × − . × − − . × − Table 1:
Mean and standard deviation of the path-wise realizations of [ X mT − e V ( n ) T ] of m = { , · · · , } withthe same setup as in Figure 3 by simulation. The second table gives the figures of E [ X mT ] with m = { , · · · , } estimated by MC simulation for clarity. An Application to λ -SABR Model By using an appropriate change of variables, the proposed polynomial expansion scheme canbe applied to a wider choice of models than what may be naively imagined. Let us consider(rescaled) λ -SABR model (SABR model with mean-reverting volatility) under an equivalentmartingale measure Q . S t = S + Z t ( S − β ) σ ( s ) ¯ Y s S βs dW s (4.1)¯ Y t = 1 + Z t (cid:16) α ( s ) ¯ Y s dB s + κ ( s )(1 − ¯ Y s ) ds (cid:17) (4.2)where W , B are one-dimensional Q -Brownian motions with d h W, B i t = ρ ( t ) dt . ( σ, ρ, κ ) areall deterministic functions, and β ∈ [0 ,
1) is a constant. Here, a factor of S − β is included tomake σ roughly equal to the at-the-money implied volatility of S .The change of variables X t := 11 − β (cid:18)(cid:16) S t S (cid:17) − β − (cid:19) (4.3) Y t := ¯ Y t − X t = Z t (cid:18) σ ( s )(1 + Y s ) dW s − β σ ( s ) b ( X s )(1 + Y s ) ds (cid:19) (4.5) Y t = Z t (cid:16) α ( s )(1 + Y s ) dB s − κ ( s ) Y s ds (cid:17) (4.6)where b ( x ) := 11 + (1 − β ) x . (4.7)The assumption on the quadratic covariation (2.4) is now satisfied for these new variables.The BSDE relevant for a European contingent claim with terminal payoff H ( X T ) atmaturity T is given by V t = H ( X T ) − Z Tt (cid:18) β σ ( s ) b ( X s )(1 + Y s ) Z s + κ ( s ) Y s Γ s (cid:19) ds − Z Tt Z s dX s − Z Tt Γ s dY s (4.8)As in the Heston’s case, we choose H ( x ) = x m , ( m = 1 , , · · · ) to obtain the moment estimateof X T and then use the Edgeworth expansion to approximate its probability density. Here,we are not claiming the Edgeworth expansion is the best choice and different basis functions(such as Laguerre polynomials) can be more appropriate.19 .2 Polynomial Expansion We now introduce ǫ to the BSDE (4.8) so that we can perform polynomial expansion V ǫt = H ( ǫX T ) − Z Tt (cid:18) β σ ( s ) b ( ǫX s )(1 + ǫY s ) Z ǫs + ǫκ ( s ) Y s Γ ǫs (cid:19) ds − Z Tt Z ǫs dX s − Z Tt Γ ǫs dY s (4.9)as V ǫt = ∞ X n =0 ǫ n V [ n ] t (4.10) Z ǫt = ∞ X n =0 ǫ n Z [ n ] t , Γ ǫt = ∞ X n =0 ǫ n Γ [ n ] t . (4.11)We have the following lemma. Lemma 2
If it exists, the polynomial solution for the expansion in (4.10) and (4.11) isuniquely given by V [ n ] t = n X m =0 m X k =0 X m − kt Y kt ( m − k )! k ! v [ n ] m − k,k ( t ) (4.12) Z [ n ] t = n X m =1 m − X k =0 X m − k − t Y kt ( m − k − k ! v [ n ] m − k,k ( t ) (4.13)Γ [ n ] t = n X m =1 m X k =1 X m − kt Y k − t ( m − k )!( k − v [ n ] m − k,k ( t ) (4.14) with the set of deterministic functions v [ n ] m − k,k ( t ) of (0 ≤ k ≤ m ≤ n ) satisfying the followingrecursive system of linear ODEs ˙ v [ n ] m − k,k ( t ) = − I (2 ≤ k ) k ( k − (cid:18) σ t v [ n ] m − k +2 ,k − ( t ) + ρ t σ t α t v [ n ] m − k +1 ,k − ( t ) + α t v [ n ] m − k,k ( t ) (cid:19) − I ( m ≤ n − , ≤ k ) k (cid:16) σ t v [ n ] m − k +2 ,k − ( t ) + 2 ρ t σ t α t v [ n ] m − k +1 ,k ( t ) + α t v [ n ] m − k,k +1 ( t ) (cid:17) − I ( m ≤ n − (cid:18) σ t v [ n ] m − k +2 ,k ( t ) + ρ t σ t α t v [ n ] m − k +1 ,k +1 ( t ) + α t v [ n ] m − k,k +2 ( t ) (cid:19) + I ( m ≤ n − , ≤ k ) kκ t v [ n − m − k,k ( t ) + I ( m ≤ n − β σ t m − k X l =0 C ( m − k, l ) ∂ lx b (0) × v [ n − l ] m − k − l +1 ,k ( t ) + I (1 ≤ k ) k v [ n − l − m − k − l +1 ,k − ( t ) + I (2 ≤ k ) k ( k − v [ n − l − m − k − l +1 ,k − ( t ) ! (4.15)20 aving the terminal conditions v [ n ] n, ( T ) = ∂ nx H (0) with all the other components zero.Proof: It can be proved in exactly the same way as Lemma 1. The derivation is givenin the Appendix A.
As in Section 3.4, let us provide several numerical examples for the estimated moments andthe comparison of the implied volatilities. The number of paths for Monte-Carlo simulationis 500 ,
000 as before. For this model, we cannot use the special relation in (3.32) and (3.33),and hence we have carried out numerical integration of the estimated density for the pricing.The styles and conventions used in each figures are the same as those in Section 3.4.Although the polynomial expansion gives similar accuracy for short maturities, its appli-cability to long maturities is rather limited compared to the previous extended Heston model.The main cause seems to be the factor k ( k −
1) appearing in the first line of the ODE givenin Lemma 2, which strongly drives { v [ n ] m,k } especially for higher moments and makes themunable to converge. This factor stems from the terms ∝ Y in the quadratic covariations. Inaddition, since the support of X t is limited to the range X t ≥ − − β , the model’s compat-ibility to the Edgeworth expansion may be lower than the Heston model. This may be oneof the reasons for somewhat unstable behavior of the implied volatilities when higher-ordercumulants are included.For completeness, we give a convergence analysis for the path-wise realizations of thetruncated approximation [ X mT − e V ( n ) T ] as before. In Table 2, the mean and standard deviationfor m = { , , · · · , } with various order of expansions are given under the same setup usedin Figure 8. In this model, the improvement of approximation stops more quickly than theprevious Heston model case. This is likely due to the smaller size of moments and possiblyother delicate model features. The quicker convergence of approximation series can also beseen from the left panel of Figure 8.Figure 7: Estimation of moments and implied volatilities. T = 0 . , σ = 0 . , α = 0 . , ρ = − . , κ =0 . , β = 0 . This is due to the performed change of parameters.
Estimation of moments and implied volatilities. T = 1 , σ = 0 . , α = 0 . , ρ = − . , κ = 0 . , β =0 . Figure 9:
Estimation of moments and implied volatilities. T = 1 , σ = 0 . , α = 0 . , ρ = 0 , κ = 0 . , β = 0 . m = 1 n = 0 n = 1 n = 2 n = 3 n = 5 n = 7 n = 10mean − . × − − . × − − . × − . × − − . × − . × − . × − stdev 0.15 1 . × − . × − . × − . × − . × − . × − m = 2 n = 0 n = 2 n = 3 n = 4 n = 5 n = 7 n = 10mean 0.024 2 . × − . × − . × − . × − . × − . × − stdev 0.037 1 . × − . × − . × − . × − . × − . × − m = 3 n = 0 n = 3 n = 4 n = 5 n = 6 n = 7 n = 10mean − . × − − . × − − . × − − . × − − . × − − . × − − . × − stdev 0.018 6 . × − . × − . × − . × − . × − . × − m = 4 n = 0 n = 4 n = 5 n = 6 n = 7 n = 8 n = 10mean 1 . × − . × − . × − . × − . × − . × − . × − stdev 9 . × − . × − . × − . × − . × − . × − . × − m = 5 n = 0 n = 5 n = 6 n = 7 n = 8 n = 9 n = 10mean − . × − − . × − − . × − − . × − − . × − − . × − − . × − stdev 5 . × − . × − . × − . × − . × − . × − . × − m = 1 m = 2 m = 3 m = 4 m = 5 E [ X mT ] − . × − . × − − . × − . × − − . × − Table 2:
Mean and standard deviation of the path-wise realizations of [ X mT − e V ( n ) T ] of m = { , · · · , } with thesetup in Figure 8 by simulation. The second table gives the figures of E [ X mT ] with m = { , · · · , } estimatedby MC simulation for clarity. Utility Optimization with Terminal Liability
European contingent claims, which we studied in the previous sections, can of course be solvedwithout resorting to a complicated BSDE formulation. The main motivation there was toget some insight about the performance of the proposed scheme by studying the two popularmodels. Now, in this section, we treat a utility-optimization problem in an incomplete marketwhere solving a BSDE becomes crucially important.Here, we adopt a simple Heston security market consists of one-risky asset with stochasticvolatility. For simplicity, we assume that the interest rate is zero. In the probability spaceof the physical measure (Ω , F , P ), the dynamics of the underlying variables is assumed to begiven by S t = S + Z t S s σ ( s ) p ¯ Y s (cid:16) dW s + ¯ θ ( s, S s , ¯ Y s ) ds (cid:17) ¯ Y t = 1 + Z t α ( s ) p ¯ Y s (cid:16) dB s + ρ ( s )¯ θ ( s, S s , ¯ Y s ) ds (cid:17) + κ ( s )(1 − ¯ Y s ) ds ! (5.1)where W, B are P -Brownian motions with d h W, B i t = ρ ( t ) dt . σ, α and κ are deterministicfunctions of time, and ¯ θ : [0 , T ] × R → R gives the risk-premium process associating with W . The risk-premium for B is implied by the ρ ¯ θ as well as the mean-reverting term of ¯ Y .Given a portfolio strategy ( π t ) t ≥ , the wealth at the terminal time T ( > t ) is given by W πT ( t, w ) = w + Z Tt π u dS u . (5.2)In the reminder of this section, we are going to study the BSDE associated with the expo-nential cost minimization: V ( t, w ) = ess inf π E " exp γ (cid:16) ¯ H ( S T , Y T ) − W πT ( t, w ) (cid:17)!(cid:12)(cid:12)(cid:12) F t (5.3)where γ is a positive constant specifying the risk averseness, and ¯ H : R → R is a smoothfunction denoting the terminal liability. Using Itˆo-Ventzell formula and the transformation V ( t ) = ln (cid:16) V ( t, w ) e γw (cid:17) (5.4)one can show that the following BSDE holds: V t = γ ¯ H ( S T , ¯ Y T ) − Z Tt (cid:26)
12 ¯ θ ( s, S s , ¯ Y s ) −
12 (1 − ρ ( s ) )¯Γ s (cid:27) ds − Z Tt ¯ Z s h dW s + ¯ θ ( s, S s , ¯ Y s ) ds i − Z Tt ¯Γ s h dB s + ρ ( s )¯ θ ( s, S s , ¯ Y s ) ds i . (5.5)It is well-known that the transformation (5.4) makes V ( t ) independent from the initial wealth w . The details and various interesting topics can be found in a comprehensive review by Mania& Tevzadze (2008) [28]. Similar qgBSDE arises in other economically important setups too,23uch as Power and HARA (hyperbolic absolute risk aversion) utilities after appropriate changeof variables. We simply study (5.5) for a demonstrative purpose of the current approximationscheme.It is important to notice that one cannot make use of Cole-Hopf transformation to convertthe qgBSDE (5.5) to a solvable linear BSDE as long as ( ¯ H , ¯ θ ) depend on both of the S and¯ Y . For example, if both of them depend only on ¯ Y , one can solve it analytically by followingthe arguments given by Zariphopoulou (2001) [37].As in Section 3, we perform the change of variables X t = ln (cid:18) S t S (cid:19) (5.6) Y t = ¯ Y t − θ ( s, X s , Y s ) := ¯ θ ( s, S s , ¯ Y s ) . (5.8)The relevant forward SDEs are now given by X t = Z t ( σ ( s ) p Y s + 1 (cid:16) dW s + θ ( s, X s , Y s ) ds (cid:17) − σ ( s ) Y s + 1) ds ) (5.9) Y t = Z t ( α ( s ) p Y s + 1 (cid:16) dB s + ρ ( s ) θ ( s, X s , Y s ) ds (cid:17) − κ ( s ) Y s ds ) . (5.10)Simple redefinition of variables yields V t = γH ( X T , Y T ) − Z Tt Z s dX s − Z Tt Γ s dY s − Z Tt (cid:26)
12 Θ( s, X s , Y s ) − α ( s ) − ρ ( s ) )(1 + Y s )Γ s + σ ( s ) Y s ) Z s + κ ( s ) Y s Γ s (cid:27) ds (5.11)where H ( X T , Y T ) := ¯ H ( S T , ¯ Y T ) and Θ( s, X s , Y s ) = θ ( s, X s , Y s ) . The control variables areconnected to those in (5.5) by¯ Z s = Z s σ ( s ) p Y s + 1 , ¯Γ s = Γ s α ( s ) p Y s + 1 . (5.12)We assume the system of the forward and backward SDEs (5.9), (5.10) and (5.11) has a well-posed solution in the reminder of the section. Although it deviates from the main subject ofthe paper, it is interesting to notice that the above BSDE has a simple exact solution in aspecial case. The details are give in Appendix C.24 .1 Polynomial Expansion In order to obtain the polynomial approximation for the system (5.9), (5.10) and (5.11), letus introduce ǫ and consider the perturbed BSDE: V ǫt = γH ( ǫX T , ǫY T ) − Z Tt Z ǫs dX s − Z Tt Γ ǫs dY s − Z Tt (
12 Θ( s, ǫX s , ǫY s ) − α ( s ) − ρ ( s ) )(1 + ǫY s )[Γ ǫs ] + σ ( s ) ǫY s ) Z ǫs + ǫκ ( s ) Y s Γ ǫs ) ds (5.13)and the associated expansion V ǫt = ∞ X n =0 ǫ n V [ n ] t (5.14) Z ǫt = ∞ X n =0 ǫ n Z [ n ] t , Γ ǫt = ∞ X n =0 ǫ n Γ [ n ] t . (5.15)The approximate solution of the original system is obtained by truncating the summation ata certain order n and putting ( ǫ = 1) as explained in Section 2.2. For this model, we havethe following result: Lemma 3
If it exists, the polynomial solution for the expansion in (5.14) and (5.15) isuniquely given by V [ n ] t = n X m =0 m X k =0 X m − kt Y kt ( m − k )! k ! v [ n ] m − k,k ( t ) (5.16) Z [ n ] t = n X m =1 m − X k =0 X m − k − t Y kt ( m − k − k ! v [ n ] m − k,k ( t ) (5.17)Γ [ n ] t = n X m =1 m X k =1 X m − kt Y k − t ( m − k )!( k − v [ n ] m − k,k ( t ) (5.18) with the set of deterministic functions v [ n ] m − k,k ( t ) of (0 ≤ k ≤ m ≤ n ) satisfying the following ecursive system of linear ODEs ˙ v [ n ] m − k,k ( t ) = − I ( m ≤ n − , ≤ k ) k (cid:18) σ t v [ n ] m − k +2 ,k − ( t ) + ρ t σ t α t v [ n ] m − k +1 ,k ( t ) + α t v [ n ] m − k,k +1 ( t ) (cid:19) − I ( m ≤ n − (cid:18) σ t v [ n ] m − k +2 ,k ( t ) + ρ t σ t α t v [ n ] m − k +1 ,k +1 ( t ) + α t v [ n ] m − k,k +2 ( t ) (cid:19) + I ( m = n ) ∂ n − kx ∂ ky Θ( t, ,
0) + I ( m ≤ n − σ t v [ n ] m − k +1 ,k ( t )+ I ( m ≤ n − , ≤ k ) k (cid:18) σ t v [ n − m − k +1 ,k − ( t ) + κ t v [ n − m − k,k ( t ) (cid:19) − I ( m ≤ n − n − X l =1 l ∧ [ m +1] X j =1 ∨ [ l +2 − n + m ] j ∧ [ k +1] X p =1 ∨ [ j − m + k ] α t ξ t C ( m − k,j − p ) C ( k,p − v [ l ] j − p,p ( t ) v [ n − l ] m − k − j + p,k − p +2 ( t ) − I (1 ≤ m ≤ n − , ≤ k ) n − X l =1 l ∧ m X j =1 ∨ [ l +2 − n + m ] j ∧ k X p =1 ∨ [ j − m + k ] α t ξ t C ( m − k,j − p ) C ( k,p ) pv [ l ] j − p,p ( t ) v [ n − l − m − k − j + p,k − p +1 ( t )(5.19) with ξ t := (1 − ρ ( t ) ) and the terminal conditions v [ n ] n − k,k ( T ) = γ∂ n − kx ∂ ky H (0 , with all theother components zero.Proof: The proof is done in a similar way to Lemma 1 and 2. The details of the deriva-tion are given in Appendix B.
For numerical examples, we shall useΘ( t, X t , Y t ) := c e − c X t ( Y t + 1) (5.20) H ( X T , Y T ) := e − g X T G ( Y T ) (5.21)where c , c , g are constants and G ( · ) a smooth function of Y . Since the parameter ofrisk-averseness γ appears only in a combination γH , the factor e − g X T can equivalently beinterpreted as a S T -dependent risk averseness.The problem analyzed in this section is intrinsically non-linear. Thus, we cannot usethe density approximation and must directly approximate the terminal payoff by a smoothfunction. In practice, however, it should not be a prohibitive limitation. Since the problemis non-linear, one has to consider the optimization in a portfolio level. Then, consideringan appropriate hedging strategy based on a smooth approximate payoff function, instead oftreating it exactly, should be reasonable.We consider the next four choices of terminal liability (except e − g x factor) in the numer-26cal examples: (1) : sin (cid:16) y + π (cid:17) (5.22)(2) : max (cid:0) , y (cid:1) (5.23)(3) : max (cid:0) , − y (cid:1) (5.24)(4) : 0 . − max (cid:0) , . − y (cid:1) (5.25)For (2) to (4), we have approximated it by a 5th-order polynomial function determined bya simple least-square method, and treat it as the true G ( y ) in the evaluation. Here, theshapes of the liability and the order of approximating polynomial function are chosen ratherarbitrary. In practice, one has to consider in a portfolio level and needs to choose a certainorder of polynomial to recover its overall shape. The impact from adding another term wouldbe easy to check directly. It is naturally expected, however, that the higher order terms playsonly a minor role otherwise it means that the firm is taking quite problematic positions andexposing it to the far-tail behavior of the underlying securities.Each of Figure 10 to 13 consists of: 1)Top left: a graph of G ( y ), 2)Top right: a graphof the truncated value function and control variables ( V ( n )0 , Z ( n )0 , Γ ( n )0 ) for each n specifiedby the horizontal axis , 3)Bottom left: a scattered plot of [ γH ( X T , Y T ) − e V ( n ) T ] for eachexpansion order, 4) Bottom right: a graph of the means as well as the standard deviations of[ γH ( X T , Y T ) − e V ( n ) T ] for (0 ≤ n ≤
10) with 100,000-path simulation, whose details are alsogiven in a table associated with each example. Note that the errors are measured relative tothe smoothly modified terminal functions in (2) to (4) cases.From the definition of the truncated approximation, one can easily see that the mean of[ γH ( X T , Y T ) − e V ( n ) T ] is equivalent to the estimate of V ( n )0 − E (cid:20) γH ( X T , Y T ) − Z T Z ( n ) t dX t − Z T Γ ( n ) t dY t − Z T n
12 Θ( t, X t , Y t ) − α − ρ )(1 + Y t )[Γ ( n ) t ] + σ Y t ) Z ( n ) t + κY t Γ ( n ) t o dt (cid:21) (5.26)by simulation. Its convergence to zero gives one consistency test for the value function atthe initial point. The scattered plots of [ γH ( X T , Y T ) − e V ( n ) T ] and the corresponding standarddeviations provide a much stronger test. They suggest that the truncated value functionsand control variables give a good path-wise approximation for the original BSDE. One canclearly observe that the deviations [ γH ( X T , Y T ) − e V ( n ) T ] at the maturity are strongly clusteringaround zero even for a relatively low expansion order n ∼
3. Of course, as one can imagine,the probability that the size of (cid:0) γH ( X T , Y T ) , X T , Y T (cid:1) becomes (meaningfully) bigger thanone should be small enough in order to obtain a converging result. This means that we needto adopt a proper “scaling” for the wealth and the other parameters to make sure the chosenutility (or cost function) stays O (1) . See, (2.11) and (2.12) for the definition of truncated variables. A similar scaling would be necessary for any risk-management in practice in the presence of non-linearities. emark For those who have checked the result of Lemma 3 by themselves, it must be clear thatderiving a closed form system of ODEs would be much harder in a realistic multi-asset setup.In that sense, the above numerical result is quite encouraging by implying that one may get areasonable approximation even by a lower order expansion, for example, n ∼
4. In this case,step-by-step derivation of the relevant ODEs following the instruction given in Section 2.2can be done without much difficulty even for a more involved BSDE. Interesting practicalapplications are left for the future research.Figure 10: T = 1 , σ = 0 . , α = 0 . , ρ = − . , κ = 0 . , c = 0 . c = 0 . , γ = 1 , g = 0 . G ( y ) = sin( y + π/ n = 0 n = 1 n = 2 n = 3 n = 4 n = 5 n = 7 n = 10mean -0.026 0.016 -0.016 -0.011 6 . × − . × − − . × − . × − stdev 0.39 0.13 0.094 0.027 0.030 0.020 0.012 9 . × − Table 3: Mean and standard deviation of [ γH T − e V ( n ) T ] for the setup in Figure 10.28igure 11: T = 1 , σ = 0 . , α = 0 . , ρ = − . , κ = 0 . , c = 0 . c = 0 . , γ = 1 , g = 0 . G ( y ) is a 5-th order polynomial approximating max(0 , y ). n = 0 n = 1 n = 2 n = 3 n = 4 n = 5 n = 7 n = 10mean 0.093 0.10 − . × − . × − − . × − − . × − . × − . × − stdev 0.29 0.14 0.040 0.030 0.026 0.011 0.011 5 . × − Table 4: Mean and standard deviation of [ γH T − e V ( n ) T ] for the setup in Figure 11.29igure 12: T = 1 , σ = 0 . , α = 0 . , ρ = − . , κ = 0 . , c = 0 . c = 0 . , γ = 1 , g = 0 . G ( y ) is a 5-th order polynomial approximating max(0 , − y ). n = 0 n = 1 n = 2 n = 3 n = 4 n = 5 n = 7 n = 10mean 0.064 0.075 -0.019 − . × − . × − − . × − . × − . × − stdev 0.17 0.089 0.044 0.042 0.041 0.011 0.012 8 . × − Table 5: Mean and standard deviation of [ γH T − e V ( n ) T ] for the setup in Figure 12.30igure 13: T = 1 , σ = 0 . , α = 0 . , ρ = − . , κ = 0 . , c = 0 . c = 0 . , γ = 1 , g = 0 . G ( y ) is a 5-th order polynomial approximating [0 . − max(0 , . − y )]. n = 0 n = 1 n = 2 n = 3 n = 4 n = 5 n = 7 n = 10mean -0.052 -0.031 0.012 0.011 − . × − − . × − . × − − . × − stdev 0.27 0.095 0.040 0.069 0.024 0.013 0.017 5 . × − Table 6: Mean and standard deviation of [ γH T − e V ( n ) T ] for the setup in Figure 13.31 Conclusions
In this paper, a polynomial scheme of asymptotic expansion for BSDEs is proposed. Wehave shown that the polynomial expansion is uniquely determined by the recursive systemof linear ODEs, which can be easily solved one-by-one by following the appropriate orderof evaluation. We have studied possible applications to the pricing of European contingentclaims as well as the exponential-utility optimization with terminal liability, each of which isprovided several illustrative numerical examples.A rigorous mathematical justification and more intensive numerical studies with realisticmodels are left for the future works. For example, a class of multi-factor Heston model pro-posed by Col et al. (2013) [8] has a nice structure of dependence to which the current schemecan be applied. Studying the BSDEs associated with the control problem with defaultablesecurities, such as those give by Pham (2010) [30], looks interesting, too.
A Proof of Lemma 2
We proceed as in the Heston’s model. Based on the dynamics (4.5) and (4.6), the forwarddynamics of the hypothesized polynomial solution is given by dV [ n ] t = n X m =0 m X k =0 X m − kt Y kt ( m − k )! k ! ( ˙ v [ n ] m − k,k ( t )+ I (2 ≤ k ) k ( k − (cid:18) σ t v [ n ] m − k +2 ,k − ( t ) + ρ t σ t α t v [ n ] m − k +1 ,k − ( t ) + α t v [ n ] m − k,k ( t ) (cid:19) + I ( m ≤ n − , ≤ k ) k (cid:16) σ t v [ n ] m − k +2 ,k − ( t ) + 2 ρ t σ t α t v [ n ] m − k +1 ,k ( t ) + α t v [ n ] m − k,k +1 ( t ) (cid:17) + I ( m ≤ n − (cid:18) σ t v [ n ] m − k +2 ,k ( t ) + ρ t σ t α t v [ n ] m − k +1 ,k +1 ( t ) + α t v [ n ] m − k,k +2 ( t ) (cid:19) ) dt + n X m =1 m − X k =0 v [ n ] m − k,k ( t ) X m − k − t Y kt ( m − k − k ! dX t + n X m =1 m X k =1 v [ n ] m − k,k ( t ) X m − kt Y k − t ( m − k )!( k − dY t (A.1)which implies the control variables given in (4.13) and (4.14).On the other hand, the n -th order part of the BSDE (4.9) is V [ n ] t = X nT n ! ∂ nx H (0) − Z Tt Z [ n ] s dX s − Z Tt Γ [ n ] s dY s − Z Tt ( β σ ( s ) n − X l =0 ∂ lx b (0) l ! X ls Z [ n − l ] s + 2 n − X l =0 ∂ lx b (0) l ! X ls Y s Z [ n − l − s + n − X l =0 ∂ lx b (0) l ! X ls Y s Z [ n − l − s ! + κ s Y s Γ [ n − s ) ds (A.2)32ubstituting the control variables by those in (4.13) and (4.14), one obtains V [ n ] t = X nT n ! ∂ nx H (0) − Z Tt Z [ n ] s dX s − Z Tt Γ [ n ] s dY s − n X m =0 m X k =0 Z Tt X m − ks Y ks ( m − k )! k ! × ( I ( m ≤ n − , ≤ k ) kκ s v [ n − m − k,k ( s ) + I ( m ≤ n − β σ s m − k X l =0 C ( m − k,l ) ∂ lx b (0) × (cid:16) v [ n − l ] m − k − l +1 ,k ( s ) + I (1 ≤ k ) k v [ n − l − m − k − l +1 ,k − ( s ) + I (2 ≤ k ) k ( k − v [ n − l − m − k − l +1 ,k − ( s ) (cid:17) ) ds . (A.3)By comparing the coefficients in the drift term, one obtains the linear ODEs as (4.15). If theforward SDE (A.1) is well-defined with the solution of (4.15), then it is clear that it gives onesolution for the BSDE (A.2). The uniqueness of the polynomial solution is clear due to thelinearity of the ODEs. B Proof of Lemma 3
The dynamics of (5.9) and (5.10) gives the forward SDEs of the assumed polynomial (5.16)as dV [ n ] t = n X m =0 m X k =0 X m − kt Y kt ( m − k )! k ! ( ˙ v [ n ] m − k,k ( t )+ I ( m ≤ n − , ≤ k ) k (cid:18) σ t v [ n ] m − k +2 ,k − ( t ) + ρ t σ t α t v [ n ] m − k +1 ,k ( t ) + α t v [ n ] m − k,k +1 ( t ) (cid:19) + I ( m ≤ n − (cid:18) σ t v [ n ] m − k +2 ,k ( t ) + ρ t σ t α t v [ n ] m − k +1 ,k +1 ( t ) + α t v [ n ] m − k,k +2 ( t ) (cid:19) ) dt + n X m =1 m − X k =0 v [ n ] m − k,k ( t ) X m − k − t Y kt ( m − k − k ! dX t + n X m =1 m X k =1 v [ n ] m − k,k ( t ) X m − kt Y k − t ( m − k )!( k − dY t (B.1)33hich then implies the control variables as in (5.17) and (5.18). On the other hand, the n -order BSDE of (5.13) is V [ n ] t = n X k =0 X n − kT Y kT ( n − k )! k ! γ∂ n − kx ∂ ky H (0 , − Z Tt Z [ n ] s dX s − Z Tt Γ [ n ] s dY s − Z Tt ( n X k =0 X n − ks Y ks ( n − k )! k ! 12 ∂ n − kx ∂ ky Θ( s, ,
0) + κ s Y s Γ [ n − s − α ( s ) − ρ ( s ) ) n − X l =1 Γ [ l ] s Γ [ n − l ] s + Y s n − X l =1 Γ [ l ] s Γ [ n − l − s ! + σ ( s ) (cid:16) Z [ n ] s + Y s Z [ n − s (cid:17) ) ds . (B.2)Substituting (5.17) and (5.18) into the above expression and reordering the summation yield V [ n ] t = n X k =0 X n − kT Y kT ( n − k )! k ! γ∂ n − kx ∂ ky H (0 , − Z Tt Z [ n ] s dX s − Z Tt Γ [ n ] s dY s − n X m =0 m X k =0 Z Tt X m − ks Y ks ( m − k )! k ! ( I ( m = n ) ∂ n − kx ∂ ky Θ( s, ,
0) + I ( m ≤ n − σ s v [ n ] m − k +1 ,k ( s )+ I ( m ≤ n − , ≤ k ) k (cid:16) σ s v [ n − m − k +1 ,k − ( s ) + κ s v [ n − m − k,k ( s ) (cid:17) − α s ξ s I ( m ≤ n − n − X l =1 l ∧ [ m +1] X j =1 ∨ [ l +2 − n + m ] j ∧ [ k +1] X p =1 ∨ [ j − m + k ] C ( m − k,j − p ) C ( k,p − v [ l ] j − p,p ( s ) v [ n − l ] m − k − j + p,k − p +2 ( s ) − α s ξ s I (1 ≤ m ≤ n − , ≤ k ) n − X l =1 l ∧ m X j =1 ∨ [ l +2 − n + m ] j ∧ k X p =1 ∨ [ j − m + k ] C ( m − k,j − p ) C ( k,p ) pv [ l ] j − p,p ( s ) v [ n − l − m − k − j + p,k − p +1 ( s ) ds (B.3)By comparing the drift terms of (B.1) and (B.3), one obtains the system of linear ODEsas in Lemma 3. It is clear that if the forward SDE (B.1) with the solution of the ODEs iswell-defined, it at least gives one solution for the BSDE of the n -th order (B.2). Due to thelinearity of the ODEs, the solution should be unique within the assumed polynomial form. C An exact solution for (5.11)
Suppose that both of the H ( x, y ) and Θ( t, x, y ) are linear functions of ( x, y ): H ( x, y ) = h x x + h y y + h (C.1)Θ( t, x, y ) = Θ x ( t ) x + Θ y ( t ) y + Θ ( t ) (C.2)34here ( h x , h y , h ) are constants and (Θ x ( t ) , Θ y ( t ) , Θ ( t )) t are some deterministic functions oftime. Then, it is almost immediate to notice that a linear value function and the deterministiccontrol variables can provide an exact solution: V t = v x ( t ) X t + v y ( t ) Y t + v ( t ) (C.3) Z t = v x ( t ) (C.4)Γ t = v y ( t ) (C.5)where ( v x ( t ) , v y ( t ) , v ( t )) are deterministic functions of time.The ODEs which fixes v x , v y , v can be obtained quite similarly to the discussed approx-imation scheme. On the one hand, the dynamics of the proposed solution becomes dV t = (cid:16) ˙ v x ( t ) X t + ˙ v y ( t ) Y t + ˙ v ( t ) (cid:17) dt + v x ( t ) dX t + v y ( t ) dY t . (C.6)On the other hand, by inserting the assumed form of control variables to (5.11), one obtains V t = γ h h x X T + h y Y T + h i − Z Tt v x ( s ) dX s − Z Tt v y ( s ) dY s − Z Tt (cid:26) (cid:16) Θ x ( s ) X s + Θ y ( s ) Y s + Θ ( s ) (cid:17) + κ ( s ) v y ( s ) Y s (C.7)+ 12 (cid:16) σ ( s ) v x ( s ) − α ( s ) (1 − ρ ( s ) ) v y ( s ) (cid:17) (1 + Y s ) (cid:27) ds . (C.8)Thus, the the system of ODEs including a Riccati-type for v y ˙ v x ( t ) = 12 Θ x ( t ) (C.9)˙ v y ( t ) = 12 (cid:16) σ ( t ) v x ( t ) − α ( t ) (1 − ρ ( s ) ) v y ( t ) (cid:17) + κ ( t ) v y ( t ) + 12 Θ y ( t ) (C.10)˙ v ( t ) = 12 (cid:16) σ ( t ) v x ( t ) − α ( t ) (1 − ρ ( s ) ) v y ( t ) (cid:17) + 12 Θ ( t ) (C.11)with the terminal conditions v x ( T ) = γh x , v y ( T ) = γh y , v ( T ) = γh (C.12)gives an exact solution if v y (and hence the others) has a finite solution for the relevant timeinterval t ∈ [0 , T ]. 35 cknowledgement This research is partially supported by Center for Advanced Research in Finance (CARF).The author is grateful to professor Takahashi for many helpful discussions and encourage-ments. The author also thanks anonymous referees whose comments significantly clarifies thepresentation of the material.