An analytical perturbative solution to the Merton Garman model using symmetries

In this paper, we introduce an analytical perturbative solution to the Merton Garman model. It is obtained by doing perturbation theory around the exact analytical solution of a model which possesses a two-dimensional Galilean symmetry. We compare our perturbative solution of the Merton Garman model to Monte Carlo simulations and find that our solutions performs surprisingly well for a wide range of parameters. We also show how to use symmetries to build option pricing models. Our results demonstrate that the concept of symmetry is important in mathematical finance.

1 Introduction calculations of option prices. This requires us to identify a "symmetric" version of the model which can easily be solved analytically. One can then reintroduce the symmetry breaking terms of the original Merton Garman differential equation and do perturbation theory around the symmetric solution thereby obtaining an approximative but analytical solution to the original Merton Garman differential equation. We then propose a new approach to model building in option pricing based on the concept of symmetry groups and representation theory. This concept has been extremely successful in modern physics. It is at the origin of 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 very specific 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 equation which describes the Merton Garman model. In section 3, we explain how to reduce the original Merton Garman to a simple, symmetrical, model. We present an exact analytical solution to the symmetrical model. We then restore the original Merton Garman by reintroducing the symmetry breaking terms and provide an analytical perturbative solution to the Merton Garman model. In section 4, we compare our solutions to different numerical solutions found in the literature. In section 5, we propose a new approach to model building in mathematical finance. Finally, we conclude in section 6.

The Merton Garman model
In the Merton Garman model, the price of an option is dependent on the time t, the price of the underlying S and the volatility V . Both S and V are taken to be time-dependent functions and thus the Merton Garman model has the potential to provide a more accurate calculation 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 and for the volatility V which resembles a stochastic, mean reverting, volatility regime. Here, ξ is the standard deviation 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, 1/2 recovers the Hull and White (see Hull and White (1987)) and Heston models, respectively. However, we do not constrain ourselves to either of these worlds. Here, α 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 a straightforward manner. Our first step is to find the associated partial differential equation which 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 in Eqs (1) and (2) is given by: where C is the price of a call.
Using this formula, one obtains ∂C ∂t + rS ∂C ∂S + 1 2 V S 2 ∂ 2 C ∂S 2 + (λ + µV ) where λ = κθ, µ = −κ. The call price C = C(S, V, t) depends on the time t, the price of the underlying S and the time dependent volatility V = V (t). In this model, there are three free parameters λ, µ and α. As we explained previously, existing solutions to this partial differential equation are numerical ones which have been obtained using Monte Carlo methods. Note that the put price P = P (S, V, t) fulfills the same differential equation, but it is obviously subject to a different boundary condition.
3 Reduction to the symmetrical model and perturbative solution to the Merton Garman model By studying the partial differential equation given in Eq. (4), it quickly becomes clear that the difficulty in finding an analytical solution to Eq. (4) is due to the lack of symmetry between the different terms of the partial differential equation. It is useful to study the dimensions of the different terms and constants in this partial differential equation. The price of the call is obviously given in a specific currency which we shall take to be the USD or $. The remaining dimensions follow from this. We have: It is instructive to see that S and V have different dimensions. Nevertheless, our goal is to treat S and V as symmetrically as possible to make a global Galilean invariance in 2+1 manifest (see section 5). This can be achieved by adequate variable transformations and by identifying the terms in the differential equation that violate this symmetry.

Symmetrical model
Our aim is to derive a differential equation that is symmetrical in S and V . With this aim in mind, let us introduce an averaged volatility σ 2 which is constant. As in the case of the Black Scholes model, different definitions for the averaged volatility are possible, the specific choice will not impact our methodology and results.
Inspecting the differential equation (4), it is clear that we need to pick α = 1 to emphasize the symmetry between S and V . We thus consider We need to keep in mind that we will need to reintroduce 1 2 V S 2 ∂ 2 C ∂S 2 , λ ∂C ∂V and the terms corresponding to deviations from 1 for α. Note that ξ 0 is different from ξ, in particular they do not have the same dimensions. Finally, we see that there is a mixed derivative term which needs to be eliminated. We thus set ρ = 0 and we will reintroduce this term as symmetry breaking term. We thus end up with: This partial differential equation can be massaged with standard substitutions into a 2+1 dimensional heat equation (see Appendix A) in which case the symmetry in S and V becomes manifest. In order to do so, we introduce and C(x, y, τ ) = Kφ(x, y, τ )ψ 0 (x, y, τ ), where K is the strike price and V 0 is some constant with units of 1/sec. The function φ(x, y, τ ) and the rescaled time τ are defined in Appendix A. Standard manipulations described in Appendix A lead to which is manifestly symmetrical in x and y. We will thus refer to the model described by the differential equation (31) as the symmetrical model. Another reason for massaging the symmetrical model into a heat equation is that this equation is easy to solve analytically.
We impose the standard boundary condition for the call price: For a put option we have

