A lognormal type stochastic volatility model with quadratic drift
aa r X i v : . [ q -f i n . M F ] A ug A lognormal type stochastic volatility model with quadratic drift ∗ Peter Carr † and Sander Willems ‡ July 17, 2019
Abstract
This paper presents a novel one-factor stochastic volatility model where the instantaneousvolatility of the asset log-return is a diffusion with a quadratic drift and a linear dispersionfunction. The instantaneous volatility mean reverts around a constant level, with a speedof mean reversion that is affine in the instantaneous volatility level. The steady-state distri-bution of the instantaneous volatility belongs to the class of Generalized Inverse Gaussiandistributions. We show that the quadratic term in the drift is crucial to avoid momentexplosions and to preserve the martingale property of the stock price process. Using a con-veniently chosen change of measure, we relate the model to the class of polynomial diffusions.This remarkable relation allows us to develop a highly accurate option price approximationtechnique based on orthogonal polynomial expansions.
The popularity of the Heston (1993) model and, more generally, affine models (see e.g., Duffie et al.(2003)) for modeling stochastic volatility is in large part due to their analytical tractability. How-ever, there is abundant empirical evidence that favours non-affine models, in particular spec-ifications with a lognormal type of diffusion for the (instantaneous) volatility. For instance,Christoffersen et al. (2010) show that absolute changes in realized volatility are positively corre-lated with the volatility level and do not follow a Gaussian distribution. In contrast, the Hestonmodel implies that (instantaneous) changes in volatility should be Gaussian and independent ofthe volatility level. Changes in the log realized volatility, on the other hand, closely resemble anormal distribution, which motivates the use of a lognormal type of diffusion component in thevolatility process. Figure 1 reproduces these results using a 5-minute sub-sampled daily realizedvolatility measure for the S&P500 index from January 2000 until June 2019 and confirms thefindings of Christoffersen et al. (2010).Lognormal type stochastic volatility models are, however, particularly prone to problems such asmoment explosions and loss of the martingale property for the asset price, see e.g. Lions and Musiela(2007) and Andersen and Piterbarg (2007). These problems are caused by the fat right tail of thevolatility distribution, which can cause large spikes in the asset price. Having finite higher order ∗ We thank Damien Ackerer for helpful comments. † New York University, Tandon School of Engineering ‡ ´Ecole Polytechnique F´ed´erale de Lausanne (EPFL) and Swiss Finance Institute With lognormal type of diffusion we mean a diffusion σ t with d[ σ, σ ] t = ν σ t d t , for some ν > See Christoffersen et al. (2010) for the S&P500 index, Andersen et al. (2001) for individual stocks in the DJIAindex, and Andersen et al. (2001) for foreign exchange markets. −
87% in order for the asset price to have a finitefourth moment. In this paper, we propose a novel non-affine one-factor stochastic volatilitymodel featuring a diffusion with a quadratic drift function and a linear dispersion function forthe volatility process. The quadratic term has a negative coefficient in our model, which allowsfor a rapid reduction following an upward spike in the volatility. The linear dispersion functionproduces lognormal type innovations in the volatility and the quadratic term in the drift controlsundesirable side effects such as moment explosions and loss of martingality. Moreover, usingthe critical moment formula of Lee (2004), we show that a nonzero quadratic term in the driftallows to control both the small strike and large strike tail of the Black-Scholes implied volatilityskew. The volatility process in our model is stationary and has a Generalized Inverse Gaussian(GIG) distribution as steady-state distribution. The GIG distribution, which contains the in-verse Gaussian, hyperbolic, gamma, and inverse-gamma as special cases, has broad empiricalsupport for modeling stochastic volatility in stock returns, see for example Barndorff-Nielsen(1997), Eberlein (2001), Eberlein and Prause (2002), and Gander and Stephens (2007).Since our model is far from affine, tractability is not straightforward. If we set the quadraticterm in the drift of the volatility to zero, then our model fits in the class of polynomial diffusions,see e.g. Filipovi´c and Larsson (2016). This class of stochastic processes, which contains all affinediffusions as special cases, is characterized by the fact that their infinitesimal generator mapspolynomials to polynomials of the same degree or less. As a consequence, all conditional momentsof the log-asset price are available in closed form and European style derivatives on the assetprice can be priced using moment-based approximation methods, see e.g. Ackerer and Filipovi´c(2019). The polynomial property is lost, however, as soon as we have a nonzero quadratic termin the drift of the volatility. We circumvent this problem by introducing a change of measureunder which the polynomial property is recovered. Under the new measure, derivative pricesare given by the expectation of the discounted payoff multiplied by the Radon-Nikodym densityof the measure change. We show how to compute all joint conditional moments of the log-assetprice and log-Radon-Nikodym density in closed form under the new measure. An orthogonalpolynomial expansion technique in the spirit of Ackerer and Filipovi´c (2019) then allows us toefficiently price European style derivatives on the asset price.The remaining of this paper is structured as follows. Section 2 describes the model dynamicsand Section 3 analyzes the steady-state distribution of the volatility process. In Section 4 westudy the problem of moment explosions. Section 5 relates our model to the class of polynomialdiffusions, which is used in Section 6 to develop a derivative pricing approximation method.Section 7 contains a numerical study of the model and Section 8 concludes. All proofs andadditional technical results are collected in the Appendix. Bakshi et al. (2006) find empirical evidence for stochastic volatility models with nonlinear drift. In particular,they find a significantly negative coefficient on the quadratic term in the drift. For applications of polynomial processes in derivative pricing, see for example Filipovi´c et al. (2016),Ackerer and Filipovi´c (2016), Filipovi´c and Willems (2017), Ackerer et al. (2018), and Willems (2019). Volatility Level A b s o l u t e C hange s (a) Absolute change in volatility -7 -6.5 -6 -5.5 -5 -4.5 -4 -3.5 -3 -2.5 -2 Log Volatility Level A b s o l u t e C hange s (b) Absolute change in log volatility -4 -3 -2 -1 0 1 2 3 4 Standard Normal Quantiles -0.05-0.04-0.03-0.02-0.0100.010.020.030.040.05 D a t a Q uan t il e s (c) Changes in volatility QQ-plot -4 -3 -2 -1 0 1 2 3 4 Standard Normal Quantiles -1.5-1-0.500.511.52 D a t a Q uan t il e s (d) Changes in log volatility QQ-plot Figure 1: The top left (right) figure shows a scatter plot of the realized (log) volatility levelagainst the absolute change one day ahead, together with a least-squares regression line in red.The bottom left (right) figure shows a quantile-quantile plot of daily changes in realized (log)volatility. Realized volatilities are obtained from Oxford-Man Institute’s realized library using5-minute sub-sampled high-frequency returns on the S&P500 index from January 2000 untilJune 2019.
We consider a financial market modeled on a filtered probability space (Ω , F , F t , Q ), where Q isa risk-neutral probability measure. Henceforth E t [ · ] denotes the F t -conditional Q -expectation.Let S t denote the stock price and assume for simplicity zero interest rates and no dividendpayments. We specify the following Q -dynamics for the log-price x t = log( S t )d x t = − σ t d t + σ t ( ρ d W t + p − ρ d B t ) , (1)d σ t = ( R + R σ t )( R − σ t ) d t + νσ t d W t , (2) Alternatively, one can also think of S t as, for example, an interest rate variable (e.g., forward rate or swaprate) and replace Q by the appropriate pricing measure. R , R ≥ R , ν, σ > ρ ∈ [ − , x ∈ R , and W t , B t independent Q -Brownianmotions. The volatility process mean-reverts around a constant level R with a stochasticspeed of mean-reversion R + R σ t . Our model specification nests many existing models. Thelognormal SABR model of Hagan et al. (2002) arises when R = R = 0, in which case σ t issimply a geometric Brownian motion without drift. In this case, however, the volatility processis not mean-reverting, which is an important empirical feature. If we set R > R = 0,then σ t becomes a mean-reverting diffusion with affine drift and linear dispersion function,which we refer to as a linear diffusion. This type of model has been studied in Lewis (2000),Karasinski and Sepp (2012), Sepp (2014, 2016), Lee et al. (2016), and Ackerer and Filipovi´c(2019). For R = 0 and R >
0, equation (2) is known as the logistic diffusion and originatedin the context of modeling constrained population growth in biology, see e.g. Tuckwell and Koziol(1987). In the context of finance, the logistic diffusion has been used, for example, in a generalequilibrium model by Merton (1975) and in a stochastic volatility model by Hull and White(1987) and Lewis (2019). The following proposition shows that the model is well defined and that zero is an unattainableboundary for σ t . Proposition 2.1.
There exists a unique strong solution ( x t , σ t ) of (1) - (2) taking values in R × (0 , ∞ ) . From the proof of Proposition 2.1, it becomes clear that the non-negativity assumption R ≥ R <
0, then σ t blows up in finite time. Remark 2.2.
Although zero is a natural lower bound for the volatility process, we can generalize (2) by adding a lower bound σ ≥ as follows d σ t = ( R + R ( σ t − σ ))( R − ( σ t − σ )) d t + ν ( σ t − σ ) d W t , with σ > σ and the same restrictions on the other parameters as before. All the results wederive in our paper are easily adjusted to accommodate this generalization. σ t In this following proposition, we explicitly derive the steady-state distribution of σ t . Proposition 3.1.
If either R > , or R = 0 and R R > ν , then the process σ t has asteady-state distribution with density function π ( x ) ∝ x ξ − exp (cid:26) − R R ν x − R ν x (cid:27) , where ξ = − R − R R ν − . If R = 0 and R R ≤ ν , then σ t → a.s. as t → ∞ . Also related is the GARCH diffusion model of Nelson (1990) and Barone-Adesi et al. (2005), where σ t ismodeled as a diffusion with affine drift and linear dispersion function. Applying Itˆo’s lemma shows that thecorresponding volatility process also has a linear dispersion function, but it does not have an affine drift. Directlymodeling volatility seems more intuitive and provides a more natural interpretation for the model parameters. The deterministic version of this SDE was developed in the early 19th century by the Belgian mathematicianPierre Fran¸ois Verhulst to model population growth. Hull and White (1987) specify d σ t = a ( σ ∗ − σ t ) σ t d t + ξσ t d W t for some parameters a, σ ∗ , ξ >
0. ApplyingItˆo’s lemma shows that σ t follows a logistic diffusion: d σ t = a σ t ( σ ∗ − ξ − σ t ) d t + ξ σ t d W t . π integrates to oneis provided in the proof (see Section A.2). We distinguish four different cases, based on thevalues for R and R :1. For R = R = 0, the process σ t becomes a geometric Brownian motion without drift,which goes to zero almost surely as t → ∞ .2. For R > R = 0, the process σ t becomes a linear diffusion and we recover the inversegamma distribution as steady-state distribution, see e.g. Barone-Adesi et al. (2005). Thefirst moment of π is equal to R , so that R can be interpreted as the long-term level ofmean reversion. Higher order moments do not always exist because the inverse gammadistribution has a right tail with polynomial decay. Remark that π (0) = 0, regardless of R and R , due to the exponential decay of of the left tail.3. For R = 0 and R >
0, the process σ t becomes a logistic diffusion. If 2 R R ≤ ν ,then a similar behavior as in the first case occurs and σ t → t → ∞ .If 2 R R > ν , we recover the gamma distribution as steady-state distribution, whichhas finite moments of any order. In particular, the first moment equals R − ν R . Ashighlighted in Merton (1975) and Ewald and Yang (2007), we can therefore no longerinterpret R as the long-term level of mean reversion. Remark that π (0) = 0 if and onlyif R R > ν .4. For R , R >
0, the steady-state distribution has a gamma tail on the right and an inversegamma tail on the left. As a consequence, π has finite moments of any order and π (0) = 0.In particular, the first moment equals (see e.g., Jorgensen (1982)) √ R R √ R K ξ +1 (4 √ R R R ν − ) K ξ (4 √ R R R ν − ) , where K ξ denotes the modified Bessel function of the second kind. In general, the firstmoment of the steady-state distribution will not be exactly equal to R , so we can notinterpret R as the long-term level of mean-reversion. However, the difference is small forstandard parameterizations. In this section, we build on the general findings of Lions and Musiela (2007) to investigatemoment explosions of the stock price in our model. In order to use (1)-(2) for option pricing,we need S t to be a Q -martingale. The following proposition derives a necessary and sufficientcondition. Proposition 4.1. S t is a Q -martingale if and only if R ≥ ρν . If R = 0, then S t is a Q -martingale if and only if ρ ≤
0, which is a well known problem withthis type of model. While equity markets generally feature a negative correlation between stockreturns and volatility, other applications might require a positive correlation. Proposition 4.1shows that our model can accommodate a positive correlation if R is sufficiently large. In In Appendix B.1 we provide a tight lower bound (based on Jensen’s inequality) for the first moment of thesteady-state distribution that does not involve any special functions. R ≥ ν , then S t is always a Q -martingale in our model, regardless of ρ . Intuitively,the quadratic drift term has a stabilizing effect on the volatility because the speed of meanreversion becomes very large at high volatility levels.Andersen and Piterbarg (2007) highlight the importance for S t to have finite moments greaterthan one for pricing contracts with super-linear payoff, which occur frequently in interest ratederivatives markets. Examples include CMS swaps, in-arrears swaps, and Eurodollar futures.Moreover, when using Monte-Carlo simulations to find the price of a derivative, the payoff needsto have finite second order moment in order to derive confidence intervals on the Monte-Carloestimator with the central limit theorem. The following proposition derives a lower bound on R such that S t has finite moments of a given order. Proposition 4.2.
Let m ∈ R \ [0 , .1. If R > ν ( ρm + √ m − m ) , then E t [ S mT ] < ∞ , ∀ T > t. If R ≥ R R , then the statement is also true for R = ν ( ρm + √ m − m ) .2. If R < ν ( ρm + √ m − m ) , then E t [ S mT ] = ∞ , ∀ T > t.
In particular, if R = 0 and m >
1, then E t [ S mT ] is finite if and only if ρ ≤ − q m − m . A negativecorrelation has a dampening effect on the moments of the return process, however it must besufficiently negative in this case for higher moments to exist. For instance, already for m = 2 werequire ρ ≤ − .
71% in case R = 0, which can be quite restrictive. Proposition 4.2 shows thatthe quadratic term in the drift of σ t can take over the role of the negative correlation to stabilizethe moments of the return process, which allows the correlation to remain a free parameter.Remark also that for R = 0 and m <
0, we have E t [ S mT ] = ∞ , regardless of ρ .The seminal work of Lee (2004) relates moment explosions to the asymptotic behaviour of theBlack-Scholes implied volatility smile as a function of log-moneyness. Specifically, define thecritical moments m + ( T ) = sup { m : S mT < ∞} , m − ( T ) = inf { m : S mT < ∞} . (3)Remark that Proposition 4.2 implies in particular that the critical moments in our model donot depend on the time horizon T , so henceforth we omit the time argument and simply write m ± . Let σ BS ( T, x ) denote the Black-Scholes implied volatility of a European call option withtime-to-maturity T and strike price S e x . Using the formulation of Keller-Ressel (2011), thecritical moment formula of Lee (2004) stateslim sup x →−∞ σ BS ( T, x ) | x | = β ( − m − ) T and lim sup x →∞ σ BS ( T, x ) | x | = β ( m + − T , (4)where we define the decreasing function β : R + → [0 , , x − √ x + x − x ). The criticalmoments in our model can directly be computed using the result of Proposition 4.2, as shownin the following corollary. Note that Black-Scholes implied volatility only makes sense if S t is a Q -martingale, so we only consider the case R ≥ ρν , cf. Proposition 4.1.6 orollary 4.3. Suppose S t is a Q -martingale, i.e., R ≥ ρν .1. If | ρ | < , then m ± = 1 − R ν ρ ± q (1 − R ν ρ ) + 4(1 − ρ ) R ν − ρ ) . (5)
2. If ρ = 1 , then m − = −∞ and m + = R R ν − ν .3. If ρ = − , then m + = ∞ and m − = R − R ν − ν . The critical moment formula (4) and Corollary 4.3 give us important information about the tailbehaviour of x σ BS ( T, x ) T in our model. If | ρ | <
1, then the critical moments are finite,which implies asymptotically linear behaviour of σ BS ( T, x ) in x for all T >
0. The slope ofthe small and large strike tail is controlled by both ρ and R ν . For R = 0, we get in particular m + = (1 − ρ ) − and m − = 0. In this case, ρ only controls the slope of large strike tail, while theslope of the small strike tail is always equal to β (0) = 2. With R as a free parameter, we cantherefore more accurately capture both the small and the large strike tail of the Black-Scholesimplied volatility skew. Remark 4.4.
As noted by Lee (2004), the critical moment formula (4) can be useful to facilitatemodel calibration. Suppose we observe a Black-Scholes implied volatility skew for a range ofstrikes and a certain maturity
T > . From the smallest and largest strike, we can approximatelyinfer m − and m + , respectively. The parameters ρ and R ν can then be calibrated to these impliedcritical moments using (5) . This approach should be seen as a way to get good initial guessesfor ρ and R ν . We end this section with an additional result on the two extreme correlation cases.
Proposition 4.5.
1. If ρ = − and R ≥ R R , then S T ≤ S t exp (cid:26) σ t ν + R R ν ( T − t ) (cid:27) , ∀ T > t.
2. If ρ = 1 , R ≥ R R , and R ≥ ν , then S T ≥ S t exp (cid:26) − σ t ν − R R ν ( T − t ) (cid:27) , ∀ T > t.
For ρ = −
1, we know from Proposition 4.1 that S t is a Q -martingale and from Proposition 4.2that E t [ S mT ] < ∞ for all m > T > t . If moreover R ≥ R R , then Proposition 4.5shows that the stock price becomes bounded form above. Remark that this additional conditionis trivially satisfied when R = 0. For ρ = 1, the stock price is a Q -martingale if and only if R ≥ ν , see Proposition 4.1. From Proposition 4.2 we have in this case E t [ S mT ] < ∞ for all m < T > t . If moreover R ≥ R R , then the stock price has a lower bound strictlylarger than zero. Remark that these results are consistent with the critical moments derived inCorollary 4.3. Remark that in our model, the critical moments do not depend on the time horizon, while the implied criticalmoments will likely not be exactly equal for different option maturities, in which case we can for example averagethe implied critical moments across maturities. A polynomial diffusion in disguise
In this section we show how our model can be related to the class of polynomial diffusions usinga conveniently chosen change of measure.Define the process y t through the following stochastic differential equation (SDE)d y t = − z σ t d t + zσ t d W t , y = 0 , (6)with z = R ν . Fix a time horizon T > Q z through thefollowing Radon-Nikodym derivatived Q z d Q = e y T = e − z R T σ t d t + z R T σ t d W t . (7)Remark Q z = Q if R = 0. The following proposition shows that the change of measure is welldefined. Proposition 5.1.
The process e y t is a Q -martingale. Henceforth E zt [ · ] denotes the F t -conditional Q z -expectation. By Girsanov’s theorem we havethat W zt = W t − z Z t σ s d s and B zt = B t are independent Q z -Brownian motions. The Q z -dynamics of σ t becomesd σ t = ( R R + σ t ( R R − R )) d t + νσ t d W zt . (8)The quadratic term in the drift of σ t vanishes and σ t becomes a polynomial diffusion under Q z .Indeed, it has an affine drift and a linear dispersion function, so that its infinitesimal generatormaps polynomials to polynomials of the same degree or less. This allows us to compute all Q z -moments of σ t in closed form, which is informative about the Q z -distribution of σ t . Forderivative pricing purposes (see Section 6 for more details), we are not particularly interestedin the Q z -distribution of σ t . Instead, we are mainly interested in the Q -distribution of x t or,equivalently, in the joint Q z -distribution of x t and y t . The process ( x t , y t , σ t ) is not a polynomialdiffusion under Q z , because the drift of x t and y t contains a quadratic term σ t :d x t = ( zρ −
12 ) σ t d t + σ t ( ρ d W zt + p − ρ d B zt ) , d y t = 12 z σ t d t + zσ t d W zt , d σ t = ( R R + σ t ( R R − R )) d t + νσ t d W zt . However, by augmenting the state with σ t , we can see that ( x t , y t , σ t , σ t ) jointly becomes apolynomial diffusion under Q z since σ t has the following dynamicsd σ t = (2 R R σ t + σ t (2 R R − R + ν )) d t + 2 νσ t d W zt . This observation makes it possible to calculate all conditional Q z -moments of ( x t , y t , σ t ) in closedform. Before we do this, we first introduce some notation. Denote for m, n ∈ N by Pol m ( R n ) the8pace of polynomials on R n with total degree at most m . Define the subspace P m ⊂ Pol ( R )of trivariate polynomials as P m = (cid:8) ( x, y, z ) p ( x, y ) q ( z ) | p ∈ Pol m ( R ) , q ∈ Pol m − deg( p )) ( R ) (cid:9) , where deg( · ) denotes the total degree of a polynomial. The following lemma provides the di-mension of P m , i.e. the number of linearly independent polynomials in P m . Lemma 5.2.
The dimension of P m is d m = dim( P m ) = 13 m + 32 m + 136 m + 1 . The following proposition provides an explicit formula for the conditional Q z -moments of ( x t , y t , σ t ),which will be the cornerstone of the derivative pricing approximation method in Section 6. Proposition 5.3.
The infinitesimal generator G of the process ( x t , y t , σ t ) under Q z leaves P m invariant. That is, there exists a matrix G m ∈ R d m × d m , such that G H m = G m H m , where H m = ( h , . . . , h d m ) ⊤ denotes a vector of polynomial basis functions for P m . As a consequence,we have for any t ≤ T E zt [ H m ( x T , y T , σ T )] = e G m ( T − t ) H m ( x t , y t , σ t ) . (9)The matrix G m is straightforward to construct in practice by choosing H m to be a monomialbasis and then collecting terms according to their exponents in the vector of polynomials G H m ,see equation (22) in the Appendix. In this section we show how European style derivatives on the stock price e x T can efficiently becomputed using the available Q z -moments of ( x T , y T ). Consider a derivative on the stock price with payoff F (e x T ) at time T >
0, for some integrablepayoff function F . The price at time 0 is given by π = E [ F (e x T )] = E z [e − y T F (e x T )] . (10)The auxiliary process y t can therefore be interpreted as a stochastic discount rate under thenew measure. The positive correlation between y t and σ t provides a dampening effect on the‘discounted’ payoff under the new measure, which is the equivalent of the dampening effect ofthe quadratic drift term of σ t that disappeared with the measure change.The conditional Q -distribution of x T is not known, but we do know all the conditional Q z -moments of ( x T , y T ) thanks to the moment formula (9). Therefore, we can approximate thederivative price by approximating the function ( x, y ) e − y F (e x ) with a polynomial p n ∈ Pol n ( R ), for some n ∈ N . We would like the polynomial approximation to be most accurate forthe values that ( x T , y T ) is most likely to take under Q z , since they contribute the most to the9ight hand side of (10). This motivates the following least-squares approach to determine theapproximating polynomial p n = arg min p ∈ Pol n ( R ) Z R ( e − y F ( e x ) − p ( x, y )) w ( x, y ) d x d y, (11)where w is an auxiliary probability density function which proxies the unknown Q z -density of( x T , y T ). Put differently, p n is the orthogonal projection of ( x, y ) e − y F (e x ) on the space ofbivariate polynomials of total degree n or less in a weighted Hilbert function space with weight w . If we denote by ~p n ∈ R d n the vector representation of p n with respect to the basis H n , theoption price approximation becomes π ≈ π n = ~p ⊤ n e G n T H n ( x , y , σ ) . (12)In Section 6.2 we show how to choose w and in Section 6.3 we solve the optimization problemin (11). It remains to choose a good auxiliary density w . We use an approach that closely resembles theGaussian mixture specification of Ackerer and Filipovi´c (2019). Conditional on the trajectory { W zt , t ≤ T } , the Q z -density function of the random variable ( x T , y T ) can be formally writtenas ( x, y ) φ M T ,V T ( x ) δ ( y − y T ) , where φ M T ,V T denotes the density function of a Gaussian distribution with mean M T and vari-ance V T , δ denotes the Dirac delta function, and M T = x + ( −
12 + zρ ) Z T σ s d s + ρ Z T σ s d W zs , V T = (1 − ρ ) Z T σ s d s The true Q z -density function of ( x T , y T ) can therefore be expressed as( x, y ) E z [ φ M T ,V T ( x ) δ ( y − y T )] . We specify the auxiliary density as w ( x, y ) = K X k =1 w ( k ) φ m ( k ) ,v ( k ) ( x ) δ ( y − y ( k ) ) , (13)where m ( k ) , y ( k ) ∈ R , v ( k ) ∈ R + , w ( k ) ∈ [0 , k = 1 , . . . , K , are constants to be determinedsubject to P K k =1 w ( k ) = 1. The quadruplets ( w ( k ) , m ( k ) , v ( k ) , y ( k ) ) represent a discretization of the Q z -distribution of ( M T , V T , y T ) in K ≥ { W zt , t ≤ T } . Specifically, we use the IJK scheme of Kahl and J¨ackel(2006) with d ≥ σ t , We assume that w is such that the double integral in (11) is finite for all p ∈ Pol n ( R ). t , V t , and y t :ˆ σ n +1 = ˆ σ n + ( R R + ˆ σ n ( R R − R ))∆ + ν ˆ σ n √ ∆ Z n +1 + 12 ν ˆ σ n (∆ Z n +1 − ∆) , ˆ M n +1 = ˆ M n +1 + ( −
12 + zρ ) ˆ σ n +1 + ˆ σ n ρ ˆ σ t n √ ∆ Z n +1 , ˆ V n +1 = ˆ V n + (1 − ρ ) ˆ σ n +1 + ˆ σ n , ˆ y n +1 = ˆ y n + 12 z ˆ σ n +1 + ˆ σ n z ˆ σ n √ ∆ Z n +1 , where ∆ = Td is the step size, ( Z , . . . , Z d ) is a d -dimensional standard normal random vari-able, and ˆ σ = σ , ˆ M = M , ˆ V = V , ˆ y = y . If we are given K weighted samples( w ( k ) , Z ( k )1 , . . . , Z ( k ) d ), k = 1 . . . , K , of the random variable ( Z , . . . , Z d ), then by plugging theminto the above scheme we obtain the quadruplets ( w ( k ) , m ( k ) , v ( k ) , y ( k ) ). As highlighted byAckerer and Filipovi´c (2019), raw Monte-Carlo simulation with w ( k ) ≡ / K requires far toomany samples to produce an accurate approximation of the distribution. Instead, deterministicdiscretizations of the d -dimensional standard normal distribution, such as the quantization tech-niques of Pag`es and Printems (2003) or Gaussian cubature rules, are preferred in order to keep K small. In the numerical study in Section 7, we use the multivariate Gauss-Hermite quadra-ture method described in J¨ackel (2005) to obtain the weighted samples ( w ( k ) , Z ( k )1 , . . . , Z ( k ) d ).The advantage of Gauss-Hermite quadrature is that the tails of the distribution are accuratelycaptured, which is important for the stability of our approximation method as n , the totalpolynomial degree of the approximation, increases. Now that we have specified the auxiliary density function, we can solve the optimization problemin (11). Denote by B n = ( b , . . . , b N n ) ⊤ , N n = (cid:0) n +22 (cid:1) , a vector of polynomial basis functions forPol n ( R ). We can rewrite (11) as c n = arg min c ∈ R Nn Z R ( e − y F ( e x ) − c ⊤ B n ( x, y )) w ( x, y ) d x d y. (14) Proposition 6.1.
The unique solution of (14) is c n = D − f , with D i,j = Z R b i ( x, y ) b j ( x, y ) w ( x, y ) d x d y, f i = Z R e − y F (e x ) b i ( x, y ) w ( x, y ) d x d y, (15) for i, j = 1 , . . . , N n . Without loss of generality, we assume that B n is a monomial basis with b i ( x, y ) = x α i y β i , forexponents α i , β i ∈ N such that α i + β i ≤ n , i = 1 , . . . , N n . Plugging (13) in the expression for D i,j in (15) gives D i,j = K X k =1 w ( k ) ( y ( k ) ) β i + β j Z R x α i + α j φ m ( k ) ,v ( k ) ( x ) d x. α i + α j )-th moment of the univivariate Gaussian distri-bution, which is known in closed form. The elements of the vector f become f i = K X k =1 w ( k ) e − y ( k ) ( y ( k ) ) β i Z R F (e x ) x α i φ m ( k ) ,v ( k ) ( x ) d x. (16)In general, the integral in (16) has to be computed numerically, for example using Gauss-Hermite quadrature. For specific payoff functions, the integral can be computed in closed form.For example, the following proposition derives a recursive formula for the case of a Europeancall option. Proposition 6.2.
Suppose F ( x ) = (e x − K ) + , for some K > . The integral I ( k ) n = Z R F (e x ) x n φ m ( k ) ,v ( k ) ( x ) d x, satisfies the following recursion for n ≥ I ( k ) n = ( m ( k ) + v ( k ) ) I ( k ) n − + v ( k ) ( n − I ( k ) n − + Kv ( k ) J ( k ) n − ,J ( k ) n = m ( k ) J ( k ) n − + v ( k ) ( n − J ( k ) n − + p v ( k ) (log( K )) n − φ ( ξ ( k ) ) , with ξ ( k ) = m ( k ) − log( K ) √ v ( k ) , φ the standard normal density, and starting values I ( k ) − = 0 , I ( k )0 = e m ( k ) + v ( k ) Φ (cid:16) ξ ( k ) + p v ( k ) (cid:17) − K Φ( ξ ( k ) ) , J ( k ) − = 0 , J ( k )0 = Φ( ξ ( k ) ) , with Φ the standard normal cumulative distribution function. In this section we investigate the numerical accuracy of the option price approximation proposedin the previous section.We set the model parameters as R = R = 5, ν = 1, R = σ = 0 . ρ = − . x = 0. Theseare realistic parameters that produce a volatility process with strong mean-reversion and a highvolatility of volatility that can cause occasional spikes, see for example Figure 2 for a simulated(under Q ) trajectory. Consider a European call option with time-to-maturity T ∈ { / , / } and log-strike log( K ) ∈ {− . , , . } . Figure 3a and 3b plot the option price approximations π n for n ranging from 1 to 10. We set d = 1, and use the Gauss-Hermite quadrature rule toobtain a discretization of the univariate standard normal distribution in K = 15 points. As abenchmark, we also run a Monte-Carlo simulation with 10 sample paths. For all strikes andmaturities considered, π n converges to within the confidence bands of the Monte-Carlo estimatorwith n ≤
10. For n <
3, the pricing error is most noticeable for the low strike option (i.e., thein-the-money call). This is not surprising, since the true log-return distribution is negatively In Appendix B.2 we provide a simple formula for the moments of the Gaussian distribution. We use a quadratic polynomial approximation of the discounted payoff as a control variate to substantiallyreduce the variance of the Monte-Carlo estimator. To determine the polynomial approximation, we perform alinear regression with the simulated trajectories. This is similar to the polynomial approximation in Section 6,where we now use the simulated empirical distribution as auxiliary distribution. t t Figure 2: Simulated trajectory for σ t with parameters R = R = 5, R = σ = 0 . ν = 1.skewed ( ρ < K , as long as it is not too small. If K is chosen very small (say, K = 3), thenthe approximation blows up for larger n . Intuitively, for small K the auxiliary distribution w hasvery thin tails and the polynomial approximation of the discounted payoff will therefore only beaccurate over a small domain. Since the true probability distribution assigns considerable weightoutside of this domain, the polynomial approximation will blow up quickly. In Figure 3c and 3d,we do the same exercise with d = 2. For the auxiliary distribution, we use the Gauss-Hermitequadrature rule with 15 points in each dimension, which gives a total of K = 15 = 225 points. Compared to the case d = 1, the approximations converges faster to the true price. However,this comes at a computational cost because the number of discretization points in the auxiliarydistribution is much larger. We have introduced a new stochastic volatility model featuring a volatility process with aquadratic drift and a linear dispersion function. We have shown that the quadratic term inthe drift is important to control moment explosions in the stock price and, in particular, thesmall strike tail of the Black-Scholes implied volatility skew. The volatility process has a sta-tionary distribution that belongs to the class of Generalized Inverse Gaussian distributions,which arises frequently in the empirical literature on volatility modeling. In order to make themodel tractable, we introduced a change of measure such that the model fits into the class ofpolynomial diffusions, which opened the door to polynomial expansion methods to accuratelyapproximate option prices. Using the pruning method described in J¨ackel (2005), we can reduce the number of discretization points to185 by omitting the ‘corner’ points that carry a very small weight. n I m p li ed v o l a t ili t y log(K)=-0.1log(K)=0log(K)=0.1 (a) T = 1 / d = 1 n I m p li ed v o l a t ili t y log(K)=-0.1log(K)=0log(K)=0.1 (b) T = 2 / d = 1 n I m p li ed v o l a t ili t y log(K)=-0.1log(K)=0log(K)=0.1 (c) T = 1 / d = 2 n I m p li ed v o l a t ili t y log(K)=-0.1log(K)=0log(K)=0.1 (d) T = 2 / d = 2 Figure 3: Black-Scholes implied volatilities of approximated European call option prices withtime-to-maturity of one and two months for varying number of terms n in the series. Solidblue lines are Monte-Carlo estimators using 10 sample paths and the dashed blue lines arethe corresponding 99% confidence intervals. We use a quadratic polynomial approximation ofthe discounted payoff as a control variate to substantially reduce the variance of the Monte-Carlo estimator. The top row uses a single time step discretization to construct the auxiliarydistribution, while the bottom row uses two equidistant time steps.14 Proofs
A.1 Proof of Proposition 2.1
We start by showing that (2) has a unique (0 , ∞ )-valued solution. We denote by µ ( x ) = ( R + R x )( R − x ) and Σ( x ) = νx the drift and dispersion function of σ t , respectively. Since µ and Σ are polynomials, they are inparticular locally Lipschitz continuous. Hence, strong uniqueness holds for solutions of (2), seee.g. (Karatzas and Shreve, 1991, Theorem 2.5). The dispersion function satisfies a linear growthcondition | Σ( x ) | ≤ K (1 + | x | ) , for K ≥ ν . The drift function does not satisfy a linear growth condition, so the classicalexistence result of Itˆo (see e.g., (Karatzas and Shreve, 1991, Theorem 2.9)) does not apply.However, since the quadratic term in µ has a negative coefficient, µ does satisfy xµ ( x ) ≤ K (1 + | x | ) , for some K ≥
0. Hence, there exists a unique global solution to (2), cf. (Kloeden and Platen,1995, Chapter 4.5, p.135). It remains to verify that the solution stays in (0 , ∞ ), which we proofusing a comparison theorem. Consider the logistic diffusiond X t = ( − R X t + ( R R − R ) X t ) d t + νX t d W t , X = σ . (17)This SDE has a unique solution given by X t = Y t R R t Y s d s , with Y t = X e ( R R − R − ν ) t + νW t . Notice that, since R ≥
0, we have X t > t ≥
0. Using a comparison theorem(Ikeda and Watanabe, 1989, Chapter VI, Theorem 1.1) and R R ≥ σ t ≥ X t > t ≥ Remark A.1.
Notice that, if R < , then X t explodes in finite time and therefore σ t as well.The assumption R ≥ is therefore crucial to guarantee existence of a global solution. Next, we show that the stochastic integrals in (1) are well defined by showing that E t (cid:20)Z Tt σ s d s (cid:21) < ∞ , ∀ T ≥ t. To this end, consider the SDEd Z t = ( R R + ( R R − R ) Z t ) d t + νZ t d W t , Z = σ , which has Z t = Y t (cid:16) R R R t Y − s d s (cid:17) as unique solution. Using a comparison theorem(Ikeda and Watanabe, 1989, Chapter VI, Theorem 1.1) and R ≥ σ t ≤ Z t for all t ≥
0. Since Z t is a polynomial diffusion, it has finite moments of any order. Therefore wehave E t (cid:20)Z Tt σ s d s (cid:21) ≤ Z Tt E t (cid:2) Z s (cid:3) d s < ∞ , ∀ T ≥ t. From the comparison arguments used in first and second part of this proof, we also obtain thefollowing pathwise bounds on σ t X t ≤ σ t ≤ Z t , a.s. (18)15 .2 Proof of Proposition 3.1 Using the Fokker-Planck equation we have that π must satisfy the following second order linearODE ν d x [ x π ( x )] = dd x [( R + R x )( R − x ) π ( x )] . Making the Ansatz π = x γ exp { α x + βx } , for some constants α , β , and γ , and collecting termsgives π ( x ) = Cx ξ − exp (cid:26) − R R ν x − R ν x (cid:27) , where C is a constant to be determined such that π integrates to one. Remark A.2.
The motivation for this Ansatz comes from the special cases R = 0 and R = 0 .If R = 0 , then σ t is a GARCH diffusion, which is known to have the inverse gamma distributionas steady-state distribution. If R = 0 , then σ t is a logistic diffusion, which is known to havethe gamma distribution as steady-state distribution if R R > ν . Therefore, the steady-statedistribution of σ t must contain the gamma and inverse gamma distribution as special cases. If R , R >
0, then π is a generalized inverse Gaussian distribution and the normalizationconstant becomes (see e.g., Jorgensen (1982)) C = (cid:16) R R R (cid:17) ξ/ K ξ (4 √ R R R ν − ) , where K ξ denotes the modified Bessel function of the second kind. If R = 0 and 2 R R > ν ,then π is a gamma density and C therefore becomes C = (cid:0) R ν (cid:1) ξ Γ( ξ ) , where Γ denotes the gamma function. Note that the condition 2 R R > ν is equivalent to ξ > R = 0 and 2 R R ≤ ν , recall from (18) that we have almost surely the following upperbound σ t ≤ Z t = σ e ( R R − ν ) t + νW t . If 2 R R ≥ ν , then e ( R R − ν ) t + νW t → t → ∞ and therefore we have σ t → t → ∞ . A.3 Proof of Proposition 4.1
By Theorem 2.4(i) in Lions and Musiela (2007), S t is a Q -martingale iflim x →∞ ρνx + ( R + R x )( R − x ) x < ∞ . x →∞ ρνx + ( R + R x )( R − x ) x = lim x →∞ ( ρν − R ) x + ( R R − R ) x + R R x = lim x →∞ ( ρν − R ) x + R R − R . Therefore, S t is a Q -martingale if R ≥ ρν . Indeed, if R > ρν the limit is −∞ and if R = ρν the limit is R R − R .Conversely, by Theorem 2.4(ii) in Lions and Musiela (2007), S t is not a Q -martingale iflim x →∞ ρνx + ( R + R x )( R − x ) φ ( x ) > , for some smooth, positive, and increasing function φ such that R ∞ ǫ φ ( x ) d x < ∞ , ǫ >
0. Choosing φ ( x ) = x gives lim x →∞ ρνx + ( R + R x )( R − x ) x = ρν − R . Therefore, S t is not a Q -martingale if R < ρν . A.4 Proof of Proposition 4.2
1. By Theorem 2.5 in Lions and Musiela (2007), we need to show that there exists an A ≥ lim x →∞ − A ν x − A [ mρνx + ( R + R x )( R − x )] − m − m x > −∞ . (19)The limit can be rewritten aslim x →∞ − A ν + x − A [ mρνx + ( R + R x )( R − x )] − m − m x = lim x →∞ (cid:18) − A ν − A [ mρν − R ] − m − m (cid:19) x − A ( R R − R ) x − AR R . (20)Define the parabola f ( u ) = − u ν − u [ mρν − R ] − m − m . If mρν − R < − ν p m − m, then f has two distinct positive roots A ± = mρν − R ± p ( mρν − R ) − ν ( m − m ) − ν . From (20) it becomes clear that if we pick A ∈ ( A − , A + ), then (19) is satisfied.If mρν − R = − ν p m − m, The paper of Lions and Musiela (2007) contains some typos that are relevant for the derivation of this proof.Specifically, in equation (26), the function β should be defined as β ( x ) = mρµ ( x ) x + b ( x ) instead of β ( x ) = mµ ( x ) x + b ( x ). In equation (28), the µ in the last term has to be replaced by m . f only has a single root A = mρν − R − ν = √ m − mν > . For all other values, f will be negative. In other words, any value other than A will make(20) equal to −∞ . It remains to check what happens to (20) for A = A lim x →∞ f ( A ) x − A ( R R − R ) x − A R R = lim x →∞ − A ( R R − R ) x − A R R . Therefore, the limit will be larger than −∞ if R ≥ R R .2. Follows directly from Theorem 2.6 in Lions and Musiela (2007) andlim x →∞ νxx = ν, lim x →∞ ( R + R x )( R − x ) x = − R . A.5 Proof of Corollary 4.3
Define the function f ( m ) = ρm + √ m − m on R \ (0 , f is increasingon [1 , ∞ ) with f (1) = ρ and decreasing on ( −∞ ,
0] with f (0) = 0.If | ρ | <
1, then standard calculations show lim m →±∞ f ( m ) = ∞ , so by Proposition 4.2 there mustbe a critical moment both in ( −∞ ,
0] and in [1 , ∞ ). In order to find the critical moments, wehave to solve the equation R ν − ρm = p m − m, m ∈ R \ (0 , . (21)Squaring both sides shows that a critical moment m has to satisfy p ( m ) = m (1 − ρ ) + (2 R ν ρ − m − R ν = 0 . If | ρ | <
1, then p is a convex parabola with p (0) = − R ν ≤ p (1) = − ρ − R ν + 2 R ν ρ ≤ R ≥ νρ . Therefore, p has tworeal roots m ± = 1 − R ν ρ ± q (1 − R ν ρ ) + 4(1 − ρ ) R ν − ρ )with m − ≤ m + ≥
1. It is directly verified that m ± solves (21).If ρ = −
1, then lim m →∞ f ( m ) = − <
0. Therefore, R ν ≥ f ( m ) for all m ≥
1, so that m + = ∞ .Since lim m →−∞ f ( m ) = ∞ , there will be a critical moment in ( −∞ ,
0] and it is given by the singleroot of p : m − = R − R ν − ν ≤ . Similarly, if ρ = 1, then lim m →−∞ f ( m ) = . Since we assume R ≥ ν in this case, we have inparticular R ν ≥ f ( m ) for all m ≤
0, so that m − = −∞ . Since lim m →∞ f ( m ) = ∞ , there will be acritical moment in [1 , ∞ ) and it is given by the single root of p : m + = R R ν − ν ≥ , where the inequality follows from the assumption R ≥ ρν = ν .18 .6 Proof of Proposition 4.5
1. Solving (1) gives S T = S t exp (cid:26) − Z Tt σ s d W s − Z Tt σ s d s (cid:27) = S t exp (cid:26) − ν (cid:18) σ T − σ t − Z Tt ( R + R σ s )( R − σ s ) d s (cid:19) − Z Tt σ s d s (cid:27) = S t exp (cid:26) − ν (cid:18) σ T − σ t − R R ( T − t ) − ( R R − R ) Z Tt σ s d s (cid:19) − (cid:18)
12 + R ν (cid:19) Z Tt σ s d s (cid:27) ≤ S t exp (cid:26) σ t ν + R R ν ( T − t ) (cid:27) .
2. Similarly as in the first part, solving (1) gives S T = S t exp (cid:26)Z Tt σ s d W s − Z Tt σ s d s (cid:27) = S t exp (cid:26) ν (cid:18) σ T − σ t − R R ( T − t ) − ( R R − R ) Z Tt σ s d s (cid:19) − (cid:18) − R ν (cid:19) Z Tt σ s d s (cid:27) ≥ S t exp (cid:26) − σ t ν − R R ν ( T − t ) (cid:27) . A.7 Proof of Proposition 5.1
The dynamics of X t := zσ t becomesd X t = ( R + R σ t )( R z − X t ) d t + νX t d W t = ( R + νX t )( R R ν − X t ) d t + νX t d W t = b ( X t ) d t + νX t d W t , where we defined the function b ( x ) = ( R + νx )( R R ν − x ). The dynamics of e y t becomesde y t = X t e y t d W t . We conclude by Theorem 2.4(i) in Lions and Musiela (2007) that e y t is a Q -martingale, sincelim x → + ∞ b ( x ) + νx x = R R − R < ∞ . A.8 Proof of Lemma 5.2
It is well known that dim(Pol m ( R n )) = (cid:0) m + nn (cid:1) . As a consequence, the dimension of the set ofpolynomials in R n with total degree exactly equal to k ∈ N is (cid:18) k + nn (cid:19) − (cid:18) k − nn (cid:19) = k + 1 , (cid:0) n − n (cid:1) = 0. We therefore getdim( P m ) = m X k =0 ( k + 1)(1 + 2( m − k ))= (2 m + 1)( m + 1) + (2 m − m X k =0 k − m X k =0 k = (2 m + 1)( m + 1) + 12 m ( m + 1)(2 m − − m ( m + 1)(2 m + 1)= 13 m + 32 m + 136 m + 1 . A.9 Proof of Proposition 5.3
Without loss of generality, we can use a monomial basis for P m . A generic element in this basiscan be represented as x α y β z γ , with α, β, γ ∈ N , α + β ≤ m and γ ≤ m − α − β ). Applyingthe Q z -generator G of ( x t , y t , σ t ) to this monomial gives G x α y β z γ = α ( zρ −
12 ) x α − y β z γ +2 + β z x α y β − z γ +2 + γR R x α y β z γ − + γ ( R R − R ) x α y β z γ + 12 α ( α − x α − y β z γ +2 + 12 β ( β − z x α y β − z γ +2 + 12 γ ( γ − ν x α y β z γ + αβzρx α − y β − z γ +2 + αγνρx α − y β z γ +1 + βγzνx α y β − z γ +1 . (22)It is readily verified by inspecting each of the above monomials that G x α y β z γ ∈ P m . A.10 Proof of Proposition 6.1
The optimization problem in (14) is a convex quadratic programming problem. The first orderconditions become2 Z R ( e − y F ( e x ) − c ⊤ B n ( x, y )) b i ( x, y ) w ( x, y ) d x d y = 0 , i = 1 , . . . , N n . Re-arranging terms we get N n X j =1 Z R c j b j ( x, y ) b i ( x, y ) w ( x, y ) d x d y = Z R e − y F ( e x ) b i ( x, y ) w ( x, y ) d x d y. In matrix notation this becomes Dc = f. Remark that the matrix D positive-definite by construction and therefore invertible.20 .11 Proof of Proposition 6.2 To lighten the notation, we suppress the superscript ( k ) throughout the proof. Using the identity xφ ( x ) = − φ ′ ( x ) and integrating by parts gives I n = Z R (e x − K ) + x n φ m,v ( x ) d x = Z R (e m + √ vx − K ) + ( m + √ vx ) n φ ( x ) d x = Z ∞− ξ (e m + √ vx − K )( m + √ vx ) n φ ( x ) d x = mI n − + √ v Z ∞− ξ (e m + √ vx − K )( m + √ vx ) n − xφ ( x ) d x = mI n − − √ v Z ∞− ξ (e m + √ vx − K )( m + √ vx ) n − φ ′ ( x ) d x = mI n − − √ v (cid:20) −√ v ( n − I n − − √ v Z ∞− ξ e m + √ vx ( m + √ vx ) n − φ ′ ( x ) d x (cid:21) = ( m + v ) I n − + v ( n − I n − + Kv Z ∞− ξ ( m + √ vx ) n − φ ( x ) d x. Define J n = R ∞− ξ ( m + √ vx ) n φ ( x ) d x . Similarly as for I n , we use integration by parts to derivethe following recursion for J n J n = mJ n − − √ v Z ∞− ξ ( m + √ vx ) n − φ ′ ( x ) d x = mJ n − − √ v (cid:20) − ( m − √ vξ ) n − φ ( − ξ ) − ( n − √ v Z ∞− ξ ( m + √ vx ) n − φ ′ ( x ) d x (cid:21) = mJ n − + v ( n − J n − + √ v (log( K )) n − φ ( ξ ) . For the starting values of the recursion, we have J = Z ∞− xi φ ( x ) d x = Φ( ξ ) ,I = e m + v Φ (cid:0) ξ + √ v (cid:1) − K Φ( ξ ) . We omit the full derivation of I since it is very similar to computing the price of a Europeancall option in the Black-Scholes model. B Auxiliary results
B.1 Lower bound on first moment of steady-state density
Suppose the assumptions of Proposition 3.1 are satisfied so that a non-trivial steady-state distri-bution exists, and suppose furthermore that R >
0. We introduce for simplicity the followingnotation α = − R R ν , β − − R ν . Z ∞ xπ ( x ) d x = R ∞ x ξ exp { α x + βx } d x R ∞ x ξ − exp { α x + βx } d x = R ∞ x ξ exp { α x } d (cid:16) exp { βx } β (cid:17)R ∞ x ξ − exp { α x + βx } d x = − ξβ + αβ R ∞ x ξ − exp { α x + βx } d x R ∞ x ξ − exp { α x + βx } d x = − ξβ + αβ Z ∞ x − π ( x ) d x, where we have used integration by parts on the integral in the numerator. Remark now that x /x is strictly convex for x >
0, so we have by Jensen’s inequality Z ∞ xπ ( x ) d x ≥ − ξβ + αβ (cid:18)Z ∞ xπ ( x ) d x (cid:19) − . Remark that the inequality is strict if and only if R = 0. If we denote µ = R ∞ xπ ( x ) d x , weobtain the following quadratic inequality for µµ + ξβ µ − αβ ≥ . By solving the roots of the parabola, we see that the above inequality can only be true if µ ≥ − ξβ + s ξ β + 4 αβ ! . B.2 Moments of the Gaussian distribution
Suppose we want to compute all moments of a univariate Gaussian distribution with mean µ ∈ R and variance σ >
0. Consider the Gaussian process X t defined through the followingSDE d X t = µ d t + σ d W t , X = 0 . where W t is a standard Brownian motion. The solution X at time 1 has a Gaussian distributionwith mean µ and variance σ . Applying the infinitesimal generator G of X t to a monomial x n gives G x n = nx n − µ + 12 n ( n − x n − σ . Therefore, if we define G n ∈ R ( n +1) × ( n +1) as G n = µ σ µ
00 3 σ µ · · · σ n ( n − nµ , G (cid:0) x x · · · x n (cid:1) ⊤ = G n (cid:0) x x · · · x n (cid:1) ⊤ . By definition of the generator, we get the following simple formula for the Gaussian mo-ments E [ (cid:0) X X · · · X n (cid:1) ⊤ ] = e G n (cid:0) · · · (cid:1) ⊤ . eferences Ackerer, D. and D. Filipovi´c (2016). Linear credit risk models. arXiv preprint arXiv:1605.07419 .Ackerer, D. and D. Filipovi´c (2019). Option pricing with orthogonal polynomial expansions.
Mathematical Finance , Forthcoming.Ackerer, D., D. Filipovi´c, and S. Pulido (2018). The Jacobi stochastic volatility model.
Financeand Stochastics 22 (3), 667–700.Andersen, L. B. and V. V. Piterbarg (2007). Moment explosions in stochastic volatility models.
Finance and Stochastics 11 (1), 29–50.Andersen, T. G., T. Bollerslev, F. X. Diebold, and H. Ebens (2001). The distribution of realizedstock return volatility.
Journal of Financial Economics 61 (1), 43–76.Andersen, T. G., T. Bollerslev, F. X. Diebold, and P. Labys (2001). The distribution of realizedexchange rate volatility.
Journal of the American Statistical Association 96 (453), 42–55.Bakshi, G., N. Ju, and H. Ou-Yang (2006). Estimation of continuous-time models with anapplication to equity volatility dynamics.
Journal of Financial Economics 82 (1), 227–249.Barndorff-Nielsen, O. E. (1997). Normal inverse Gaussian distributions and stochastic volatilitymodelling.
Scandinavian Journal of statistics 24 (1), 1–13.Barone-Adesi, G., H. Rasmussen, and C. Ravanelli (2005). An option pricing formula for theGARCH diffusion model.
Computational Statistics & Data Analysis 49 (2), 287–310.Christoffersen, P., K. Jacobs, and K. Mimouni (2010). Volatility dynamics for the S&P500:evidence from realized volatility, daily returns, and option prices.
The Review of FinancialStudies 23 (8), 3141–3189.Duffie, D., D. Filipovi´c, W. Schachermayer, et al. (2003). Affine processes and applications infinance.
The Annals of Applied Probability 13 (3), 984–1053.Eberlein, E. (2001). Application of generalized hyperbolic L´evy motions to finance. In
L´evyprocesses , pp. 319–336. Springer.Eberlein, E. and K. Prause (2002). The generalized hyperbolic model: financial derivatives andrisk measures. In
Mathematical Finance - Bachelier Congress 2000 , pp. 245–267. Springer.Ewald, C.-O. and Z. Yang (2007). Geometric mean reversion: formulas for the equilibriumdensity and analytic moment matching.
Available at SSRN 999561 .Filipovi´c, D., E. Gourier, and L. Mancini (2016). Quadratic variance swap models.
Journal ofFinancial Economics 119 (1), 44–68.Filipovi´c, D. and M. Larsson (2016). Polynomial diffusions and applications in finance.
Financeand Stochastics 20 (4), 931–972.Filipovi´c, D. and S. Willems (2017). A term structure model for dividends and interest rates.
Working Paper . 24ander, M. P. and D. A. Stephens (2007). Stochastic volatility modelling in continuous timewith general marginal distributions: Inference, prediction and model selection.
Journal ofStatistical Planning and Inference 137 (10), 3068–3081.Hagan, P. S., D. Kumar, A. S. Lesniewski, and D. E. Woodward (2002). Managing smile risk.
The Best of Wilmott 1 , 249–296.Heston, S. L. (1993). A closed-form solution for options with stochastic volatility with applica-tions to bond and currency options.
The Review of Financial Studies 6 (2), 327–343.Hull, J. and A. White (1987). The pricing of options on assets with stochastic volatilities.
TheJournal of Finance 42 (2), 281–300.Ikeda, N. and S. Watanabe (1989).
Stochastic Differential Equations and Diffusion Processes (2nd ed.), Volume 24. Elsevier.J¨ackel, P. (2005). A note on multivariate gauss-hermite quadrature.
Technical report .Jorgensen, B. (1982).
Statistical Properties of the Generalized Inverse Gaussian Distribution ,Volume 9. Springer Science & Business Media.Kahl, C. and P. J¨ackel (2006). Fast strong approximation monte carlo schemes for stochasticvolatility models.
Quantitative Finance 6 (6), 513–536.Karasinski, P. and A. Sepp (2012). Beta stochastic volatility model.
Risk , 66–71.Karatzas, I. and S. Shreve (1991).
Brownian Motion and Stochastic Calculus (2nd ed.). Springer-Verlag.Keller-Ressel, M. (2011). Moment explosions and long-term behavior of affine stochastic volatil-ity models.
Mathematical Finance: An International Journal of Mathematics, Statistics andFinancial Economics 21 (1), 73–98.Kloeden, P. E. and E. Platen (1995).
Numerical Solution of Stochastic Differential Equations ,Volume 23. Springer Science & Business Media.Lee, G., Z. Zhu, et al. (2016). Switching to non-affine stochastic volatility: A closed-formexpansion for the Inverse Gamma model. Technical report.Lee, R. W. (2004). The moment formula for implied volatility at extreme strikes.
MathematicalFinance 14 (3), 469–480.Lewis, A. L. (2000).
Option Valuation Under Stochastic Volatility: With Mathematica Code .Finance Press.Lewis, A. L. (2019). Exact solutions for a GBM-type stochastic volatility model having astationary distribution. arXiv preprint arXiv:1809.08635 .Lions, P.-L. and M. Musiela (2007). Correlations and bounds for stochastic volatility models. In
Annales de l’Institut Henri Poincare (C) Non Linear Analysis , Volume 24, pp. 1–16. Elsevier.Merton, R. C. (1975). An asymptotic theory of growth under uncertainty.
The Review ofEconomic Studies 42 (3), 375–393. 25elson, D. B. (1990). ARCH models as diffusion approximations.
Journal of Econometrics 45 (1-2), 7–38.Pag`es, G. and J. Printems (2003). Optimal quadratic quantization for numerics: the Gaussiancase.
Monte Carlo Methods and Applications 9 (2), 135–165.Sepp, A. (2014). Empirical calibration and minimum-variance delta under log-normal stochasticvolatility dynamics.
Available at SSRN 2387845 .Sepp, A. (2016). Log-normal stochastic volatility model: Affine decomposition of moment gen-erating function and pricing of vanilla options.Tuckwell, H. C. and J. A. Koziol (1987). Logistic population growth under random dispersal.
Bulletin of Mathematical Biology 49 (4), 495–506.Willems, S. (2019). Asian option pricing with orthogonal polynomials.