An analytical perturbative solution to the Merton Garman model using symmetries
AAn analytical perturbative solution to the MertonGarman model using symmetries
Xavier Calmet a and Nathaniel Wiesendanger Shaw a,ba School of Mathematical and Physical Sciences, University of Sussex, Brighton, BN1 9QH,United Kingdom b Business School, University of Sussex, Falmer, Brighton BN1 9RH, United Kingdom
Abstract
In this paper, we introduce an analytical perturbative solution to the Merton Garmanmodel. It is obtained by doing perturbation theory around the exact analytical solutionof a model which possesses a two-dimensional Galilean symmetry. We compare ourperturbative solution of the Merton Garman model to Monte Carlo simulations and findthat our solutions performs surprisingly well for a wide range of parameters. We alsoshow how to use symmetries to build option pricing models. Our results demonstratethat the concept of symmetry is important in mathematical finance.
Acknowledgments: The work of XC is supported in part by the Science and Technol-ogy Facilities Council (grant number ST/P000819/1). XC is very grateful to MITP fortheir generous hospitality during the academic year 2017/2018 where most of this work wascompleted. The authors would like to thank Andreas Kaeck for useful discussions.JEL Classification: C02, G10. Corresponding author, e-mail: [email protected], phone number: +44 (0)1273 877029, Fax: +44(0)1273 678097. E-mail: [email protected] a r X i v : . [ q -f i n . P R ] S e p Introduction
Calculating the price of an option is an important challenge in mathematical finance. Thefirst attempts in that direction are attributed to Louis Bachelier who during in his Doctoralthesis, Th´eorie de la sp´eculation, published in 1900, considered a mathematical model ofBrownian motion and its use for valuing options. This work provided the foundations forthe Black Scholes model (Black, Scholes (1973)). However, while the Black Scholes modelwas a breakthrough in the field, it is widely accepted that it has limitations. In particular,the volatility is treated as a constant which is not very realistic.Since the seminal works of Black, Scholes (1973) and Merton (1973), more sophisticatedmodels with a time dependent volatility have been proposed. For example, the affine Hestonmodel (See Heston (1993)), which assumes a time-dependent volatility, with a stochasticprocess involving the square-root of the stochastic volatility, and a leverage effect, has beenimplemented in a large number of empirical studies (Andersen et al. (2002), Bakshi et al.(1997), Bates (2000), Bates (2006), Chernov et al. (2003), Huan and Wu (2004), Pan (2002),Eraker (2004) to name a few). Such models have however limitations and are often modifiedartificially by combining them with models of jumps in returns and/or in volatility (suchas in Jones (2003) and Benzoni (2002)). As a consequence, there is a substantial strandof literature devoted to non-affine volatility models, which note that the popular square-root stochastic volatility model is not very realistic (see e.g. (Eraker et al. (2003) Duanand Yeh (2010), A¨ıt-Sahalia, Kimmel (2007) Christoffersen et al. (2010), Chourdakis andDotsis (2011) and Kaeck and Alexander (2012)) to name a few). However, the issue withsuch models is a general lack of closed form characteristic function, which makes pricingmuch more challenging. As stated in Chourdakis and Dotsis (2011) when regarding theplace of non-affine models and the debate of their tractability against affine models: “doesanalytically tractability come at the cost of empirical misspecification?”. It is a usefulendeavor to study non-affine model as we propose in this paper, if an analytical solution forthe option pricing formula can be found.A well-known example of such models is the Merton Garman model (Garman (1976)and Merton (1973)) which is indeed a more realistic model as it allows for a time-dependentvolatility and it is not restricted to an affine model for the volatility. However, solving non-affine models is time consuming, as it involves numerical methods. Thus, many practitionersare still using the Black Scholes formula to obtain a fast, albeit not necessarily very reliable,price quote for an option.The aim of our work is twofold. We will derive an analytical approximative solution tothe partial differential equation describing the Merton Garman model which enables one fast1alculations of option prices. This requires us to identify a “symmetric” version of the modelwhich can easily be solved analytically. One can then reintroduce the symmetry breakingterms of the original Merton Garman differential equation and do perturbation theory aroundthe symmetric solution thereby obtaining an approximative but analytical solution to theoriginal Merton Garman differential equation. We then propose a new approach to modelbuilding in option pricing based on the concept of symmetry groups and representationtheory. This concept has been extremely successful in modern physics. It is at the originof all successful models in physics, e.g., in particle physics, cosmology or solid state physics.We note that perturbation theory has been used in option pricing models (Baaquie (1997),Baaquie (2003), Blazhyevskyi and Yanishevsky (2011), Aguilar (2017), Kleinert and Korbel(2016), Utama and Purgon (2016)) but here we organize perturbation theory around a veryspecific solution, namely that of the symmetrical model which we will introduce in this paper.This paper is organized as follows. In section 2, we derive the partial differential equa-tion which describes the Merton Garman model. In section 3, we explain how to reduce theoriginal Merton Garman to a simple, symmetrical, model. We present an exact analyticalsolution to the symmetrical model. We then restore the original Merton Garman by rein-troducing the symmetry breaking terms and provide an analytical perturbative solution tothe Merton Garman model. In section 4, we compare our solutions to different numericalsolutions found in the literature. In section 5, we propose a new approach to model buildingin mathematical finance. Finally, we conclude in section 6.
In the Merton Garman model, the price of an option is dependent on the time t , the priceof the underlying S and the volatility V . Both S and V are taken to be time-dependentfunctions and thus the Merton Garman model has the potential to provide a more accuratecalculation of an option price than e.g. the Black Scholes model.We start from the stochastic differential equations for the price of the underlying S andfor the volatility V dS = rSdt + √ V SdW S , (1) dV = κ ( θ − V ) dt + ξV α dW V , (2)which resembles a stochastic, mean reverting, volatility regime. Here, ξ is the standarddeviation of the volatility and κ is the speed of mean reversion to the long run variance θ . The interest rate r is assumed to be constant. The model described by Eqs (1) and (2)covers many well-known stochastic volatility models, for instance setting α = 1 , / α can take arbitrary values.We will denote the correlation between the two Brownian motions W S and W V by ρ .We shall first consider a call option, but our results can be extended to a put option in astraightforward manner. Our first step is to find the associated partial differential equationwhich describes this model. We do so by applying the Feynman-Kac formula, see e.g Hull(1997), which states that for the price of a call option, as defined by the model dynamics inEqs (1) and (2) is given by: ∂C∂t + (cid:88) i =1 µ i ( t, x ) ∂C∂x i + 12 (cid:88) i,j =1 ρ ij σ i ( t, x ) σ j ( t, x ) ∂ C∂x i ∂x j − rC = 0 , (3)where C is the price of a call.Using this formula, one obtains ∂C∂t + rS ∂C∂S + 12 V S ∂ C∂S + ( λ + µV ) ∂C∂V + ρξV / α S ∂ C∂S∂V + ξ V α ∂ C∂V − rC = 0 , (4)where λ = κθ , µ = − κ . The call price C = C ( S, V, t ) depends on the time t , the priceof the underlying S and the time dependent volatility V = V ( t ). In this model, there arethree free parameters λ , µ and α . As we explained previously, existing solutions to thispartial differential equation are numerical ones which have been obtained using Monte Carlomethods. Note that the put price P = P ( S, V, t ) fulfills the same differential equation, butit is obviously subject to a different boundary condition.
By studying the partial differential equation given in Eq. (4), it quickly becomes clear thatthe difficulty in finding an analytical solution to Eq. (4) is due to the lack of symmetrybetween the different terms of the partial differential equation. It is useful to study thedimensions of the different terms and constants in this partial differential equation. Theprice of the call is obviously given in a specific currency which we shall take to be the USDor $. The remaining dimensions follow from this. We have: • [ C ] = $ • [ ∂C∂t ] = $ / time 3 [ rS ∂C∂S ] = [ r ]$ thus [ r ] = 1 / time • [ V S ∂ C∂S ] = [ V ]$ thus [ V ] = 1 / time • [( λ + µV ) ∂C∂V ] = ([ λ ] + [ µ ]1 / time) thus [ λ ] = 1 / time and [ µ ] = time • [ ρξV / α S ∂ C∂S∂V ] = [ ρ ][ ξ ](1 / time) / α $ time = $ / time thus [ ρ ][ ξ ] = time α − / • [ ξ V α ∂ C∂V ] = [ ξ ] (1 / time) α − $ = $ / time thus [ ξ ] = time α − / and ρ is dimensionless.It is instructive to see that S and V have different dimensions. Nevertheless, our goal isto treat S and V as symmetrically as possible to make a global Galilean invariance in 2+1manifest (see section 5). This can be achieved by adequate variable transformations and byidentifying the terms in the differential equation that violate this symmetry. Our aim is to derive a differential equation that is symmetrical in S and V . With this aimin mind, let us introduce an averaged volatility σ which is constant. As in the case of theBlack Scholes model, different definitions for the averaged volatility are possible, the specificchoice will not impact our methodology and results.Inspecting the differential equation (4), it is clear that we need to pick α = 1 to emphasizethe symmetry between S and V . We thus consider ∂C∂t + rS ∂C∂S + 12 σ S ∂ C∂S + µV ∂C∂V + ρξ V / S ∂ C∂S∂V + ξ V ∂ C∂V = rC. (5)We need to keep in mind that we will need to reintroduce V S ∂ C∂S , λ ∂C∂V and the termscorresponding to deviations from 1 for α . Note that ξ is different from ξ , in particular theydo not have the same dimensions. Finally, we see that there is a mixed derivative term whichneeds to be eliminated. We thus set ρ = 0 and we will reintroduce this term as symmetrybreaking term. We thus end up with: ∂C∂t + rS ∂C∂S + 12 σ S ∂ C∂S + µV ∂C∂V + ξ V ∂ C∂V = rC. (6)This partial differential equation can be massaged with standard substitutions into a 2+1dimensional heat equation (see Appendix A) in which case the symmetry in S and V becomesmanifest. In order to do so, we introduce x = log( S/K ) , (7) y = log( V /V ) , (8)4nd C ( x, y, τ ) = Kφ ( x, y, τ ) ψ ( x, y, τ ) , (9)where K is the strike price and V is some constant with units of 1 / sec. The function φ ( x, y, τ )and the rescaled time τ are defined in Appendix A. Standard manipulations described inAppendix A lead to ∂ψ ∂τ = ∂ ψ ∂x + ∂ ψ ∂y . (10)which is manifestly symmetrical in x and y . We will thus refer to the model described bythe differential equation (31) as the symmetrical model. Another reason for massaging thesymmetrical model into a heat equation is that this equation is easy to solve analytically.We impose the standard boundary condition for the call price: C ( S, V, T ) = (cid:32) S ( T ) − K (cid:33) + . (11)For a put option we have P ( S, V, T ) = (cid:32) K − S ( T ) (cid:33) + . (12) Details of the derivation of the analytical solution of the symmetrical model, i.e., of the 2+1dimensional heat equations, are given in Appendix B. We find C ( S, V, t ) = S N ( d ) − Ke − r ( T − t ) N ( d ) , (13)where N ( d ) = 1 √ π (cid:90) d −∞ exp (cid:18) − z (cid:19) dz, (14)and d = x √ τ + √ τ R + 1) = log( S/K ) + ( r + σ / T − t ) σ √ T − t , (15) d = x √ τ + √ τ R −
1) = d − σ √ T − t. (16)Remarkably, because of the boundary condition that only depends on S , it is identical tothe Black Scholes solution. 5or a put option, the very same procedure leads to P ( S, V, t ) = C ( S, V, t ) − S + Ke − r ( T − t ) . (17)In the next subsection, we shall restore the symmetry breaking terms and discuss the fullMerton Garman model. We are now in a position to solve the full Merton Garman model using perturbation theoryaround the symmetrical solution C ( S, V, t ). We organize perturbation theory as an expan-sion in terms the coefficients of the symmetry breaking terms. We first need to restore thefull model by re-introducing the symmetry breaking terms ∂C∂t + rS ∂C∂S + 12 σ S ∂ C∂S + c S (cid:32) V − σ (cid:33) ∂ C∂S + µV ∂C∂V + c λ ∂C∂V + ξ V ∂ C∂V + c (cid:32) ξ V α − ξ V (cid:33) ∂ C∂V + c ρξV α +1 / S ∂ C∂S∂V − rC = 0 . (18)Note that we have introduced dimensionless coefficients c i which denote the strength of thesymmetry breaking terms. In the limit c i = 1 one recovers the original Merton Garmanmodel. These coefficients are simply introduced as a bookkeeping trick to keep track ofwhich terms correspond to a deviation of the 2+1 Galilean invariant theory. In the end ofthe day, we set c i = 1. We now do perturbation theory around the symmetrical solution C ( S, V, t ) and obtain C ( S, V, t ) = − K (cid:0) SK (cid:1) − rσ e ( SK ) + ( r + σ ) t − T )28 σ t − T ) √ π (cid:16) √ γσξ + 1 (cid:17) (cid:112) σ ( T − t ) (19) × (cid:32) σ (cid:32) √ γσξ + 1 (cid:33) ( t − T ) + V (cid:18) e σ (cid:16) √ γσξ +1 (cid:17) ( T − t ) − (cid:19)(cid:33) , where we have set c = 1. We expect that our approximation should work well when λ and ρ are small, when α is close to one and when the variation of V around is average value σ is not too large. In the limit when V is large, σ is large as well and we expect that, as inthe Black Scholes case, the price of the call becomes the price of the underlying S . Detailsof the derivation can be found in Appendix C.6t may appear surprising that the leading order correction does not depend on the sym-metry breaking terms parametrized by c , c and c . It can easily be shown (see AppendixC) that the boundary condition (43) insures that only the contribution from the c termsurvives. The boundary condition implies that the contributions of c , c and c vanish toleading order in the perturbation theory. These symmetry breaking terms will, however,contribute to higher order corrections. Higher precision, if required, can be obtained bygoing to higher order in perturbation theory. Option prices can be calculated extremelyrapidly using this formalism. Note that, in principle, if we resummed perturbation theoryto all order in c i , the dependence on σ and ξ would vanish. It is also worth noticing thatour results are independent on V which is only introduced to match the dimension of V .It is straightforward to show that we obtain the same result for a put option P ( S, V, t ) = − K (cid:0) SK (cid:1) − rσ e ( SK ) + ( r + σ ) t − T )28 σ t − T ) √ π (cid:16) √ γσξ + 1 (cid:17) (cid:112) σ ( T − t ) (20) × (cid:32) σ (cid:32) √ γσξ + 1 (cid:33) ( t − T ) + V (cid:18) e σ (cid:16) √ γσξ +1 (cid:17) ( T − t ) − (cid:19)(cid:33) . The prices obtained to leading order in perturbation theory for a call and put option arethus given by C ( S, V, t ) = C ( S, V, t ) + C ( S, V, t ) , (21)and P ( S, V, t ) = P ( S, V, t ) + P ( S, V, t ) . (22)In the next section we shall compare our results to exact solutions obtained numerically. In this section we investigate how the approximative solution compares to a Monte Carlosimulation of the full Merton Garman model. It is well-known that the Merton Garmanmodel is not solvable analytically, however it can be solved using numerical methods. Thefirst step is a comparison of static cross sections of options. Then we compare using asimulated time series calibration exercise. 7 .1 Static Cross-Section Comparison
The first step in evaluating the performance of the leading order perturbative solution isto compare to a multitude of simulated data of the Merton Garman model to ensure thatthe approximation is sufficient to fit a range of different options at one time. We start bydescribing the data simulation of the Merton Garman model, then move onto the calibrationprocedure for the approximative solution and discuss the results.We choose a standard Monte Carlo framework, using stratified sampling and antitheticvariables, simulating seven million paths with a time step of one-tenth of a day for optionmaturities. We choose a spot underlying price of S = $100 and a strike range of K ∈ [90 , K/S ∈ [0 . , .
1] to simulate call options throughoutthe spectrum of moneyness . We simulate the Merton Garman model with the structuralparameter vector: Θ MG = { . , . , . , − . , } and initial volatility V ( t = 0) ∈ [10 , c , c and c . The solution is thus independent of the parameters: ρ, θ . However, as a by-product of the perturbation theory we introduced the followingparameters: ξ , σ . From inspection we fix ξ using ξ = ξσ α − which guarantees that ithas the right dimensions. While σ remains to be determined by calibration. Yielding theparameter vector: Θ pert. = { κ, ξ, α, σ } .The parameter σ is determined by calibration. When fitting σ to these simulations, thereis a risk of overfitting expensive out of the money (OTM) options. For that reason, it isbest to consider the implied volatility objective function (this is noted in Christoffersen etal. (2014)) which is given by: IV RM SE = (cid:118)(cid:117)(cid:117)(cid:116) N N (cid:88) i =1 ( IV MC i − IV pert. i ) , (23)where IV MC i stands for the implied volatility of the i th -option simulated using the MonteCarlo and IV pert. i is the implied volatility of the i th -option calculated using the leadingorder perturbative solution. The calibration exercise is extremely fast as our formula foroption prices is an approximative analytical solution. Figure (1) and Table 1 demonstratesthe results of the static calibration exercise for a 30 day maturity horizon . It should benoted that in Figure (1) the price Panels are the difference of the log, this is needed to In this exercise we simulate call options only, as put prices can be calculated from the put-call parity. We also simulate for α ∈ [0 . , .
5] and κ ∈ [1 . , , ρ ∈ [ − . , − . , θ ∈ [0 . , . While we simulated time horizons between 5-100 days, this is representative of our results and we dropthe other results for brevity.
K/S < C MC /C pert. ) ∼
0. Also apparent is that at some moneyness level theleading order perturbative solution will over price options, the level of moneyness at whichthis occurs is inversely proportional to volatility. Lastly, the range of under to over pricingis also inversely proportional to volatility. However, as the prices from both methods arevery similar it is more informative to look at the implied volatility cross-section in the evenPanels of Figure (1) along with the IVRMSE column of Table 1. These show throughoutthe range of volatility regimes the leading order perturbative solution is able to approximatesuccessfully a low IVRMSE, with a maximum occurring from the low volatility regime (Panel8) of 1 . C / ( C + C ) as a function of moneyness, K/S , where C is the contribution to the priceof the symmetrical solution and C is the leading order correction in perturbation theory.Figure (2) shows these ratios for several cases. Clearly C (cid:28) C even when the volatilityis large. This demonstrates nicely the validity of the perturbative expansion even in thecase where the volatility is large. While this exercise confirms that for fixed scenarios theperturbative solution is a very good approximation to the actual Merton Garman solution,it is essentially a multiple curve fitting exercise, a more thorough analysis is needed to beable to gauge the reliability of the perturbative approximation. This is what we shall focuson next. The second step in evaluating the performance of the leading order perturbative solutionis to estimate it against a time series simulation of the Merton Garman model. The timeseries uses 100 different Monte Carlo paths to simulate the asset price and variance pathsincluding 6 unique maturities within [7 , θ, ρ has on performance, as the other parameters in the Θ pert. vector will haveto attempt to absorb the information contained in the absent parameters. Secondly, it alsoprovides a first estimate into the applicability of the perturbative solution to different typesof options markets. We simulate four different data sets with varying parameter vectors,described below. • data set 1 : Θ MG = { . , . , . , − . , . } , with a negative correla-9ion which is a reasonable choice for modeling equity options, such as the S&P 500. • data set 2 : Θ MG = { . , . , . , . , . } , with a correlation of zero. • data set 3 : Θ MG = { . , . , . , +0 . , . } , with a positive correla-tion coefficient this is used to gain inference about modeling VIX options, see Park(2015). • data set 4 : Θ MG = { . , . , . , − . , . } , with a high(er) centraltendency we investigate an equity style option market with a central tendency whichis significantly higher than the initial variance value of: 0 .
08 and thus tests how theleading order perturbative solution handles significant change in the variance path.For the following simulated calibration exercises, it is imperative to note the difference inIVRMSE and parameters between the data sets as this will highlight the following: firstly,parameter regions where the leading order perturbative solution might breakdown. Sec-ondly, potential difficulties in estimating certain parameters of the model. Thirdly, whereinformation contained in θ, ρ might be absorbed. These results can be found in Tables 2-3.Table 2 contains the results of the parameter vector estimates for data sets 1-4 along withsummary statistics comparing to the Merton Garman parameter vector. Table 3 containsresults of the IVRMSE and standard deviations for each data set. The results of the datasets are described below: • Data set 1 : from Table 2 Panel 1 parameters κ, ξ appear to be challenging to estimate,being significantly larger and with quite high standard deviations, while α appearsto be stable. While Θ Pert. does not contain θ a significant amount of this missinginformation is absorbed by σ and ξ . Table 3 demonstrates the leading perturbativesolution does very well in approximating the Merton Garman model with an IVRMSEof 1 . . • Data set 2 : from Table 2 Panel 2 it starts to become clear some of the informationcontained in the correlation coefficient is absorbed by both ξ, α , particularly the latter.With the value of α reducing significantly while the standard deviation approximatelydoubles (relative to data set 1). Although it does appear that this regime is slightlyeasier to estimate κ, ξ . Table 3 reports a significant decrease (relative to data set 1) inIVRMSE with a moderate decrease in standard deviation. • Data set 3 : Table 2 Panel 3 demonstrates the absence of the information containedin the correlation coefficient has an effect on the ability to estimate κ, ξ, α in a similar10anor to that of Panel 1, observing very similar biases and standard deviations. Per-haps suggesting that it is more the non-zero nature of the correlation coefficient whichthe leading perturbative solution struggles with. Furthermore, Table 3 demonstratesvery similar errors to data set 1. • Data set 4 : from Table 2 Panel 4 the increase in central tendency also clearly has animpact in ξ, α the two variance related parameters of the leading perturbative solution.This is due to the increased difference in magnitude between the central tendency andinitial variance. It suggests that both parameters absorb missing information containedin θ . This difference manifests itself in Table 3 with a resulting IVRMSE of 1 . κ could certainly be a challenge to estimate, althoughthis is not unique to our approach. For substantial difference between variance and centraltendency it appears that α, ξ could also be a challenge, other than this estimating α is gener-ally inconsiderable. Regarding the matter of information absorption we note that it is clearthat σ absorbs significant amount of the information contained in θ , with contribution from α for large disparity between initial variance and central tendency. The parameters ξ, α alsoseem to share the majority of the information contained in ρ . Table 3 indicates that acrossdata sets the leading perturbative solution does well in approximating the Merton Garmanmodel with fairly consistent errors, given the standard deviations, with an approximate errorrange of 1 . − . In this section, we shall first discuss Galilean invariance, see e.g. (Bose (1995) and L´evy-Leblond (1967)), in the context of mathematical finance before explaining how new optionpricing models can be constructed using the concept of symmetries. We shall assume thatthe dimensionless option price ψ ( x, y, t ), from which the usual price of the option is derived,is the fundamental quantity. If we posit that ψ ( x, y, t ) is a measurable quantity, it shouldnot depend on the coordinate system P that is being used to measure it. We could use P parametrized by ( x, y, t ) or P (cid:48) parametrized by ( x (cid:48) , y (cid:48) , t (cid:48) ) and obtain the same dimensionlessprice, assuming that the two coordinate systems are related by a transformation whichwe shall take to be a Galilean transformation, knowing that ψ ( x, y, t ) is a solution to the2+1 heat equation. A Galilean transformation can be decomposed as the composition of arotation, a translation and a uniform motion in the space ( x, y, t ) where (cid:126)x = ( x, y ) represents11 point in two-dimensional space, and t a point in one-dimensional time. A general pointin this space can be ordered by a pair ( (cid:126)x, t ). Let us first consider the case where the twocoordinate systems are moving away from each other in a uniform motion fixed by a two-dimensional constant vector (cid:126)w . A uniform motion with vector (cid:126)w , is given by( (cid:126)x, t ) (cid:55)→ ( (cid:126)x + t (cid:126)w, t ) . (24)A translation is given by ( (cid:126)x, t ) (cid:55)→ ( (cid:126)x + (cid:126)a, t + s ) , (25)where (cid:126)a is two-dimensional vector and s a real number. Finally a rotation is given by( (cid:126)x, t ) (cid:55)→ ( G(cid:126)x, t ) , (26)where G is 2 × ψ ( x, y, τ ) (cid:18) ∂∂τ − ∂ ∂x − ∂ ∂y (cid:19) ψ ( x, y, τ ) = 0 , (27)is covariant under two-dimensional Galilean transformations Gal(2).Our Monte Carlo investigation of the Merton Garman model demonstrates that thesymmetry Gal(2) is a good symmetry of the model as the leading perturbation theory aroundthe symmetrical solution works very well. This demonstrates the importance of the symmetrybetween S and V which is manifest when expressed in the appropriate variables. To a verygood approximation and for a range of parameters relevant to financial applications, theMerton Garman model possesses a hidden Gal(2) symmetry that is only softly broken. Thisis very clearly illustrated by the results obtained in Figure (2) which demonstrates thatthe larger the volatility fluctuations are, the more the symmetry breaking terms becomeimportant. In panel 4 where the volatility is very close to being constant, the leadingorder correction is essentially vanishing and the symmetrical solution is very close to thatobtained with the original Merton Garman model. This is suggestive of an alternative wayfor building option pricing models. Instead of starting from stochastic processes, we couldsimply have derived the symmetrical model by positing that the option price should dependon the price of the underlying, a time dependent volatility and time. By requiring that thedimensionless option price follows a differential equation that is Gal(2) covariant. We wouldimmediately have obtained the 2+1 dimensional heat equation. The very small deviationsfrom the Gal(2) covariance, can be accounted for by symmetry breaking terms as explainedabove. This led us to an approximative perturbative analytical solution of the originalMerton Garman differential equation. 12here is another interesting consequence of having a symmetry group as a fundamentalbuilding block. Namely, one can classify how the different objects in the model will transformunder this symmetry. This is a very well developed field of mathematics called representa-tion theory. In the case above, ψ ( x, y, t ) is a scalar under Galilean transformations. A scalarrepresentation of a symmetry group means that it is invariant under the group transforma-tions. The differential equation is covariant under such transformations. Besides the scalartransformations, there are vector representations such as e.g. ∂ψ ( x, y, t ) /∂x , ∂ψ ( x, y, t ) /∂y or ∂ψ ( x, y, t ) /∂t which are nothing but the Greeks. They appear from a very differentperspective than is usually the case in finance.These ideas open up new directions in model building for option pricing. One could forexample consider a 3+1 dimensional model ( x, y, z, t ) where the z -direction could describea time dependent interest rate or additionally this could be used to model a two-factorvolatility process with stochastic central tendency, see e.g. Bardgett et al. (2018). Thedifferential equation for the price would be of the type ∂ψ ( x, y, t ) ∂t = ∂ ψ ( x, y, z, t ) ∂x + ∂ ψ ( x, y, z, t ) ∂y + ∂ ψ ( x, y, z, t ) ∂z , (28)for which it is easy to find solutions: G ( x, y, τ ; x (cid:48) , y (cid:48) , z (cid:48) , t ) = Θ( τ − t ) 1(4 π ( τ − t )) / exp (cid:18) − ( x − x (cid:48) ) + ( y − y (cid:48) ) + ( z − z (cid:48) ) τ − t ) (cid:19) . (29)One could also consider “relativistic” extensions of the Merton Garman model treating timeon the same footing as the underlying price and the volatility: ∂ ψ ( x, y, t ) ∂t = ∂ ψ ( x, y, t ) ∂x + ∂ ψ ( x, y, t ) ∂y . (30)It is well known that this model is covariant under Lorentz transformations. Clearly iden-tifying the right symmetry group for a given financial system is of paramount importance.Making use of symmetries to model physical system has been extremely successful in allfields of physics. Applying these ideas to option pricing models opens up new perspec-tives for model building in finance using the concept of symmetry groups and representationtheory. In this paper we have introduced a perturbative method to obtain analytical approximativesolutions to models such as the Merton Garman model. The key idea consists in treating the13rice of the underlying and the volatility in a symmetrical way. This leads to a model whichhas an exact Galilean invariance in two-dimensions as it is described by the two-dimensionalheat equation which has an analytical solution. By folding this solution with the boundarycondition leading to the correct price at maturity for a call option, we obtained an analyticalsymmetrical solution for this model which corresponds to the Black Scholes solution, despitebeing derived from a very different perspective and in a framework with a time dependentvolatility.The Merton Garman model is recovered by introducing symmetry breaking terms andwe have calculated the leading order correction to the symmetrical solution. We have shownthat our perturbative solution works very well by comparing it to a Monte Carlo simulationof the Merton Garman model for a range of parameters which are relevant from a financialpoint of view. The moneyness curves of the two prices are so overlapping that we had toplot implied volatility curves to be able to discuss in a quantitative manner the differencesbetween the two solutions.We argued that the fact that the symmetrical model works so well is a sign that theMerton Garman model has a hidden two-dimensional Galilean symmetry which is softlybroken for the relevant parameter ranges. We have explained that the concept of symmetry,groups and representation theory could be extremely useful in building pricing models infinancial mathematics. This clearly needs to be explored further. From this point of view,our work is opening up a new perspective on model building in mathematical finance.
Data Availability Statement:
Data sharing is not applicable to this article as no new datawere created or analyzed in this study.
A Heat equation
In this Appendix, we show how to massage the partial differential equation correspondingto the symmetrical model: ∂C∂t + rS ∂C∂S + 12 σ S ∂ C∂S + µV ∂C∂V + ξ V ∂ C∂V = rC, (31)into the heat equation in 2+1 dimensions. To do so, we now make the standard substitutionsfor the underlying and variance, transforming them to dimensionless variables: x = log( S/K ) , (32) y = log( V /V ) , (33)14here K is the strike price and V is some constant with units of 1 / sec. Re-casting Eq. (31)in terms of x and y , we obtain ∂C∂t + ω ∂C∂x + γ ∂C∂y + 12 σ ∂ C∂x + ξ ∂ C∂y = rC, (34)where, ω = r − σ and γ = µ − ξ .In order to remove the constant term, σ in front of the second derivative with respectto x , we make the time transformation: τ = σ ( T − t ), yielding: − ∂C∂τ + ( R − ∂C∂x + 2 γσ ∂C∂y + ∂ C∂x + η ∂ C∂y − R C = 0 , (35)where R = rσ and η = ξ σ . We then proceed with the further substitution: y = √ η y ,transforming the coefficient of the second derivative with respect to y to unity: − ∂C∂τ + ( R − ∂C∂x + √ γσξ ∂C∂y + ∂ C∂x + ∂ C∂y − R C = 0 . (36)Next, we make the price transformation: C ( x, y, τ ) = Kφ ( x, y, τ ) ψ ( x, y, τ ) , (37)with: φ ( x, y, τ ) = e ax + by + cτ . (38)The constants a, b and c are chosen by inspection, after substitution into Eq. (36) we seethat the choice: a = −
12 ( R − , (39) b = −
12 ( R − , (40) c = − (cid:32) ( R − + ( R + 1) (cid:33) , (41)where: R = 1 + √ γ/σξ , leads to the heat equation in 2+1 dimensions: ∂ψ ∂τ = ∂ ψ ∂x + ∂ ψ ∂y . (42)There are well known techniques to solve this partial differential equation analytically. Itis also well known that the heat equation has a symmetry based on the Galilean group in2+1 dimensions. This symmetry is now manifest. In particular, the variable x and y are15nterchangeable as advertised previously. Before we can solve the heat equation, we need tospecify an appropriate boundary equation.Guided by our intuition formed by the Black Scholes model, we propose a boundarycondition of the form for a call option: ψ ( x, y,
0) = e ( R − y (cid:32) e ( R +1) x/ − e ( R − x/ (cid:33) + . (43)To verify that this boundary condition makes sense from a option pricing point of view, wework out the implied boundary conditions for the call price in the original variables: C ( x, y,
0) = Ke ax + by ψ ( x, y, , (44)and thus C ( x, y,
0) = K (cid:32) e x − (cid:33) + . (45)Substitution back to original variable, using: S = Ke x , we find C ( S, V, T ) = (cid:32) S ( T ) − K (cid:33) + . (46)This is the standard payoff of a call option. We have thus found an appropriate boundarycondition for the heat equation and can thus proceed to solving this partial differentialequation. For a put option we have ψ ( x, y,
0) = − e ( R − y (cid:32) e ( R +1) x/ − e ( R − x/ (cid:33) + , (47)which leads to P ( S, V, T ) = (cid:32) K − S ( T ) (cid:33) + . (48) B Solution of the symmetrical model
In this Appendix, we show how to solve the heat equation in 2+1 dimensions: ∂ψ ∂τ = ∂ ψ ∂x + ∂ ψ ∂y . (49)So far the function ψ ( x, y, τ ) is only defined for τ >
0, however by introducing the Heavisidefunction Θ( τ ) we may extend the definition domain to the range τ < ψ ( x, y, τ ) = Θ( τ ) ψ ( x, y, τ ) , (50)16hus we get an inhomogeneous differential equation: (cid:32) ∂∂τ − ∂ ∂x − ∂ ∂y (cid:33) ¯ ψ ( x, y, τ ) = ¯ ψ ( x, y, δ ( τ ) . (51)This equation is solved by¯ ψ ( x, y, τ ) = (cid:90) ¯ ψ ( x (cid:48) , y (cid:48) , G ( x, y, τ | x (cid:48) , y (cid:48) , dx (cid:48) dy (cid:48) . (52)This is the partial differential equation for the Gaussian propagator of the heat equation in2+1 dimensions: G ( x, y, τ | X, Y,
0) = 14 πτ e − ( x − X )24 τ − ( y − Y )24 τ Θ( τ ) . (53)Combining the above two results, the solution can be written in the form¯ ψ ( x, y, τ ) = 14 πτ (cid:90) ¯ ψ ( X, Y, e − ( x − X )24 τ − ( y − Y )24 τ dXdY, (54)¯ ψ ( x, y, τ ) = 14 πτ (cid:90) e ( R − Y/ (cid:32) e ( R +1) X/ − e ( R − X/ (cid:33) + e − ( x − X )24 τ − ( y − Y )24 τ dXdY, (55)which leads to ψ ( x, y, τ ) = 14 πτ (cid:90) ψ ( X, Y, e − ( x − X )24 τ − ( y − Y )24 τ dXdY, (56)and ψ ( x, y, τ ) = 14 πτ (cid:90) e ( R − Y/ (cid:32) e ( R +1) X/ − e ( R − X/ (cid:33) + e − ( x − X )24 τ − ( y − Y )24 τ dXdY. (57)Note that the two integrals can be separated: ψ ( x, y, τ ) = 1 √ πτ (cid:90) ∞ dX (cid:32) e ( R +1) X/ − e ( R − X/ (cid:33) + e − ( x − X )24 τ × √ πτ (cid:90) ∞−∞ dY e ( R − Y/ e − ( y − Y )24 τ , (58)where the first of these integrals precisely corresponds to the 1+1 dimensional Black Scholesmodel. We thus finally obtain ψ ( x, y, τ ) = e ( R − ( τ ( R − y ) (cid:34) e ( R +1) x/ R +1) τ/ N ( d ) − e ( R − x/ R − τ/ N ( d ) (cid:35) , (59)17here N ( d ) = 1 √ π (cid:90) d −∞ exp (cid:18) − z (cid:19) dz, (60)and d = x √ τ + √ τ R + 1) = log( S/K ) + ( r + σ / T − t ) σ √ T − t , (61) d = x √ τ + √ τ R −
1) = d − σ √ T − t. (62)We have obtained an analytical solution to the symmetrical model. Remarkably, because ofthe boundary condition that only depends on S , it is identical to the Black Scholes solution.Going back to the original variables we find: C ( S, V, t ) = S N ( d ) − Ke − r ( T − t ) N ( d ) . (63) C Symmetry Breaking terms and solution to the Mer-ton Garman model
In this Appendix, we give details of the derivation of the perturbative solution to the fullMerton Garman. We first need to restore the full model by re-introducing the symmetrybreaking terms ∂C∂t + rS ∂C∂S + 12 σ S ∂ C∂S + c S (cid:32) V − σ (cid:33) ∂ C∂S + µV ∂C∂V + c λ ∂C∂V + ξ V ∂ C∂V + c (cid:32) ξ V α − ξ V (cid:33) ∂ C∂V + c ρξV α +1 / S ∂ C∂S∂V − rC = 0 . (64)Note that we have introduced dimensionless coefficients c i which denote the strength of thesymmetry breaking terms. In the limit c i = 1 one recovers the original Merton Garmanmodel. These coefficients are simply introduced as a bookkeeping trick to keep track ofwhich terms correspond to a deviation of the 2+1 Galilean invariant theory. In the end ofthe day, we set c i = 1. We can now apply the same variables transformations to Eq. (64)that we had applied in the symmetric case and obtain (cid:18) ∂∂τ − ∂ ∂x − ∂ ∂y + D ( x, y ) (cid:19) ψ ( x, y, τ ) = 0 , (65)18here the operator D ( x, y ) is defined as: D ( x, y ) = (cid:18) c (cid:0) V e y − σ (cid:1) (cid:18) ( a − a ) + (2 a − ∂∂x + ∂ ∂x (cid:19) (66)+ c λV e − y (cid:18) ∂∂y + b (cid:19) + c (cid:0) ξ V α − e (2 α − y − ξ (cid:1) (cid:18) ( b − b ) + (2 b − ∂∂y + ∂ ∂y (cid:19) + c ξρV α − e ( α − ) y (cid:18) ab + b ∂∂x + a ∂∂y + ∂ ∂x∂y (cid:19)(cid:19) , where a and b are respectively given in Eq. (39) and Eq. (40). Note that D ( x, y ) is τ independent.We now do perturbation theory around the symmetrical solution ψ which was given inEq. (59). To leading order in c i , we write ψ = ψ + ψ where ψ is of order c i . Keeping inmind that D is order c i , we find (cid:18) ∂∂τ − ∂ ∂x − ∂ ∂y (cid:19) ψ ( x, y, τ ) = −D ( x, y ) ψ ( x, y, τ ) . (67)This equation can be solved by the Green’s function method, we obtain ψ ( x, y, τ ) = − (cid:90) τ dt (cid:90) ∞−∞ dx (cid:48) (cid:90) ∞−∞ dy (cid:48) G ( x, y, τ ; x (cid:48) , y (cid:48) , t ) D ( x (cid:48) , y (cid:48) ) ψ ( x (cid:48) , y (cid:48) , t ) , (68)where G ( x, y, τ ; x (cid:48) , y (cid:48) , t ) = 14 π ( τ − t ) exp (cid:18) − ( x − x (cid:48) ) + ( y − y (cid:48) ) τ − t ) (cid:19) . (69)These integrals can be performed analytically. Our result provides an approximative andanalytical solution to the Merton Garman model. We find: ψ ( x, y, τ ) = ψ ( x, y, τ ) + ψ ( x, y, τ ) , (70)with ψ ( x, y, τ ) = c (cid:18) σ τ (cid:0) √ γ + σξ (cid:1) − σξ V e y (cid:18) e √ γτσξ + τ − (cid:19)(cid:19) e (cid:18) γ τ σξ − x τ + γy √ σξ (cid:19) √ πτ (cid:0) √ γ + σξ (cid:1) , (71)to leading order. In the original variables, we find: C ( S, V, t ) = − K (cid:0) SK (cid:1) − rσ e ( SK ) + ( r + σ ) t − T )28 σ t − T ) √ π (cid:16) √ γσξ + 1 (cid:17) (cid:112) σ ( T − t ) (72) × (cid:32) σ (cid:32) √ γσξ + 1 (cid:33) ( t − T ) + V (cid:18) e σ (cid:16) √ γσξ +1 (cid:17) ( T − t ) − (cid:19)(cid:33) , c = 1.It may appear surprising that the leading order correction does not depend on the sym-metry breaking terms parametrized by c , c and c . To understand what is happening,one can organize the perturbation theory slightly different, but mathematically equivalent,fashion by looking at corrections to the Green’s function G ( x, y, τ | x (cid:48) , y (cid:48) , (cid:18) ∂∂τ − ∂ ∂x − ∂ ∂y + D (cid:19) G ( x, y, τ | x (cid:48) , y (cid:48) ,
0) = δ ( τ ) δ ( x − x (cid:48) ) δ ( y − y (cid:48) ) . (73)Perturbation theory is organized by expanding around the Green’s function of the symmet-rical symmetry: G ( x, y, τ | x (cid:48) , y (cid:48) ,
0) = G ( x, y, τ | x (cid:48) , y (cid:48) ,
0) + G ( x, y, τ | x (cid:48) , y (cid:48) ,
0) to leading order.One obtains: (cid:18) ∂∂τ − ∂ ∂x − ∂ ∂y (cid:19) G ( x, y, τ | x (cid:48) , y (cid:48) ,
0) = δ ( τ ) δ ( x − x (cid:48) ) δ ( y − y (cid:48) ) . (74)The solution to this partial differential equation was given above. The correction G ( x, y, τ | x (cid:48) , y (cid:48) , (cid:18) ∂∂τ − ∂ ∂x − ∂ ∂y (cid:19) G ( x, y, τ | x (cid:48) , y (cid:48) ,
0) = −D ( x, y ) G ( x, y, τ | x (cid:48) , y (cid:48) , , (75)which can be solved easily. One finds: G ( x, y, τ | x (cid:48) , y (cid:48) ,
0) = − (cid:90) τ dt (cid:48) (cid:90) ∞−∞ dx (cid:48)(cid:48) (cid:90) ∞−∞ dy (cid:48)(cid:48) G ( x, y, τ | x (cid:48)(cid:48) , y (cid:48)(cid:48) , t (cid:48) ) D ( x (cid:48)(cid:48) , y (cid:48)(cid:48) ) G ( x (cid:48)(cid:48) , y (cid:48)(cid:48) , t (cid:48) | x (cid:48) , y (cid:48) , . (76)It is easy to show that this correction depends on all four symmetry breaking terms. Wefind G ,c ( x, y, τ | x (cid:48) , y (cid:48) ,
0) = c πτ / e − ( x − x (cid:48) )2+( y − y (cid:48) )24 τ (cid:0) ( − R ) τ + 2 τ ( − R ( x − x (cid:48) )) + ( x − x (cid:48) ) (cid:1) × (cid:18) σ √ τ − √ πe τ y − y (cid:48) )2+2 τ ( y − y (cid:48) )4 τ V (cid:18) Erf (cid:18) τ + y − y (cid:48) √ τ (cid:19) + Erf (cid:18) τ − y + y (cid:48) √ τ (cid:19)(cid:19)(cid:19) , (77) G ,c ( x, y, τ | x (cid:48) , y (cid:48) ,
0) = c λ πV τ e − ( x − x (cid:48) )2+ y y (cid:48) τ ( y − y (cid:48) )4 τ (cid:32) e yy (cid:48) τ ( e y − e y (cid:48) ) τ + e ( τ + y ) +2 τy (cid:48) + y (cid:48) τ √ π ( R − τ / (cid:18) Erf (cid:18) τ − y + y (cid:48) √ τ (cid:19) − Erf (cid:18) − τ − y + y (cid:48) √ τ (cid:19)(cid:19) (cid:33) , (78)20 ,c ( x, y, τ | x (cid:48) , y (cid:48) ,
0) = c τ π (1 − α ) e − ( x − x (cid:48) )24 τ (79) × (cid:32) e − ( y − y (cid:48) )24 τ − x − x (cid:48) ) ( − α − ξ e y + y (cid:48) ) (cid:0)(cid:0) R − (cid:1) τ + 2 τ ( R y − R y (cid:48) −
1) + ( y − y (cid:48) ) (cid:1) + ξ V α − (cid:16) − e αy + y (cid:48) ) (cid:17) ( − τ ( α + R − − y + y (cid:48) ) − ξ V α − e y + αy (cid:48) ) (2 τ (3 α + R −
3) + y − y (cid:48) ))+ √ πξ τ / V α − (2 α + R − α + R − e ( α − α − τ + y + y (cid:48) ) (cid:18) Erf (cid:18) α − τ + y − y (cid:48) √ τ (cid:19) − Erf (cid:18) − α − τ + y − y (cid:48) √ τ (cid:19)(cid:19) (cid:33) , and G ,c ( x, y, τ | x (cid:48) , y (cid:48) ,
0) = c ξρV α − (( R − τ + x − x (cid:48) )32 πτ (1 − α ) e − ατ y ( τ + y )+( x − x (cid:48) )2+ y (cid:48) τ (80) (cid:32) (cid:16) e αy + y (cid:48) − e y + αy (cid:48) (cid:17) e (cid:16) ατ + y (cid:16) y (cid:48) τ − (cid:17) − y (cid:48) (cid:17) −√ πτ (2 α + 2 R − (cid:18) Erf (cid:18) − ατ + τ + 2 y − y (cid:48) √ τ (cid:19) − Erf (cid:18) ατ − τ + 2 y − y (cid:48) √ τ (cid:19)(cid:19) exp (cid:18) α τ + 8 ατ ( y + y (cid:48) ) + ( τ − y (cid:48) ) + 4 y τ (cid:19) (cid:33) . It is easy to see that when folding these corrections to the Green’s function with theboundary condition (43) that only the contribution from the c term survives and one recoversthe result of (71). The boundary condition thus implies that the contributions of c , c and c vanish to leading order in the perturbation theory. These symmetry breaking terms will,however, contribute to higher order corrections. Higher precision, if required, can be obtainedby going to higher order in perturbation theory. Option prices can be calculated extremelyrapidly using this formalism. Note that, in principle, if we resummed perturbation theoryto all order in c i , the dependence on σ and ξ would vanish. It is also worth noticing thatour results are independent on V which is only introduced to match the dimension of V .21 Time Series Simulation
The Merton Garman model is defined by the coupled two-dimensional SDE: dS t = rS t dt + (cid:112) V t S t dW St , (81) dV t = κ ( θ − V t ) dt + ξV αt dW Vt , (82)with the two Brownian motion components: W St , W Vt having correlation ρ . We use an Eulerdiscretization scheme for the asset price and variance process, with a full truncation to tacklethe issue of negative variances. Conditional on a time s for t > s the discretization schemefor the asset price and variance processes are: S t = S s + rS t ∆ t + S s (cid:113) V + t ∆ tz S , (83) V t = V s − κ ( θ − V + s )∆ t + ξ ( V + s ∆ t ) α z V , (84)where ∆ t = t − s , z V ∼ N (0 ,
1) and z S = ρz V + √ − ρz with z ∼ N (0 , t = 7 days. At each observation time we simulate six uniquematurities within [7 , K/S ∈ [0 . , .
1] acrossten strikes. Each option price is computed using the Monte Carlo framework with 50 , / th of a trading day.The calibration process is done using the objective function defined in Eq. (23). It shouldbe noted that as initial conditions we start with the true parameter vector, i.e Θ pert.initial = Θ MG true and for σ we start with the initial variance. We also pass the variance path at each timestep. 22igure 1: Comparative fit of four different initial volatilities. Odd numbered panels (1,3,5,7)display the natural logarithm difference of the prices, while even numbered panels (2,4,6,8)display the implied volatility curves. For the implied volatility panels the solid black linerepresents the Monte Carlo implied volatility and the dashed line is that of the leading orderorder perturbative solution. We pick an option maturity of 30 days.23igure 2: Test of the validity of the perturbation theory. These panels depict the ratios C / ( C + C ) in % where C is the contribution to the price of the symmetrical solutionand C the leading order correction in perturbation theory. Clearly C (cid:28) C even when thevolatility is large. This demonstrates the validity of the perturbative expansion. The fourcases correspond respectively to panels (1,3,5,7) of Figure 1.24able 1: Implied Volatility Root Mean Squared Error (IVRMSE, defined in Eq. (23)) whencalibrated the leading order order perturbative solution and compared to the Monte Carlosimulation of the Merton Garman model for four different initial volatilities for a time horizonof 30 days. Also reported is the calibrated values of σ for each scenario. V σ
IVRMSE0.3500 0.3254 0.00550.2500 0.2361 0.00370.1800 0.1730 0.00700.1000 0.1069 0.015725able 2: Simulated calibration exercise for the leading order order perturbative solution with100 paths, per data set. The MG row displays the true parameter vector. The leading orderperturbative solution (Pert.) row displays the mean calibrated results of the leading orderorder perturbative solution. The std. error row displays the stand deviation of the leadingorder order perturbative solution parameters.Panel 1. data set 1 κ θ ξ ρ α σ MG 1.1768 0.0823 0.3000 -0.5459 1.0000 -Pert. 4.1905 - 0.8772 - 0.9698 0.0855Bias 3.0137 - 0.5772 - 0.0302 -Std. Error 3.6179 - 0.4163 - 0.1717 0.0120Panel 2. Dataset 2 κ θ ξ ρ α σ MG 1.1768 0.0823 0.3000 0.0000 1.0000 -Pert. 2.5454 - 0.6326 - 0.6178 0.0810Bias 1.3686 - 0.3326 - 0.3822 -Std. Error 3.3240 - 0.4805 - 0.3584 0.0133Panel 3. Dataset 3 κ θ ξ ρ α σ MG 1.1768 0.0823 0.3000 0.5459 1.0000 -Pert. 4.4657 - 0.8806 - 0.9585 0.0851Bias 3.2889 - 0.5806 - 0.0415 -Std. Error 3.3320 - 0.3094 - 0.1448 0.0120Panel 4. Dataset 4 κ θ ξ ρ α σ MG 1.1768 0.1250 0.3000 -0.5459 1.0000 -Pert. 5.4417 - 1.0104 - 1.6696 0.1022Bias 4.2649 - 0.7104 - 0.6696 -Std. Error 4.0231 - 0.4529 - 1.4916 0.014526able 3: Displays IVRMSE between the perturbative solution (Pert.) and the solution ofthe Merton Garman model for each data set. IVRMSE is calculated using 23. The std. errorrow denotes the standard deviation of the IVRMSE.data set 1 data set 2 data set 3 data set 4IVRMSE 0.0130 0.0117 0.0132 0.0172Std. Error 0.0035 0.0031 0.0036 0.006127 eferences
Aguilar, J.-P.
A series representation for the black-scholes formula. arXiv preprintarXiv:1710.01141 (2017).
Aı, Y., Kimmel, R., et al.
Maximum likelihood estimation of stochastic volatilitymodels.
Journal of financial economics 83 , 2 (2007), 413–452.
Ajaib, M. A.
Introducing spin in 2d quantum tunneling. arXiv preprint arXiv:1705.10409 (2017).
Andersen, T. G., Benzoni, L., and Lund, J.
An empirical investigation of continuous-time equity return models.
The Journal of Finance 57 , 3 (2002), 1239–1284.
Baaquie, B. E.
A path integral approach to option pricing with stochastic volatility:some exact results.
Journal de Physique I 7 , 12 (1997), 1733–1753.
Baaquie, B. E., Coriano, C., and Srikant, M.
Quantum mechanics, path integralsand option pricing: Reducing the complexity of finance. In
Nonlinear Physics: Theory andExperiment II . World Scientific, 2003, pp. 333–339.
Bachelier, L.
Th´eorie de la sp´eculation . Gauthier-Villars, 1900.
Bakshi, G., Cao, C., and Chen, Z.
Empirical performance of alternative option pricingmodels.
The Journal of finance 52 , 5 (1997), 2003–2049.
Bardgett, C., Gourier, E., and Leippold, M.
Inferring volatility dynamics and riskpremia from the s&p 500 and vix markets.
Journal of Financial Economics (2018).
Bates, D. S.
Post-’87 crash fears in the s&p 500 futures option market.
Journal ofEconometrics 94 , 1-2 (2000), 181–238.
Bates, D. S.
Maximum likelihood estimation of latent affine processes.
The Review ofFinancial Studies 19 , 3 (2006), 909–965.
Benzoni, L.
Pricing options under stochastic volatility: an empirical investigation. Tech.rep., working paper, University of Minnesota, 2002.
Black, F., and Scholes, M.
The pricing of options and corporate liabilities.
Journalof political economy 81 , 3 (1973), 637–654.
Blazhyevskyi, L., and Yanishevsky, V.
The path integral representation kernel ofevolution operator in merton-garman model. arXiv preprint arXiv:1106.5143 (2011).28 ose, S.
Representations of the (2+ 1)-dimensional galilean group.
Journal of Mathemat-ical Physics 36 , 2 (1995), 875–890.
Chernov, M., Gallant, A. R., Ghysels, E., and Tauchen, G.
Alternative modelsfor stock price dynamics.
Journal of Econometrics 116 , 1-2 (2003), 225–257.
Chourdakis, K., and Dotsis, G.
Maximum likelihood estimation of non-affine volatilityprocesses.
Journal of Empirical Finance 18 , 3 (2011), 533–545.
Christoffersen, P., Feunou, B., Jacobs, K., and Meddahi, N.
The economicvalue of realized volatility: Using high-frequency returns for option valuation.
Journal ofFinancial and Quantitative Analysis 49 , 3 (2014), 663–697.
Christoffersen, P., Jacobs, K., and Mimouni, K.
Volatility dynamics for thes&p500: evidence from realized volatility, daily returns, and option prices.
The Review ofFinancial Studies 23 , 8 (2010), 3141–3189.
Duan, J.-C., and Yeh, C.-Y.
Jump and volatility risk premiums implied by vix.
Journalof Economic Dynamics and Control 34 , 11 (2010), 2232–2244.
Eraker, B.
Do stock prices and volatility jump? reconciling evidence from spot andoption prices.
The Journal of Finance 59 , 3 (2004), 1367–1403.
Eraker, B., Johannes, M., and Polson, N.
The impact of jumps in volatility andreturns.
The Journal of Finance 58 , 3 (2003), 1269–1300.
Garman, M. B., et al.
A general theory of asset valuation under diffusion state pro-cesses. Tech. rep., University of California at Berkeley, 1976.
Heston, S. L.
A closed-form solution for options with stochastic volatility with applica-tions to bond and currency options.
The review of financial studies 6 , 2 (1993), 327–343.
Huang, J.-z., and Wu, L.
Specification analysis of option pricing models based ontime-changed l´evy processes.
The Journal of Finance 59 , 3 (2004), 1405–1439.
Hull, J., and White, A.
The pricing of options on assets with stochastic volatilities.
The journal of finance 42 , 2 (1987), 281–300.
Hull, J. C.
Options, futures and other derivatives. -prentice hall, upper saddle river, newjersey.
Jones, C. S.
The dynamics of stochastic volatility: evidence from underlying and optionsmarkets.
Journal of econometrics 116 , 1-2 (2003), 181–224.29 aeck, A., and Alexander, C.
Volatility dynamics for the s&p 500: Further evidencefrom non-affine, multi-factor jump diffusions.
Journal of Banking & Finance 36 , 11 (2012),3110–3121.
Kleinert, H., and Korbel, J.
Option pricing beyond black–scholes based on double-fractional diffusion.
Physica A: Statistical Mechanics and its Applications 449 (2016),200–214.
L´evy-Leblond, J.-M.
Nonrelativistic particles and wave equations.
Communications inMathematical Physics 6 , 4 (1967), 286–311.
Merton, R. C.
Theory of rational option pricing.
The Bell Journal of economics andmanagement science (1973), 141–183.
Pan, J.
The jump-risk premia implicit in options: Evidence from an integrated time-seriesstudy.
Journal of financial economics 63 , 1 (2002), 3–50.
Park, Y.-H.
Price dislocation and price discovery in the s&p 500 options and vix deriva-tives markets. Tech. rep., Working Paper). Federal Reserve Board, 2015.
Utama, B., and Purqon, A.
Feynman path integral application on deriving black-scholes diffusion equation for european option pricing. In