Solution of the symmetrical model
Details of the derivation of the analytical solution of the symmetrical model, i.e., of the 2+1 dimensional heat equations, are given in Appendix B. We find where and Remarkably, because of the boundary condition that only depends on S, it is identical to the Black Scholes solution.
For a put option, the very same procedure leads to In the next subsection, we shall restore the symmetry breaking terms and discuss the full Merton Garman model.

Symmetry Breaking terms and solution to the Merton Garman model
We are now in a position to solve the full Merton Garman model using perturbation theory around the symmetrical solution C 0 (S, V, t). We organize perturbation theory as an expansion in terms the coefficients of the symmetry breaking terms. We first need to restore the full model by re-introducing the symmetry breaking terms Note that we have introduced dimensionless coefficients c i which denote the strength of the symmetry breaking terms. In the limit c i = 1 one recovers the original Merton Garman model. These coefficients are simply introduced as a bookkeeping trick to keep track of which terms correspond to a deviation of the 2+1 Galilean invariant theory. In the end of the day, we set c i = 1. We now do perturbation theory around the symmetrical solution C 0 (S, V, t) and obtain where we have set c 1 = 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 σ 2 is not too large. In the limit when V is large, σ 2 is large as well and we expect that, as in the Black Scholes case, the price of the call becomes the price of the underlying S. Details of the derivation can be found in Appendix C.
It may appear surprising that the leading order correction does not depend on the symmetry breaking terms parametrized by c 2 , c 3 and c 4 . It can easily be shown (see Appendix C) that the boundary condition (43) insures that only the contribution from the c 1 term survives. The boundary condition implies that the contributions of c 2 , c 3 and c 4 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 obtained by going to higher order in perturbation theory. Option prices can be calculated extremely rapidly using this formalism. Note that, in principle, if we resummed perturbation theory to all order in c i , the dependence on σ and ξ 0 would vanish. It is also worth noticing that our results are independent on V 0 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 The prices obtained to leading order in perturbation theory for a call and put option are thus given by and In the next section we shall compare our results to exact solutions obtained numerically.

Comparison with numerical simulations
In this section we investigate how the approximative solution compares to a Monte Carlo simulation of the full Merton Garman model. It is well-known that the Merton Garman model is not solvable analytically, however it can be solved using numerical methods. The first step is a comparison of static cross sections of options. Then we compare using a simulated time series calibration exercise.

