Computing the CEV option pricing formula using the semiclassical approximation of path integral
CComputing the CEV option pricing formula usingthe semiclassical approximation of path integral
Axel A. Araneda ∗ and Marcelo J. Villena †‡ Last version: March 29, 2018
Abstract
The Constant Elasticity of Variance (CEV) model significantly outperforms the Black-Scholes (BS)model in forecasting both prices and options. Furthermore, the CEV model has a marked advantage incapturing basic empirical regularities such as: heteroscedasticity, the leverage effect, and the volatilitysmile. In fact, the performance of the CEV model is comparable to most stochastic volatility models,but it is considerable easier to implement and calibrate. Nevertheless, the standard CEV model solution,using the non-central chi-square approach, still presents high computational times, specially when: i) thematurity is small, ii) the volatility is low, or iii) the elasticity of the variance tends to zero. In this paper,a new numerical method for computing the CEV model is developed. This new approach is based onthe semiclassical approximation of Feynman’s path integral. Our simulations show that the method isefficient and accurate compared to the standard CEV solution considering the pricing of European calloptions.
Keywords: Option pricing, constant elasticity of variance model, path integral, numerical methods. ∗ Faculty of Engineering & Sciences, Universidad Adolfo Ibáñez, Avda. Diagonal las Torres 2640, Peñalolén, 7941169 Santiago,Chile. Email: [email protected]. † Universidad Adolfo Ibáñez, Address: Diagonal Las Torres 2640, Peñalolén, Santiago, Chile. Telephone: 56-223311491,Email: [email protected]. ‡ Aid from the Fondecyt Program, project Nº 1131096, is grateful acknowledged by Marcelo. a r X i v : . [ q -f i n . C P ] M a r Introduction
One of the most significant limitations of the Black-Scholes (BS) [1] model is the assumption of constantvolatility, which ignores some well-known empirical regularities such as: the leverage effect [2, 3], and thevolatility smile [4, 5]. These shortcomings have inspired several non-constant volatility models in continuoustime , considering ‘stochastic volatility’ or ‘level-dependent volatility’ models [10]. In the former, both theasset and the volatility have their own diffusion processes. In the level-dependent volatility models only theasset is governed by a diffusion process, and its volatility is modeled in function of the asset level. In thispaper, the analysis will be focused on the the constant elasticity of variance (CEV) model, proposed by J.Cox, the most known level-dependent volatility approach [11, 12].Furthermore, the CEV model has a marked advantage in capturing basic empirical regularities such as:heteroscedasticity, the leverage effect, and the volatility smile[13–16]. As a consequence, the CEV modelsignificantly outperforms the Black-Scholes (BS) model in forecasting both prices and options [17–21]. Fur-thermore, the performance of the CEV model is comparable to most stochastic volatility models, but it isconsiderable easier to implement and calibrate [22].In terms of the option pricing using the CEV model, the exact formula for a vanilla European optioninvolves a complex computation of an infinite series of incomplete gamma functions [11]. Subsequently, [14]matched the Cox pricing formula with the non-central chi-square distribution. Schroder also provides asimple approximated method for its computation, see [23] for a detailed derivation of the two methodologies.Since then, the use of the non-central chi-square distribution becomes the most widely used method of pricingfor options under the CEV model. Besides, several alternative methods for its implementation have beendeveloped [24].Nevertheless, the standard CEV model solutions, using the non-central chi-square approach, still presentsconsiderably high computational times, specially when: i) the maturity is small, ii) the volatility is low, or iii)the elasticity of the variance tends to zero[14, 25, 26]. In order to deal with these problems, many approacheshave been reported for the European-vanilla type option pricing. These approaches include numerical schemes[25, 27–29], Montecarlo simulations [30], perturbation theory model [31], and analytical approximations tothe transition density [32] or to the hedging strategy [33], among others.In this paper, a new numerical method for computing the CEV model is developed. This new approachis based on the semiclassical approximation of Feynman’s path integral model. In financial literature, pathintegral techniques have already been used in the option pricing problem, see [34–39]. Nevertheless, the mainfocus has been in theoretical issues rather than in practical applications. On the other hand, the application ofthe semiclassical approximation of Feynman’s path integral technique on financial problems is rather limited,see for example [40–42].[42] points out that, just as in the case of quantum mechanics, the path integral approach in finance isneither a panacea, nor is it intended to yield fundamentally new results, but in some cases it provides clarityand insight into old problems. In this paper, we analyze the possibility that the path integral approachcould also be an interesting computational tool to solve complex problems in quantitative finance. In thiscontext, this research could be important not only because it develops a novel and efficient technique for thesolution of the renowned CEV model, but also because it could open the door to computational applicationsof such methods of quantum mechanics. Indeed, our simulations show that computing the CEV optionpricing formula using the semiclassical approximation of path integral is efficient and accurate comparedto the standard CEV solution considering the pricing of European call options. Additionally, the proposedapproximation reduces execution times importantly and keep the simplicity of the traditional solution.Thus, the main idea of the paper is twofold, firstly to use ideas from quantum mechanics to deal withapplied finance problems, and secondly, to develop practical methodologies and to test them numericallyin specific case studies, while discussing their practical advantage and limitations. The structure of thepaper is the following. Firstly, Feynman’s path integral formulation is revisited. Secondly, the path integralapproximation is applied to the basic BS model. Thirdly, the path integral approximation is applied to For discrete-time approaches to modeling volatility, see refs. [6, 7] A comprehensive review for stochastic volatility models can be found in [8] and [9] . A.k.a. ‘local volatility’
The path integral formalism was developed by Richard P. Feynman [43], introducing the action principlefrom classical mechanics to quantum mechanics. Nowadays Feynman’s path integral is a well-known tool inquantum mechanics, and statistical and mathematical physics, with applications in many branches of physicssuch as: optics, thermodynamics, nuclear physics, atomic and molecular physics, cosmology, polymer scienceand other interdisciplinary areas [44, 45].In the following lines, we describe the fundamentals of the path integral methodology. The starting pointis the Schrödinger equation: i (cid:126) ∂ Ψ ∂ ˜ t = ˆH BS Ψ (1)where Ψ is the wave function and ˆH the Hamiltonian quantum operator (for this instance we consider a timeindependent Hamiltonian).Considering Ψ ( x ) as the initial value of Ψ (i.e., Ψ( x, t = 0) = Ψ ( x )), the general solution of 20 is givenin term of the unitary evolution operator:Ψ( x, t ) = e − i ˆH t / (cid:126) Ψ ( x ) (2)Equivalently, using convolution properties, the value at time t of the wave function is represented by:Ψ( x, t ) = ˆ ∞−∞ e − i ˆH t/ (cid:126) δ ( x − x ) Ψ ( x ) d x = ˆ ∞−∞ K ( x, t | x ,
0) Ψ ( x ) d x (3)where K ( x, t | x ,
0) = < x | e − i tˆH / (cid:126) | x > is called the propagator.Feynman concentrated on a previous work of Dirac [46], related with the proportionality between theexponential of the action over the classical path (which come from the Lagrangian formalism) and thepropagator in quantum mechanics: K ( x, t | x , ∝ e ( i/ (cid:126) ) A [ x cl ] where A is the action functional, defined as the time integral Lagrangian: A [ x ( t )] = ˆ t L ( x, ˙ x, t )d t A [ x cl ] indicates that the action is evaluated over the classical trajectory from x to x .Feynman reformulated Dirac formulation and described the propagator as the contributions of the allvirtual paths, not only the classical ones: K ( x, t | x ,
0) = X All Pathsfrom x to x ˜ N e ( i/ (cid:126) ) A [ x ( t )] N is an appropriated normalization for K .Thus, using the Riemann integral for each path (see ref. [44]), the propagator is defined as: K ( x, t | x ,
0) = ˆ D [ x ( t )]e ( i/ (cid:126) ) A [ x ( t )] (4)The functional integral of the right-hand side of the Eq. 4 is defined as a ‘Path Integral’, and the measureof the integration is given by D [ x ( t )] which means the integrations over all trajectories.The computation of the path integral is done via the time slicing scheme [43, 44], which is not a straight-forward procedure. Nevertheless, there is an alternative and popular method used in physics called the‘semiclassical approximation’, which approximates the argument of the path integral into a Gaussian func-tion, arriving this way to a solution in terms of the classical path, see [42, 47, 48]. Since the aim of the paper isto find a more efficient numerical solution to a complex problem, this avenue seems plausible and attractive,see [42] for a presentation of the semiclassical approximation to option pricing. The general procedure isexplained below.First, we write the path that links the points x ( t ) = x with x ( t ) = x as the classical trajectories asthe main contribution plus the fluctuations around it: x ( t ) = x ( t ) cl + δx ( t ) (5)with the fixed conditions (extremality condition): δx ( t ) = δx ( t ) = 0 (6)Later, we can expand the action around to x cl ( t ) using a functional Taylor series [49]: A [ x ( t ) cl + δx ] = A [ x ( t )] (cid:12)(cid:12)(cid:12)(cid:12) x cl ( t ) + ˆ t t d t δA [ x ( t )] δx ( t ) (cid:12)(cid:12)(cid:12)(cid:12) x cl ( t ) δx ( t ) (7)+ 12 ˆ t t d t d t δ Aδx ( t ) δx ( t ) δx ( t ) δx ( t ) (cid:12)(cid:12)(cid:12)(cid:12) x cl ( t ) (8)+ 13! ˆ t t d t d t d t δ Aδx ( t ) δx ( t ) δx ( t ) δx ( t ) δx ( t ) δx ( t ) (cid:12)(cid:12)(cid:12)(cid:12) x cl ( t ) + O (4) (9)The semiclassical approximation consist in truncated up to the quadratic terms the expansion 7: A [ x ( t )] ≈ A [ x cl ( t )] + 12 ˆ t t δ Aδx ( t ) δx ( t ) δx ( t ) δx ( t ) (cid:12)(cid:12)(cid:12)(cid:12) x cl ( t ) where the linear term is vanished due to the extremality condition.Thus, the propagator in the semiclassical limit becomes: K SC ( x , t | x , t ) = e ( i/ (cid:126) ) A [ x cl ( t )] ˆ xx T D [ χ ( t )] e ( i/ (cid:126) ) ´ t t δ Aδx ( t ) δx ( t ) δx ( t ) δx ( t ) (cid:12)(cid:12)(cid:12)(cid:12) xcl ( t ) = e ( i/ (cid:126) ) A [ x cl ( t )] N (10)where N is a normalization constant which incorporates the contribution of the second order term, definedby a Gaussian path integral . An analytical expression was developed for it in ref. [50] as the necessarycondition to maintain the unitary measure of the probability amplitudes [51], and it’s equal to:4 = r − M π (11)where M is the the van Vleck-Pauli-Morette determinant [50, 53], computed as: M = ∂ A class ∂y ∂y T (12)Finally, in the semiclassical regime, the propagator becomes : K ( x, t | x ,
0) = r − M π e i (cid:126) A [ x cl ] (13)The only necessary condition to get a solution for Eq. 13 is to have an analytical expression for the actionover the classical path. This can be achieved via the Hamilton equations (or Euler-Lagrange equation) usingthe classical Hamiltonian related to the quantum Hamiltonian defined in 1.Finally, two important notes must be considered in relation to the semiclassical approximation [45]:i) It is exact if the Lagrangian is quadratic.ii) It satisfy the Schrödinger equation up to terms of order (cid:126) .In the next section, we apply the semiclassical approximation of path integral to the European-vanilla typeoption pricing, arriving to the famous Black-Scholes model. We assume stochastic spot prices S t , governing by a standard geometric Brownian motion under the physical P -measure of the form: d S t S t = u d t + σ d ˜ W t (14)where W t is a standard Gauss-Wiener process with variance t . The parameters u and σ are the drift and thevolatility of the return, respectively. At this stage, we set these parameters as constants.Given the risk-free rate r , and defining the market price of risk: λ = µ − rσ we can describe the diffusion process under the unique risk-neutral measure ( Q -measure) instead of thephysical measure ( P -measure) using the Girsanov’s theorem (see [54] for a detail explanation). In short, wedefine a new Brownian motion under the Martingale measure of the form:d W t = λ d t + d ˜ W t and replacing into Eq. 14, the price dynamics is described under the risk neutral measure, and it is givenby : A.k.a Morette-Van Hove determinant. See ref. [52] for details The Eq. 13 is called the Pauli formula [45] Also called equivalent martingale measure (EMM) The Girsanov’s theorem ensure a equivalent measure in which W t is a Wiener process and S t is a martingale (risk-neutral) S t S t = r d t + σ d ˜ W t (15)By Itô’s calculus, is possible to rewrite the Eq. 14 into:d (ln S t ) = (cid:18) r − σ (cid:19) d t + σ d W t (16)and labeling x t = ln S t : d x t = (cid:18) r − σ (cid:19) d t + σ d W t (17)The probability density P ( x t , t, x , t ) for the random variable x t evolves according to the Fokker-Planck(or forward Kolmogorov) equation [55]: ∂P∂t = − ∂∂x (cid:20)(cid:18) r − σ (cid:19) P (cid:21) + 12 ∂ ∂x (cid:0) σ P (cid:1) = 12 σ ∂ P∂x − (cid:18) r − σ (cid:19) ∂P∂x (18)with initial condition: P ( x t , t = 0) = δ ( x )Using the following simple transformation: c = e − rt P and rewriting x in terms of S ( x = ln S ), the Eq. 18 yields to the Black-Scholes equation in it standard form[1]: ∂c∂t = 12 σ S ∂ c∂S + rs ∂c∂x − rc (19)Using the wick rotation (˜ t = it ), the evolution of the probability density P (Eq. 18) can be mapped tothe Schödringer equation: i (cid:126) ∂ Ψ ∂ ˜ t = ˆH BS Ψ (20)where the wave function Ψ represents the probability P , and the quantum Hamiltonian ˆH BS , namely for thisinstance the Black-Scholes Hamiltonian, is given by [56]:ˆH BS = 12 σ ∂ ∂x − (cid:18) r − σ (cid:19) ∂∂x In order to ensure the compatibility between the Eqs. 18 and 20, we need to set (cid:126) = 1.Given the momentum operator ˆp = − i (cid:126) ∂∂x = − i ∂∂x , the Hamiltonian can be expressed as :6H BS = − σ ˆp − i (cid:18) r − σ (cid:19) ˆpConsidering Ψ ( x ) as the initial value of Ψ (i.e., Ψ( x, t = 0) = Ψ ( x )), the general solution of 20 is givenby (see [37]): Ψ( x, t ) = e − i ˆ H BS ˜ t Ψ ( x )= e ˆ H BS t Ψ ( x )However, the known conditions in the option pricing context (i.e., contract function) is set at time T .Thus, defining the backward time τ = T − t and considering a final term value of the wave function Ψ( x, T ) =Ψ T ( x T ), the solution becomes: Ψ( x, t ) = e − ˆ H BS τ Ψ T ( x ) (21)Equivalently, using the convolution properties, the value at time t of the wave function is represented by:Ψ( x, t ) = ˆ ∞−∞ e − ˆ H BS τ δ ( x − x T ) Ψ T ( x T ) d x T = ˆ ∞−∞ K BS ( x, τ | x T ,
0) Ψ T ( x T ) d x T where K BS ( x, τ | x ,
0) is the propagator, which admits the following path integral representation in euclideantime [44]: K BS ( x, τ | x T ,
0) = ˆ D x ( τ ) e − S BS [ x ( τ )] being S BS [ x ( τ )] the euclidean classical action along all the paths x ( t ) which link the points x ( T ) = x T and x ( t ) = x ; defined by: S [ x ( t )] = ˆ Tt L BS dτ with L BS the Lagrangian.In order to obtain an expression for the propagator (Eqs. 4-5) we request the classical action evaluatedover the classical path. This can be obtained using the classical Hamiltonian mechanics.The classical Hamiltonian H BS associated to the operator ˆH BS is: H BS = − σ p − i (cid:18) r − σ (cid:19) p − i ˙ x = ∂ H BS ∂p − i ˙ p = − ∂ H BS ∂x or explicitly: p = iσ (cid:20) ˙ x − (cid:18) r − σ (cid:19)(cid:21) (22)˙ p = 0 (23)Then, the Lagrangian is given via the Legendre transformation: L BS = − ip ˙ x − H BS = − ip ˙ x + 12 σ p + i (cid:18) r − σ (cid:19) p = p (cid:20) σ p − i (cid:18) ˙ x − r + σ (cid:19)(cid:21) Using the values that solves the Hamilton’s equation (Eqs. 22-23), the Lagrangian is: L BS = 12 σ (cid:20) ˙ x − (cid:18) r − σ (cid:19)(cid:21) (24)Later, the Euler-Lagrange equation: dd t (cid:18) ∂ L BS ∂ ˙ x (cid:19) − ∂ L BS ∂x = 0 (25)yields to the free particle Newton equation: ¨ x = 0 (26)which leads to: ˙ x = C (27) x = Ct + D (28)8he values for C and D are obtained using the border conditions (fixed values) for x : x (0) = x T x ( τ ) = x Thus, the classical path, with 0 ≤ τ ≤ T , is described by:˙ x ( t ) = x T − x τ (29) x ( t ) = x T − x τ τ + x (30)Then, using 29 and 30, the corresponding classical action over the classical path: A [ x class ( t )] = ˆ τ σ (cid:20) ˙ x ( τ ) − (cid:18) r − σ (cid:19)(cid:21) d ( τ )= 12 σ τ (cid:20) x T − x − τ (cid:18) r − σ (cid:19)(cid:21) (31)Now, we are in conditions to compute the propagator. According to Eqs. 11 and 12, for this case: N = 1 p πσ ( T − t )and the semiclassical approximation for the propagator becomes: K SCBS ( x, τ | x T ,
0) = e − S BS [ x class ( τ )] p πσ ( T − t )= 1 √ πσ τ e − σ τ (cid:2) x T − x − τ (cid:0) r − σ (cid:1)(cid:3) Then, the wave function solution is reduced to:Ψ( x, t ) = 1 p πσ ( T − t ) ˆ ∞−∞ e − S BS [ x class ( τ )] Ψ T ( x T ) d x T (32)= 1 √ πσ τ ˆ ∞−∞ e − σ T − t ) (cid:2) x T − x − ( T − t ) (cid:0) r − σ (cid:1)(cid:3) Ψ T ( x T ) d x T (33)which is equal to the convolution between the propagator and the contract function:9 ( x, t ) = K BSSC ∗ Ψ T ( x T )The solution of Eq. 32 depends on the border condition Ψ T (contract function). We analyze the case ofan European call option, i.e.,: Ψ T ( x T ) = e − rτ max { S T − E, } = e − rτ max { e x T − E, } being E the strike price.Then, the wave function for this case is:Ψ( x, t ) = e − rτ √ πσ τ ˆ ∞ ln K e − σ τ (cid:2) x T − x − τ (cid:0) r − σ (cid:1)(cid:3) ( e x T − E ) d x T = e − rτ √ πσ τ ˆ ∞ ln E e x T e − σ τ (cid:2) x T − x − τ (cid:0) r − σ (cid:1)(cid:3) d x T − Ee − rτ √ πσ τ ˆ ∞ ln E e − σ τ (cid:2) x T − x − τ (cid:0) r − σ (cid:1)(cid:3) d x T = I − I Developing I , we have: I = e − rτ √ πσ τ ˆ ∞ ln E e − σ τ h x T − x T (cid:0) x + τ (cid:0) r − σ (cid:1)(cid:1) + (cid:0) x + τ (cid:0) r − σ (cid:1)(cid:1) i + x T d x T = e x √ πσ τ ˆ ∞ ln E e − σ τ (cid:2) x T − (cid:0) x + τ (cid:0) r + σ (cid:1)(cid:1)(cid:3) d x T Carrying out the change of variable u = h − x T + x + τ (cid:16) r + σ (cid:17)i / √ σ τ , and replacing x = ln S , wehave: I = − S " ˆ −∞ x − ln E + τ (cid:0) r + σ (cid:1) e − v d u = S N ( d )where N ( · ) is the standard normal cumulative function and d = ln ( S E ) + τ (cid:0) r + σ (cid:1) √ σ τ .For to solve I we use the change of variable v = − h x T − x − τ (cid:16) r − σ (cid:17)i / √ σ τ , so:10 = − Ee − rτ √ π ˆ −∞ x − ln E + τ (cid:0) r − σ (cid:1) e − v d u = Ke − rτ N ( d )being d = ln ( S E ) + τ (cid:0) r − σ (cid:1) √ σ τ = d − √ σ τ .Finally, the price of a call option at time t , using the path integral formulæis given by:Ψ ( S , t ) = S N ( d ) − Ke − rτ N ( d )which is exactly the same value obtained by Black-Scholes [1] for an European call option . In the CEV model, under the risk-neutral measure, the asset is governed by the following stochastic differ-ential equation [11, 12]: d S ( S, t ) = rS d t + σS α d W (34)being r the constant risk-free of interest, σ and α taking constant values and W a standard Wiener process,whit d W ∼ N (0 , d t ). In its paper, Cox imposed the domain for α in the range [0 , α <
0, the volatility unrealistically goes to zero as S increases [59]. Then, the sameCox’s condition for α is assumed in this paper.The process described by the Eq. 34 can be interpreted as a generalization of the standard geometricBrownian motion used in the Black-Scholes model [1], but considering a non-constant local volatility functionequals to σS α − . In fact, for the limit case α = 2, the Eq. 34 is degenerated to the BS case. Also, the CEVmodel has correspondence with other approaches: For α = 1, it becomes a square root process, addressed byCox and Ross [60]; and for α = 0 , S follows an Ornstein-Uhlenbeck type process [61].The CEV model described in Eq. 34 owes its name to the fact that the variance of the return is given by: v = var (cid:18) d SS (cid:19) = var (cid:16) r d t + σS α − d W (cid:17) = σ S α − d t and then, the elasticity of the variance with respect to the spot: As noted previously (On page 5), the semiclassical approximation is exact if the Lagrangian is quadratic, as in the case ofthe B-S Lagrangian (Eq. 24) v/v d S/S = α − y ( S, t ) = S − α and by the Itô’s Lemma, Eq. 34 can be rewrite as:d y = (2 − α ) (cid:20) ry + 12 (1 − α ) σ (cid:21) d t + (2 − α ) σ √ y d W The Fokker-Planck equation rules the transition probability P ( Y, t ) of the variable Y . Thus: ∂P∂t = 12 ∂ ∂y h (2 − α ) σ yP i − ∂∂y (cid:20) (2 − α ) (cid:18) ry + 12 (1 − α ) σ (cid:19) P (cid:21) = 12 β σ y ∂ P∂y + βr [ γ − y ] ∂P∂y − βrP (35)being β and γ constant values (parameters), defined as: β = 2 − αγ = 3 − α r σ The relationship 35 can be interpreted as the Schrödinger equation in Euclidean (Wick-rotated) time,with (cid:126) = 1: ∂ Ψ ∂t = ˆHΨwhere the wave function Ψ is equivalent to the probability P and the Hamiltonian operator ˆH is given by:ˆH = 12 β σ y ∂ ∂y + βr [ γ − y ] ∂∂y − βr p = − i ∂∂Y , the Hamiltonian goes to:ˆH = − β σ y ˆ p + iβr [ γ − y ] ˆ p − βr Later, we consider a final term condition (contract function) of the form:Ψ( y, t = T ) = Ψ( y T )The wave function Ψ, can be written in terms of it propagator K :Ψ ( y, t ) = ˆ ∞−∞ K ( y, τ | y T ,
0) Ψ( y T )d y T where τ = T − t , is the backward time and: K ( y, τ | y T ,
0) = < y | e − τ ˆH | y T > = e − ˆ Hτ δ ( y − y T )On the other hand, the propagator can be estimated using the path integral: K ( y T , T | y ,
0) = ˆ D [ y ( t )] e − S [ y ( t )] being D [ y ( τ )] the infinitesimal contribution of all the paths y ( τ ) that satisfies the boundary conditions y ( t = T ) = y T and y ( t = 0) = y ; and S the euclidean classical action functional over y ( t ).Using semiclassical arguments the propagator becomes: K ( y T , | y , τ ) = e − A [ y class ( t )] r π M The classical path is obtained as the solution of the Hamilton equations. The classical Hamiltonian H related to ˆH is: H = − β σ yp + iβr [ γ − y ] p − βr where p represents the classical momentum. Considering the Hamilton equation in Euclidean time, themomentum can be written in terms of y and ˙ y : p = i ˙ y + βr [ γ − y ] β σ y (36)So, using Eq. 36, the Lagrangian takes the form: 13 = − i ˙ yp − H (37)= { ˙ y + βr [ γ − y ] } β σ y + Ar The unique classical trajectory is which obeys the Euler-Lagrange equation:dd t (cid:18) ∂ L ∂ ˙ y (cid:19) − ∂ L ∂y = 0 (38)Computing the derivatives: ∂ L ∂ ˙ y = ˙ y + βr ( γ − y ) β σ y dd t (cid:18) ∂ L ∂ ˙ y (cid:19) = ¨ yy − ˙ y − βγr ˙ yβ σ y ∂ L ∂y = − [ ˙ y + βr ( γ + y )] [ ˙ y + βr ( γ − y )]2 β σ y and replacing into the Eq. 38, we have a second order differential equation that rules the classical behaviorof y ( t ): 2 y ¨ y − ˙ y + β r (cid:0) γ − y (cid:1) = 0 (39)Then, solving the Eq. 39, the classical path is given by: y class ( t ) = (cid:0) C + 2 C e − rtβ (cid:1) − γ C e − rtβ (40)being C and C constants given by the fixed values of the path at time t = 0 and t = T : y ( T ) = y T y ( t ) = y which yields to: C = (cid:0) e rτβ + 1 (cid:1) q γ ( e rτβ − + 4 y y T e rτβ − e rτβ ( y + y T )( e rτβ − (41) C = (cid:20) y T e rτβ + y − q γ ( e rτβ − + 4 y y T e rτβ (cid:21) ( e rτβ − (42)14ater, using Eq. 40, the Lagrangian over the classical path is: L class = L [ y class ]= r (cid:0) C + 2 C e rtβ + γ (cid:1) ( γ − C ) σ C e rtβ ( C + 2 C e rtβ − γ ) + βr (43)Thus the classical action is obtained by time integration of the Eq. 43 : A class = t = T ˆ t = t L class dt = rσ ( βσ t − γrt + 2 γβ ln (cid:2) γ − (cid:0) C + 2 C e rtβ (cid:1)(cid:3) + (cid:0) γ − C (cid:1) Aβ e rtβ )(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) t = Tt =0 = βrτ − γr σ τ + 2 γrβσ ln " γ − (cid:0) C + 2 C e rτβ (cid:1) γ − ( C + 2 C ) + (cid:0) γ − C (cid:1) βC e rτβ (cid:0) − e rτβ (cid:1) (44)So, using the Eq. 15, the van Vleck determinant (Eq. 12) is computed as: M = 2 γr [ γ − ( C + 2 C )] βσ [ γ − ( C + 2 C e rτβ )] ( − ∂ ∂x ∂x T (cid:0) C + 2 C e rτβ (cid:1) γ − ( C + 2 C ) − h ∂∂x ( C + 2 C ) i h ∂∂x T (cid:0) C + 2 C e rτβ (cid:1)i + h ∂∂x T ( C + 2 C ) i h ∂∂x (cid:0) C + 2 C e rτβ (cid:1)i [ γ − ( C + 2 C )] + 2 (cid:2) γ − (cid:0) C + 2 C e rτβ (cid:1)(cid:3) h ∂∂x T ( C + 2 C ) i h ∂∂x ( C + 2 C ) i [ γ − ( C + 2 C )] + (cid:2) γ − (cid:0) C + 2 C e rτβ (cid:1)(cid:3) h ∂ ∂x ∂x T ( C + 2 C ) i [ γ − ( C + 2 C )] ) − ( γr [ γ − ( C + 2 C )] h − ∂∂x (cid:0) C + 2 C e rτβ (cid:1)i βσ [ γ − ( C + 2 C e rτβ )] × − ∂∂x T (cid:0) C + 2 C e rτβ (cid:1) γ − ( C + 2 C ) + (cid:2) γ − (cid:0) C + 2 C e rτβ (cid:1)(cid:3) h ∂∂x T ( C + 2 C ) i [ γ − ( C + 2 C )] ) + 2 γr h − ∂∂x ( C + 2 C ) i βσ [ γ − ( C + 2 C e rτβ )] ( (cid:2) γ − (cid:0) C + 2 C e rτβ (cid:1)(cid:3) h ∂∂x T ( C + 2 C ) i [ γ − ( C + 2 C )] − ∂∂x T (cid:0) C + 2 C e rτβ (cid:1) γ − ( C + 2 C ) ) + r (cid:0) e βrτ − (cid:1) h(cid:16) ∂C ∂x (cid:17) (cid:16) ∂C ∂x T (cid:17) + ∂ C ∂x ∂x T i βσ C e rτβ r (cid:0) e βrτ − (cid:1) h(cid:16) ∂C ∂x T (cid:17) (cid:16) ∂C ∂x (cid:17) + (cid:16) ∂C ∂x (cid:17) (cid:16) ∂C ∂x T (cid:17) − (cid:0) γ − C (cid:1) (cid:16) ∂ C ∂x ∂x T (cid:17)i βσ C e rτβ − r (cid:0) e βrτ − (cid:1) (cid:0) γ − C (cid:1) (cid:16) ∂C ∂x T (cid:17) (cid:16) ∂C ∂x (cid:17) βσ C e rτβ Then is possible to compute the semiclassical propagator, through the Euclidean form of the Pauli’sformula (Eq. 5): K = e − A [ y class ( t )] r π M Finally, the value of the wave function at time t , is given by:Ψ ( y, t ) = r π ˆ ∞−∞ √M e − A class Ψ( y T )d y T Coming back to the option pricing problem, if we consider an European call option, with strike E andmaturity T , the value of the option at time t under the CEV model will be: C ( S, t ) = r π ˆ ∞ E / (2 − α ) √M e − A [ y class ( t )] (cid:0) y − αT − E (cid:1) d y T (45)which unfortunately is not possible to evaluate analytically, but it can be easily computed numerically forany conventional integration method. We compute numerically, using an standard method (global adaptive quadrature [62]), the integral definedin Eq. 45. We also compute the pricing for the same European call option using the Schroder approach [14]that consider the non-central chi-square distribution, and set it as the benchmark.We examine the results of both models, in terms of the pricing and the running time of each computation;considering several volatilities and elasticities of variances. Besides we test our results for short-time matu-rities ( T = { . , . } ) and long time maturities ( T = { , } ). In all the experiments we assume r = 0 . S = 100 and E = 110 ,Firstly, we consider a maturity equal to six months. In Table 1 both the pricing and computational timeare reported. We can see that the path integral method has similar pricing values but with a clear advantagein the running time. The times observed in the Table 1 for the proposed method of path integral are alwayslower to 0.008 seconds; however for the non-central chi-square approach the times are at least greater in oneorder of magnitude with a increase when the elasticity parameter is higher.16 α Path Integral BenchmarkPricing ($) Running Time (s) Pricing ($) Running Time (s)20% 1 4.4289e-08 0.0079 4.6567e-08 0.15951.45 0.0580 0.0064 0.0600 0.08861.9 1.8505 0.0060 1.8706 0.210150% 1 0.0259 0.0078 0.0583 0.02751.45 1.3437 0.0079 1.4181 0.04801.9 8.0777 0.0059 8.2636 0.060390% 1 0.3847 0.0059 0.4148 0.03071.45 3.9003 0.0077 4.2358 0.02361.9 16.4965 0.0074 17.1870 0.0413Table 1: Comparison of pricing and computational time for a Call option for some values of σ and α , using T = 0 . S = 100 and E = 110For a clearer and complete view, we present the continuous results in Figs. 1, 2 and 3 The pricing andthe running time are showed for both models sweeping on values of α . The figures confirm the observedconcussions in Table 1 in the sense that the running times of the proposed method of path integral aresignificantly lower (right-hand side figures) than the traditional solution methodology for the CEV model,especially when α tends to 2 where the time of the benchmark method rises considerably. In terms of theaccuracy, we can see that the path integral method fits very well in all cases. In order to have an estimationof the path integral approach, in the Fig. 4 the absolute and relative errors are shown for several values of α and σ . Always, the relative relative error is no longer that 10% for the assumed parameters. α P r i ce ( $ ) Call Price for the CEV model (for σ =0.2) Path integral ApproachBenchmark (a) Option Pricing α T i m e ( s ) Running time for a Call option pricing under the CEV model (for σ =0.2) Path integral ApproachBenchmark (b) Running time
Figure 1: Pricing and computational time for a Call option using σ =20%, T = 0 . r = 0 . S = 100 and E = 110 17 α P r i ce ( $ ) Call Price for the CEV model (for σ =0.5) Path integral ApproachBenchmark (a) Option Pricing α T i m e ( s ) Running time for a Call option pricing under the CEV model (for σ =0.5) Path integral ApproachBenchmark (b) Running time
Figure 2: Pricing and computational time for a Call option using σ =50%, T = 0 . r = 0 . S = 100 and E = 110 α P r i ce ( $ ) Call Price for the CEV model (for σ =0.9) Path integral ApproachBenchmark (a) Option Pricing α T i m e ( s ) Running time for a Call option pricing under the CEV model (for σ =0.9) Path integral ApproachBenchmark (b) Running time
Figure 3: Pricing and computational time for a Call option using σ =90%, T = 0 . r = 0 . S = 100 and E = 110 18 A b s o l u t e e rr o r α σ (a) Absolute error A b s o l u t e e rr o r α σ (b) Relative error Figure 4: Absolute and relative error of of the path integral approach with T = 0 . r = 0 . S = 100 and E = 110If we use a lower time to maturity (three months) the results in terms of computational time are verysimilar to the case T = 0 .
5, and the fit is still good too. In fact, for lower maturity the path integral methodperforms better because the error is no greater than 2%. This is showed from Fig. 5 to Fig. 8. α P r i ce ( $ ) Call Price for the CEV model (for σ =0.2) Path integral ApproachBenchmark (a) Option Pricing α T i m e ( s ) Running time for a Call option pricing under the CEV model (for σ =0.2) Path integral ApproachBenchmark (b) Running time
Figure 5: Pricing and computational time for a Call option using σ =20%, T = 0 . r = 0 . S = 100 and E = 110 19 α P r i ce ( $ ) Call Price for the CEV model (for σ =0.5) Path integral ApproachBenchmark (a) Option Pricing α T i m e ( s ) Running time for a Call option pricing under the CEV model (for σ =0.5) Path integral ApproachBenchmark (b) Running time
Figure 6: Pricing and computational time for a Call option using σ =50%, T = 0 . r = 0 . S = 100 and E = 110 α P r i ce ( $ ) Call Price for the CEV model (for σ =0.9) Path integral ApproachBenchmark (a) Option Pricing α T i m e ( s ) Running time for a Call option pricing under the CEV model (for σ =0.9) Path integral ApproachBenchmark (b) Running time
Figure 7: Pricing and computational time for a Call option using σ =90%, T = 0 . r = 0 . S = 100 and E = 110 20 A b s o l u t e e rr o r
1 1.7354 0.8 0.9 α σ (a) Absolute error R e l a t i ve e rr o r ( % ) α σ (b) Relative error Figure 8: Absolute and relative error of of the path integral approach with T = 0 . r = 0 . S = 100 and E = 110For greater times to maturity, we have a change in the results, indicating the limits of the semiclassicalapproximation. For a two years maturity Figs. 9-12 we find results very similar to that of the previous cases,but with a more deviation in the pricing (absolute error). Still, the relative error remains lower than 12%.However the figures show an increase in the running time, being comparable the times of both approaches,specially when the volatility rises and the elasticity is low. This fact is confirmed when we use a maturityequals to 4 years. Indeed, the computational cost increase, being the proposed method still competitive forhigher α . In the same way, the pricing error goes up, despite the fact that the relative error remains under20%. α P r i ce ( $ ) Call Price for the CEV model (for σ =0.2) Path integral ApproachBenchmark (a) Option Pricing α T i m e ( s ) Running time for a Call option pricing under the CEV model (for σ =0.2) Path integral ApproachBenchmark (b) Running time
Figure 9: Pricing and computational time for a Call option using σ =20%, T = 2, r = 0 . S = 100 and E = 110 21 α P r i ce ( $ ) Call Price for the CEV model (for σ =0.5) Path integral ApproachBenchmark (a) Option Pricing α T i m e ( s ) Running time for a Call option pricing under the CEV model (for σ =0.5) Path integral ApproachBenchmark (b) Running time
Figure 10: Pricing and computational time for a Call option using σ =50%, T = 2, r = 0 . S = 100 and E = 110 α P r i ce ( $ ) Call Price for the CEV model (for σ =0.9) Path integral ApproachBenchmark (a) Option Pricing α T i m e ( s ) Running time for a Call option pricing under the CEV model (for σ =0.9) Path integral ApproachBenchmark (b) Running time
Figure 11: Pricing and computational time for a Call option using σ =90%, T = 2, r = 0 . S = 100 and E = 110 22 A b s o l u t e e rr o r α σ (a) Absolute error R e l a t i ve e rr o r ( % ) α σ (b) Relative error Figure 12: Absolute and relative error of of the path integral approach with T = 2, r = 0 . S = 100 and E = 110 α P r i ce ( $ ) Call Price for the CEV model (for σ =0.2) Path integral ApproachBenchmark (a) Option Pricing α T i m e ( s ) Running time for a Call option pricing under the CEV model (for σ =0.2) Path integral ApproachBenchmark (b) Running time
Figure 13: Pricing and computational time for a Call option using σ =20%, T = 4, r = 0 . S = 100 and E = 110 23 α P r i ce ( $ ) Call Price for the CEV model (for σ =0.5) Path integral ApproachBenchmark (a) Option Pricing α T i m e ( s ) Running time for a Call option pricing under the CEV model (for σ =0.5) Path integral ApproachBenchmark (b) Running time
Figure 14: Pricing and computational time for a Call option using σ =50%, T = 4, r = 0 . S = 100 and E = 110 α P r i ce ( $ ) Call Price for the CEV model (for σ =0.9) Path integral ApproachBenchmark (a) Option Pricing α T i m e ( s ) Running time for a Call option pricing under the CEV model (for σ =0.9) Path integral ApproachBenchmark (b) Running time
Figure 15: Pricing and computational time for a Call option using σ =90%, T = 4, r = 0 . S = 100 and E = 110 24 A b s o l u t e e rr o r α σ (a) Absolute error R e l a t i ve e rr o r ( % )
15 0.91.7354 0.8 α
20 0.71.61405 σ (b) Relative error Figure 16: Absolute and relative error of the path integral approach with T = 4, r = 0 . S = 100 and E = 110 In this paper, a new numerical method for computing the CEV model was developed. In particular, thisnew approach was based on the semiclassical approximation of Feynman’s path integral model. This formu-lation dealt with some of the limitations of the conventional approach based on the non-central chi-squareddistribution.The experimental results showed a good fit between the new proposed method and the traditional method-ology (setting the former as benchmark), and also a lower computational cost, measured as the running timeof each model.We analyze several hypothetical scenarios, using different maturities, volatilities and elasticities. In mostcases, the running time is one order of magnitude lower than the benchmark, but if the elasticity tends to one,this difference is higher. As an accuracy measure the absolute and relative error are computed. For the range10% < σ < . < α < .
97 the relative error is below than 20% in all the cases. Nevertheless, forshort maturities and lower volatilities, the error decreases considerably, coming to be less than 10% for smallmaturities (under 2% for T=0.25!) and for σ < eferences [1] Fischer Black and Myron Scholes. The pricing of options and corporate liabilities.
Journal of PoliticalEconomy , 81(3):637–654, 1973.[2] Geert Bekaert and Guojun Wu. Asymmetric volatility and risk in equity markets.
The review of financialstudies , 13(1):1–42, 2000.[3] Tim Bollerslev, Julia Litvinova, and George Tauchen. Leverage and volatility feedback effects in high-frequency data.
Journal of Financial Econometrics , 4(3):353–384, 2006.[4] Mark Rubinstein. Nonparametric tests of alternative option pricing models using all reported trades andquotes on the 30 most active cboe option classes from august 23, 1976 through august 31, 1978.
TheJournal of Finance , 40(2):455–480, 1985.[5] Jens Carsten Jackwerth and Mark Rubinstein. Recovering probability distributions from option prices.
The Journal of Finance , 51(5):1611–1631, 1996.[6] Robert F Engle. Autoregressive conditional heteroscedasticity with estimates of the variance of unitedkingdom inflation.
Econometrica: Journal of the Econometric Society , pages 987–1007, 1982.[7] Tim Bollerslev, Ray Y Chou, and Kenneth F Kroner. Arch modeling in finance: A review of the theoryand empirical evidence.
Journal of econometrics , 52(1-2):5–59, 1992.[8] Eric Ghysels, Andrew C. Harvey, and Eric Renault. Stochastic volatility. In G.S. Maddala and C.R.Rao, editors,
Statistical Methods in Finance , volume 14 of
Handbook of Statistics , chapter 5, pages 119– 191. Elsevier, 1996.[9] Neil Shephard and Torben G Andersen. Stochastic volatility: origins and overview. In Thomas Mikosch,Jens-Peter Kreiß, Richard A. Davis, and Torben Gustav Andersen, editors,
Handbook of Financial TimeSeries , pages 233–254. Springer, Berlin, Heidelberg, 2009.[10] David G Hobson and Leonard CG Rogers. Complete models with stochastic volatility.
MathematicalFinance , 8(1):27–48, 1998.[11] John C Cox. Notes on option pricing i: Constant elasticity of variance diffusions. Working paper,Stanford University, 1975.[12] John C Cox. The constant elasticity of variance option pricing model.
The Journal of Portfolio Man-agement , 23(5):15–17, 1996.[13] Stan Beckers. The constant elasticity of variance model and its implications for option pricing.
TheJournal of Finance , 35(3):661–673, 1980.[14] Mark Schroder. Computing the constant elasticity of variance option pricing formula. the Journal ofFinance , 44(1):211–219, 1989.[15] Dmitry Davydov and Vadim Linetsky. Pricing and hedging path-dependent options under the cevprocess.
Management science , 47(7):949–965, 2001.[16] Min Chen. Equity volatility smile and skew under a cev-based structural leverage model.
AppliedMathematics & Information Sciences , 9(3):1095–1104, 2015.[17] James D MacBeth and Larry J Merville. Tests of the black-scholes and cox call option valuation models.
The Journal of Finance , 35(2):285–301, 1980.[18] David C Emanuel and James D MacBeth. Further results on the constant elasticity of variance calloption pricing model.
Journal of Financial and Quantitative Analysis , 17(4):533–554, 1982.2619] Alan L Tucker, David R Peterson, and Elton Scott. Tests of the black-scholes and constant elasticity ofvariance currency call option valuation models.
Journal of Financial Research , 11(3):201–214, 1988.[20] Beni Lauterbach and Paul Schultz. Pricing warrants: An empirical study of the black-scholes model andits alternatives.
The Journal of Finance , 45(4):1181–1209, 1990.[21] VK Singh and N Ahmad. Forecasting performance of constant elasticity of variance model: Empiricalevidence from india.
International Journal of Applied Economics and Finance , 5(1):87–96, 2011.[22] KC Yuen, H Yang, and KL Chu. Estimation in the constant elasticity of variance model.
BritishActuarial Journal , 7(2):275–292, 2001.[23] Ying-Lin Hsu, TI Lin, and CF Lee. Constant elasticity of variance (cev) option pricing model: Integrationand detailed derivation.
Mathematics and Computers in Simulation , 79(1):60–71, 2008.[24] Manuela Larguinho, José Carlos Dias, and Carlos A Braumann. On the computation of option pricesand greeks under the cev model.
Quantitative Finance , 13(6):907–917, 2013.[25] Nawdha Thakoor, Désiré Yannick Tangman, and Muddun Bhuruth. A new fourth-order numericalscheme for option pricing under the cev model.
Applied Mathematics Letters , 26(1):160–164, 2013.[26] Nawdha Thakoor, Désiré Yannick Tangman, and Muddun Bhuruth. Fast valuation of cev americanoptions.
Wilmott , 2015(75):54–61, 2015.[27] Richard Lu and Yi-Hwa Hsu. Valuation of standard options under the constant elasticity of variancemodel.
International Journal of Business and Economics , 4(2):157, 2005.[28] Rajae Aboulaich, Abdelilah Jraifi, and Ibtissam Medarhri. Stochastic runge kutta methods with the con-stant elasticity of variance (cev) diffusion model for pricing option.
International journal of mathematicalanalysis , 8(18):849–856, 2014.[29] Aricson Cruz and José Carlos Dias. The binomial cev model and the greeks.
Journal of Futures Markets ,37(1):90–104, 2017.[30] AE Lindsay and DR Brecher. Simulation of the cev process and the local martingale property.
Mathe-matics and Computers in Simulation , 82(5):868–878, 2012.[31] Sang-Hyeon Park and Jeong-Hoon Kim. Asymptotic option pricing under the cev diffusion.
Journal ofMathematical Analysis and Applications , 375(2):490–501, 2011.[32] Stefano Pagliarani and Andrea Pascucci. Analytical approximation of the transition density in a localvolatility model.
Central European Journal of Mathematics , 10(1):250–270, 2012.[33] Vladislav Krasin, Ivan Smirnov, and Alexander Melnikov. Approximate option pricing and hedging inthe cev model via path-wise comparison of stochastic processes.
Annals of Finance , pages 1–15, 2017.[34] Vadim Linetsky. The path integral approach to financial modeling and options pricing.
ComputationalEconomics , 11(1):129–163, 1997.[35] Guido Montagna, Oreste Nicrosini, and Nicola Moreni. A path integral way to option pricing.
PhysicaA: Statistical Mechanics and its Applications , 310(3):450–466, 2002.[36] G Bormetti, Guido Montagna, N Moreni, and Oreste Nicrosini. Pricing exotic options in a path integralapproach.
Quantitative Finance , 6(1):55–66, 2006.[37] Belal E Baaquie.
Quantum finance: Path integrals and Hamiltonians for options and interest rates .Cambridge University Press, 2007. 2738] D Lemmens, M Wouters, J Tempere, and S Foulon. Path integral approach to closed-form optionpricing formulas with applications to stochastic volatility and interest rate models.
Physical Review E ,78(1):016101, 2008.[39] Jeroen PA Devreese, Damiaan Lemmens, and Jacques Tempere. Path integral approach to asian optionsin the black–scholes model.
Physica A: Statistical Mechanics and its Applications , 389(4):780–788, 2010.[40] Mauricio Contreras, Rely Pellicer, Marcelo Villena, and Aaron Ruiz. A quantum model of option pricing:When black–scholes meets schrödinger and its semi-classical limit.
Physica A: Statistical Mechanics andits Applications , 389(23):5447–5459, 2010.[41] Mauricio Contreras. Stochastic volatility models at ρ = ± Physica A: Statistical Mechanics and its Applications , 405:289–302, 2014.[42] Zura Kakushadze. Path integral and asset pricing.
Quantitative Finance , 15(11):1759–1771, 2015.[43] Richard Phillips Feynman. Space-time approach to non-relativistic quantum mechanics.
Reviews ofModern Physics , 20(2):367, 1948.[44] Richard P Feynman and Albert R Hibbs.
Quantum mechanics and path integrals . Courier Corporation,Mineola, NY, 2010.[45] Christian Grosche and Frank Steiner.
Handbook of Feynman Path Integrals , volume 145 of
SpringerTracts in Modern Physics . Springer, 1998.[46] Paul A. M. Dirac. The lagrangian in quantum mechanics.
Phys. Z. Sowjetunion , 3(64):64–72, 1933.[47] Rf Rajaraman. Some non-perturbative semi-classical methods in quantum field theory (a pedagogicalreview).
Physics Reports , 21(5):227–313, 1975.[48] Thijs Koeling and RA Malfliet. Semi-classical approximations to heavy ion scattering based on thefeynman path-integral method.
Physics Reports , 22(4):181–213, 1975.[49] Masud Chaichian and A Demichev.
Path Integrals in Physics: Stochastic Processes and QuantumMechanics . Institute of Physics, 2001.[50] Cécile Morette. On the definition and approximation of feynman’s path integrals.
Physical Review ,81(5):848, 1951.[51] Cécile DeWitt-Morette. Feynman path integrals.
Communications in Mathematical Physics , 37(1):63–81, 1974.[52] Ph Choquard and F Steiner. The story of van vleck’s and morette van hove’s determinants.
HelveticaPhysica Acta , 69(5-6):636–654, 1996.[53] John H Van Vleck. The correspondence principle in the statistical interpretation of quantum mechanics.
Proceedings of the National Academy of Sciences , 14(2):178–188, 1928.[54] Rangarajan K Sundaram. Equivalent martingale measures and risk-neutral pricing: An expository note.
The Journal of Derivatives , 5(1):85–98, 1997.[55] Hannen Risken.
The Fokker–Planck equation , volume 18 of
Springer series in synergetics . Springer, firstedition, 1984.[56] Belal E Baaquie, Claudio Coriano, and Marakani Srikant. Hamiltonian and potentials in derivative pric-ing models: exact results and lattice simulations.
Physica A: Statistical Mechanics and its Applications ,334(3):531–557, 2004. 2857] Dirk Veestraeten. On the multiplicity of option prices under cev with positive elasticity of variance.
Review of Derivatives Research , 20(1):1–13, 2017.[58] Freddy Delbaen and Hiroshi Shirakawa. A note on option pricing for the constant elasticity of variancemodel.
Asia-Pacific Financial Markets , 9(2):85–99, 2002.[59] Matthew Lorig. The exact smile of certain local volatility models.
Quantitative Finance , 13(6):897–905,2013.[60] John C Cox and Stephen A Ross. The valuation of options for alternative stochastic processes.
Journalof financial economics , 3(1-2):145–166, 1976.[61] George E Uhlenbeck and Leonard S Ornstein. On the theory of the brownian motion.
Physical review ,36(5):823, 1930.[62] Lawrence F Shampine. Vectorized adaptive quadrature in matlab.