Static Cross-Section Comparison
The first step in evaluating the performance of the leading order perturbative solution is to compare to a multitude of simulated data of the Merton Garman model to ensure that the approximation is sufficient to fit a range of different options at one time. We start by describing the data simulation of the Merton Garman model, then move onto the calibration procedure for the approximative solution and discuss the results.
We choose a standard Monte Carlo framework, using stratified sampling and antithetic variables, simulating seven million paths with a time step of one-tenth of a day for option maturities. We choose a spot underlying price of S = $100 and a strike range of K ∈ [90, 110] to give a moneyness range of K/S ∈ [0.9, 1.1] to simulate call options throughout the spectrum of moneyness 3 . We simulate the Merton Garman model with the structural parameter vector: Θ M G = {1.5, 0.08, 1.5, −0.5, 1} 4 and initial volatility V (t = 0) ∈ [10, 35]%.
The leading order perturbative solution is independent on the the symmetry breaking terms, characterized by c 2 , c 3 and c 4 . The solution is thus independent of the parameters: ρ, θ. However, as a by-product of the perturbation theory we introduced the following parameters: ξ 0 , σ. From inspection we fix ξ 0 using ξ 0 = ξσ 2(α−1) which guarantees that it has the right dimensions. While σ remains to be determined by calibration. Yielding the parameter vector: Θ pert. = {κ, ξ, α, σ}.
The parameter σ is determined by calibration. When fitting σ to these simulations, there is a risk of overfitting expensive out of the money (OTM) options. For that reason, it is best to consider the implied volatility objective function (this is noted in Christoffersen et al. (2014)) which is given by: where IV M C i stands for the implied volatility of the i th -option simulated using the Monte Carlo and IV pert. i is the implied volatility of the i th -option calculated using the leading order perturbative solution. The calibration exercise is extremely fast as our formula for option prices is an approximative analytical solution. Figure (1) and Table 1 demonstrates the results of the static calibration exercise for a 30 day maturity horizon 5 . It should be noted that in Figure (1) the price Panels are the difference of the log, this is needed to 3 In this exercise we simulate call options only, as put prices can be calculated from the put-call parity. 4 We also simulate for α ∈ [0.75, 1.5] and κ ∈ [1.5, 5], ρ ∈ [−0.5, −0.9], θ ∈ [0.08, 0.15]. However, the results of the calibration exercise are represented by the choice of parameters made above. 5 While we simulated time horizons between 5-100 days, this is representative of our results and we drop the other results for brevity.
observe any difference in prices as the two methods produce prices which are very similar. However, from these Panels it is clear that smaller moneyness, i.e K/S < 1 see extremely small errors where log(C M C /C pert. ) ∼ 0. Also apparent is that at some moneyness level the leading order perturbative solution will over price options, the level of moneyness at which this occurs is inversely proportional to volatility. Lastly, the range of under to over pricing is also inversely proportional to volatility. However, as the prices from both methods are very similar it is more informative to look at the implied volatility cross-section in the even Panels of Figure (1) along with the IVRMSE column of Table 1. These show throughout the range of volatility regimes the leading order perturbative solution is able to approximate successfully a low IVRMSE, with a maximum occurring from the low volatility regime (Panel 8) of 1.57%.
Another way to test the consistency of the perturbative expansion is to consider the ratio C 1 /(C 0 + C 1 ) as a function of moneyness, K/S , where C 0 is the contribution to the price of the symmetrical solution and C 1 is the leading order correction in perturbation theory. Figure (2) shows these ratios for several cases. Clearly C 1 C 0 even when the volatility is large. This demonstrates nicely the validity of the perturbative expansion even in the case where the volatility is large. While this exercise confirms that for fixed scenarios the perturbative 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 be able to gauge the reliability of the perturbative approximation. This is what we shall focus on next.

Simulated Time Series Calibration
The second step in evaluating the performance of the leading order perturbative solution is to estimate it against a time series simulation of the Merton Garman model. The time series uses 100 different Monte Carlo paths to simulate the asset price and variance paths including 6 unique maturities within [7,180] days, for details see Appendix D.
The benefits of the stress test is two fold: firstly it is particularly pertinent to run a number of different simulations for different parameter values, specifically investigating the effect different θ, ρ has on performance, as the other parameters in the Θ pert. vector will have to attempt to absorb the information contained in the absent parameters. Secondly, it also provides a first estimate into the applicability of the perturbative solution to different types of options markets. We simulate four different data sets with varying parameter vectors, described below.
• data set 4 : Θ M G = {1.1768, 0.1250, 0.3000, −0.5459, 1.0000}, with a high(er) central tendency we investigate an equity style option market with a central tendency which is significantly higher than the initial variance value of: 0.08 and thus tests how the leading order perturbative solution handles significant change in the variance path.
For the following simulated calibration exercises, it is imperative to note the difference in IVRMSE and parameters between the data sets as this will highlight the following: firstly, parameter regions where the leading order perturbative solution might breakdown. Secondly, potential difficulties in estimating certain parameters of the model. Thirdly, where information 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 with summary statistics comparing to the Merton Garman parameter vector. Table 3 contains results of the IVRMSE and standard deviations for each data set. The results of the data sets 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 α appears to be stable. While Θ Pert. does not contain θ a significant amount of this missing information is absorbed by σ and ξ. Table 3 demonstrates the leading perturbative solution does very well in approximating the Merton Garman model with an IVRMSE of 1.2977% and standard deviation of 0.3474%.
• Data set 2 : from Table 2 Panel 2 it starts to become clear some of the information contained in the correlation coefficient is absorbed by both ξ, α, particularly the latter.
With the value of α reducing significantly while the standard deviation approximately doubles (relative to data set 1). Although it does appear that this regime is slightly easier to estimate κ, ξ. Table 3 reports a significant decrease (relative to data set 1) in IVRMSE with a moderate decrease in standard deviation.
• Data set 3 : Table 2 Panel 3 demonstrates the absence of the information contained in the correlation coefficient has an effect on the ability to estimate κ, ξ, α in a similar manor to that of Panel 1, observing very similar biases and standard deviations. Perhaps suggesting that it is more the non-zero nature of the correlation coefficient which the leading perturbative solution struggles with. Furthermore, Table 3 demonstrates very similar errors to data set 1.
• Data set 4 : from Table 2 Panel 4 the increase in central tendency also clearly has an impact in ξ, α the two variance related parameters of the leading perturbative solution. This is due to the increased difference in magnitude between the central tendency and initial variance. It suggests that both parameters absorb missing information contained in θ. This difference manifests itself in Table 3 with a resulting IVRMSE of 1.7199%, the largest across all data sets.
In summary, on the estimation side κ could certainly be a challenge to estimate, although this is not unique to our approach. For substantial difference between variance and central tendency it appears that α, ξ could also be a challenge, other than this estimating α is generally inconsiderable. Regarding the matter of information absorption we note that it is clear that σ absorbs significant amount of the information contained in θ, with contribution from α for large disparity between initial variance and central tendency. The parameters ξ, α also seem to share the majority of the information contained in ρ. Table 3 indicates that across data sets the leading perturbative solution does well in approximating the Merton Garman model with fairly consistent errors, given the standard deviations, with an approximate error range of 1.2 − 1.7%.

Model building in finance, symmetries and group theory
In this section, we shall first discuss Galilean invariance, see e.g. (Bose (1995) and Lévy-Leblond (1967)), in the context of mathematical finance before explaining how new option pricing models can be constructed using the concept of symmetries. We shall assume that the 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 should not depend on the coordinate system P that is being used to measure it. We could use P parametrized by (x, y, t) or P parametrized by (x , y , t ) and obtain the same dimensionless price, assuming that the two coordinate systems are related by a transformation which we shall take to be a Galilean transformation, knowing that ψ(x, y, t) is a solution to the 2+1 heat equation. A Galilean transformation can be decomposed as the composition of a rotation, a translation and a uniform motion in the space (x, y, t) where x = (x, y) represents a point in two-dimensional space, and t a point in one-dimensional time. A general point in this space can be ordered by a pair ( x, t). Let us first consider the case where the two coordinate systems are moving away from each other in a uniform motion fixed by a twodimensional constant vector w. A uniform motion with vector w, is given by A translation is given by where a is two-dimensional vector and s a real number. Finally a rotation is given by where G is 2 × 2 orthogonal transformation. It is easy to see that the partial differential equation obeyed by ψ(x, y, τ ) is covariant under two-dimensional Galilean transformations Gal(2).
Our Monte Carlo investigation of the Merton Garman model demonstrates that the symmetry Gal(2) is a good symmetry of the model as the leading perturbation theory around the symmetrical solution works very well. This demonstrates the importance of the symmetry between S and V which is manifest when expressed in the appropriate variables. To a very good approximation and for a range of parameters relevant to financial applications, the Merton Garman model possesses a hidden Gal(2) symmetry that is only softly broken. This is very clearly illustrated by the results obtained in Figure (2) which demonstrates that the larger the volatility fluctuations are, the more the symmetry breaking terms become important. In panel 4 where the volatility is very close to being constant, the leading order correction is essentially vanishing and the symmetrical solution is very close to that obtained with the original Merton Garman model. This is suggestive of an alternative way for building option pricing models. Instead of starting from stochastic processes, we could simply have derived the symmetrical model by positing that the option price should depend on the price of the underlying, a time dependent volatility and time. By requiring that the dimensionless option price follows a differential equation that is Gal(2) covariant. We would immediately have obtained the 2+1 dimensional heat equation. The very small deviations from the Gal(2) covariance, can be accounted for by symmetry breaking terms as explained above. This led us to an approximative perturbative analytical solution of the original Merton Garman differential equation.
There is another interesting consequence of having a symmetry group as a fundamental building block. Namely, one can classify how the different objects in the model will transform under this symmetry. This is a very well developed field of mathematics called representation theory. In the case above, ψ(x, y, t) is a scalar under Galilean transformations. A scalar representation of a symmetry group means that it is invariant under the group transformations. The differential equation is covariant under such transformations. Besides the scalar transformations, 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 different perspective than is usually the case in finance.
These ideas open up new directions in model building for option pricing. One could for example consider a 3+1 dimensional model (x, y, z, t) where the z-direction could describe a time dependent interest rate or additionally this could be used to model a two-factor volatility process with stochastic central tendency, see e.g. Bardgett et al. (2018). The differential equation for the price would be of the type for which it is easy to find solutions: (29) One could also consider "relativistic" extensions of the Merton Garman model treating time on the same footing as the underlying price and the volatility: It is well known that this model is covariant under Lorentz transformations. Clearly identifying 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 all fields of physics. Applying these ideas to option pricing models opens up new perspectives for model building in finance using the concept of symmetry groups and representation theory.

Conclusions
In this paper we have introduced a perturbative method to obtain analytical approximative solutions to models such as the Merton Garman model. The key idea consists in treating the price of the underlying and the volatility in a symmetrical way. This leads to a model which has an exact Galilean invariance in two-dimensions as it is described by the two-dimensional heat equation which has an analytical solution. By folding this solution with the boundary condition leading to the correct price at maturity for a call option, we obtained an analytical symmetrical solution for this model which corresponds to the Black Scholes solution, despite being derived from a very different perspective and in a framework with a time dependent volatility. The Merton Garman model is recovered by introducing symmetry breaking terms and we have calculated the leading order correction to the symmetrical solution. We have shown that our perturbative solution works very well by comparing it to a Monte Carlo simulation of the Merton Garman model for a range of parameters which are relevant from a financial point of view. The moneyness curves of the two prices are so overlapping that we had to plot implied volatility curves to be able to discuss in a quantitative manner the differences between the two solutions.
We argued that the fact that the symmetrical model works so well is a sign that the Merton Garman model has a hidden two-dimensional Galilean symmetry which is softly broken 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 in financial 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 data were created or analyzed in this study.

A Heat equation
In this Appendix, we show how to massage the partial differential equation corresponding to the symmetrical model: into the heat equation in 2+1 dimensions. To do so, we now make the standard substitutions for the underlying and variance, transforming them to dimensionless variables: where K is the strike price and V 0 is some constant with units of 1/sec. Re-casting Eq. (31) in terms of x and y, we obtain where, ω = r − 1 2 σ 2 and γ = µ − ξ 2 0 .
The constants a, b and c are chosen by inspection, after substitution into Eq. (36) we see that the choice: where: R 2 = 1 + √ 2γ/σξ 0 , leads to the heat equation in 2+1 dimensions: There are well known techniques to solve this partial differential equation analytically. It is also well known that the heat equation has a symmetry based on the Galilean group in 2+1 dimensions. This symmetry is now manifest. In particular, the variable x and y are interchangeable as advertised previously. Before we can solve the heat equation, we need to specify an appropriate boundary equation.
Guided by our intuition formed by the Black Scholes model, we propose a boundary condition of the form for a call option: To verify that this boundary condition makes sense from a option pricing point of view, we work out the implied boundary conditions for the call price in the original variables: C(x, y, 0) = Ke ax+by ψ 0 (x, y, 0), and thus Substitution back to original variable, using: S = Ke x , we find This is the standard payoff of a call option. We have thus found an appropriate boundary condition for the heat equation and can thus proceed to solving this partial differential equation. For a put option we have which leads to

B Solution of the symmetrical model
In this Appendix, we show how to solve the heat equation in 2+1 dimensions: So far the function ψ 0 (x, y, τ ) is only defined for τ > 0, however by introducing the Heaviside function Θ(τ ) we may extend the definition domain to the range τ < 0 Thus we get an inhomogeneous differential equation: This equation is solved bȳ This is the partial differential equation for the Gaussian propagator of the heat equation in 2+1 dimensions: Combining the above two results, the solution can be written in the form which leads to and ψ 0 (x, y, τ ) = 1 4πτ e (R 2 −1)Y /2 e (R 1 +1)X/2 − e (R 1 −1)X/2 dXdY. (57) Note that the two integrals can be separated: where the first of these integrals precisely corresponds to the 1+1 dimensional Black Scholes model. We thus finally obtain where and We have obtained an analytical solution to the symmetrical model. Remarkably, because of the boundary condition that only depends on S, it is identical to the Black Scholes solution.
Going back to the original variables we find:

C Symmetry Breaking terms and solution to the Merton Garman model
In this Appendix, we give details of the derivation of the perturbative solution to the full Merton Garman. We first need to restore the full model by re-introducing the symmetry breaking terms Note that we have introduced dimensionless coefficients c i which denote the strength of the symmetry breaking terms. In the limit c i = 1 one recovers the original Merton Garman model. These coefficients are simply introduced as a bookkeeping trick to keep track of which terms correspond to a deviation of the 2+1 Galilean invariant theory. In the end of the 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 where the operator D(x, y) is defined as: 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 ψ 0 which was given in Eq. (59). To leading order in c i , we write ψ = ψ 0 + ψ 1 where ψ 1 is of order c i . Keeping in mind that D is order c i , we find This equation can be solved by the Green's function method, we obtain where These integrals can be performed analytically. Our result provides an approximative and analytical solution to the Merton Garman model. We find: ψ(x, y, τ ) = ψ 0 (x, y, τ ) + ψ 1 (x, y, τ ), with to leading order. In the original variables, we find: where we have set c 1 = 1.
It may appear surprising that the leading order correction does not depend on the symmetry breaking terms parametrized by c 2 , c 3 and c 4 . 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 , y , 0). The differential equation to solve is given by: Perturbation theory is organized by expanding around the Green's function of the symmetrical symmetry: G(x, y, τ |x , y , 0) = G 0 (x, y, τ |x , y , 0) + G 1 (x, y, τ |x , y , 0) to leading order. One obtains: The solution to this partial differential equation was given above. The correction G 1 (x, y, τ |x , y , 0) is obtained by solving: which can be solved easily. One finds: dy G 0 (x, y, τ |x , y , t )D(x , y )G 0 (x , y , t |x , y , 0).
(76) It is easy to show that this correction depends on all four symmetry breaking terms. We find G 1,c 1 (x, y, τ |x , y , 0) = c 1 64πτ 5/2 e − (x−x ) 2 +(y−y ) 2 G 1,c 2 (x, y, τ |x , y , 0) = c 2 λ 16πV 0 τ 2 e − (x−x ) 2 +y 2 +y 2 +4τ (y−y ) 4τ 2e yy 2τ (e y − e y )τ + e (τ +y) 2 +2τ y +y 2 4τ √ π(R 2 − 2)τ 3/2 It is easy to see that when folding these corrections to the Green's function with the boundary condition (43) that only the contribution from the c 1 term survives and one recovers the result of (71). The boundary condition thus implies that the contributions of c 2 , c 3 and c 4 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 obtained by going to higher order in perturbation theory. Option prices can be calculated extremely rapidly using this formalism. Note that, in principle, if we resummed perturbation theory to all order in c i , the dependence on σ and ξ 0 would vanish. It is also worth noticing that our results are independent on V 0 which is only introduced to match the dimension of V .

D Time Series Simulation
The Merton Garman model is defined by the coupled two-dimensional SDE: with the two Brownian motion components: W S t , W V t having correlation ρ. We use an Euler discretization scheme for the asset price and variance process, with a full truncation to tackle the issue of negative variances. Conditional on a time s for t > s the discretization scheme for the asset price and variance processes are: where ∆t = t − s, z V ∼ N (0, 1) and z S = ρz V + √ 1 − ρz with z ∼ N (0, 1). This scheme is used to generate 100 different sample paths of weekly returns and latent variance over one year, i.e 52 observations with ∆t = 7 days. At each observation time we simulate six unique maturities within [7,180] days maturity, with a moneyness range of K/S ∈ [0.9, 1.1] across ten strikes. Each option price is computed using the Monte Carlo framework with 50, 000 simulations and a time-step of 1/20 th of a trading day.
The calibration process is done using the objective function defined in Eq. (23). It should be noted that as initial conditions we start with the true parameter vector, i.e Θ pert. initial = Θ M G true and for σ we start with the initial variance. We also pass the variance path at each time step.
Figure 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 line represents the Monte Carlo implied volatility and the dashed line is that of the leading order order perturbative solution. We pick an option maturity of 30 days. Figure 2: Test of the validity of the perturbation theory. These panels depict the ratios C 1 /(C 0 + C 1 ) in % where C 0 is the contribution to the price of the symmetrical solution and C 1 the leading order correction in perturbation theory. Clearly C 1 C 0 even when the volatility is large. This demonstrates the validity of the perturbative expansion. The four cases correspond respectively to panels (1,3,5,7) of Figure 1.