Anomalous diffusions in option prices: connecting trade duration and the volatility term structure
AAnomalous diffusions in option prices: connecting trade durationand the volatility term structure ∗ Antoine Jacquier † Lorenzo Torricelli ‡ April 13, 2020
Abstract
Anomalous diffusions arise as scaling limits of continuous-time random walks (CTRWs)whose innovation times are distributed according to a power law. The impact of a non-exponential waiting time does not vanish with time and leads to different distribution spreadrates compared to standard models. In financial modelling this has been used to accom-modate for random trade duration in the tick-by-tick price process. We show here thatanomalous diffusions are able to reproduce the market behaviour of the implied volatilitymore consistently than usual L´evy or stochastic volatility models. Two distinct classes ofunderlying asset models are analyzed: one with independent price innovations and waitingtimes, and one allowing dependence between these two components. These models capturethe well-known paradigm according to which shorter trade duration is associated with higherreturn impact of individual trades. We fully describe these processes in a semimartingalesetting leading to no-arbitrage pricing formulae, study their statistical properties, and inparticular observe that skewness and kurtosis of asset returns do not tend to zero as timegoes by. We finally characterize the large-maturity asymptotics of Call option prices, andfind that the convergence rate to the spot price is slower than in standard L´evy regimes,which in turn yields a declining implied volatility term structure and a slower time decay ofthe skew.
Keywords : Anomalous diffusions, volatility skew term structure, derivative pricing, CTRWs,inverse L´evy subordinators, time changes, L´evy processes, subdiffusions, Beta distribution, tri-angular arrays.
In quantitative finance, models of asset returns typically evolve according to Itˆo diffusions orL´evy-type models. From a microstructural point of view, these can be seen as scaling limits ∗ The authors would like to thank Mark Meerschaert and Peter Straka for the valuable suggestions. † Department of Mathematics, Imperial College London, and Alan Turing Institute. ‡ Department of Economics and Management. Email: [email protected]. a r X i v : . [ q -f i n . P R ] A p r f continuous-time random walks (CTRWs) with exponentially distributed inter-arrival times.Instead, time changing CTRWs to a renewal process whose waiting times obey a power lawyields, in the scaling limit, an anomalous diffusion , namely a space-time propagation processwhere the particle spreads at a rate different from linear, which is observed in the classicaldiffusive case. The use of anomalous diffusions in financial models was pioneered by Mainardiet al. (2000) and Scalas et al. (2000), and they have proved useful to capture memory effects,trade idle time, and other microstructural price features exhibited by high-frequency time series.However, applications of anomalous diffusions for continuous-time option pricing have sofar been scarce. The sub-diffusive Black-Scholes model was introduced in Magdziarz (2009) tocapture asset staleness and periods of trade inactivity, but implications on option pricing andimplied volatilities were not illustrated. Cartea and Meyer-Brandis (2010) analysed the volatilitysurface of a CTRW whose innovation times are distributed according to a Mittag-Leffler hazardfunction, produced explicit option pricing formulae, and provided evidence that the long-termskewness and smile can be captured.We show here how anomalous diffusions in equity returns can also capture the long-termbehaviour of the implied volatility surface. Specifically, we argue that the persistence of a slowlydecaying volatility skew can be explained by postulating the survival of trade durations effects atlonger maturities. We consider returns and innovation time random walks which converge in thescaling limit to a pair of L´evy processes, one of which is a subordinator. According to Becker-Kern et al. (2004); Meerschaert and Scheffler (2008, 2010); Straka and Henry (2011); Jurlewiczet al. (2012), the associated CTRW time-changed with the renewal process of the innovationtimes converges to an anomalous diffusion which can be represented as a time-changed L´evyprocess. One appealing feature is that both analytical formulae for the Laplace transforms (inthe time variable) of the characteristic function of this limit and integral expressions for thedensity functions (in terms of the L´evy measures) are known .We analyse two distinct classes of anomalous diffusion models. The first is the purely subd-iffusive L´evy model (SL), where the CTRW limiting diffusion consists of a L´evy process time-changed by an independent inverse-stable subordinator. Several instances of such models havebeen already investigated in the literature. In terms of the generating fractional Fokker-Planckequations such a class has been investigated in Cartea and del-Castillo-Negrete (2007). The par-ticular case where the parent L´evy process is a Brownian motion was introduced in Magdziarz(2009); the compound Poisson case in Cartea and Meyer-Brandis (2010). A general treatmentwhen the driving noise is a generic L´evy process has been recently provided in Torricelli (2020).Moreover the classical models in Mainardi et al. (2000), Scalas et al. (2000) also admit a rep-resentation of the SL form: we revisit such models as stochastic time changes, well suited toolsfor option pricing purposes.The time change representation of subdiffusive models also paves the way for our secondsecond class of models, developing an idea that originally appeared in (Becker-Kern et al.,2004, Example 2.8) and (Jurlewicz et al., 2012, Examples 5.4). This asset price evolution2ealistically incorporates the dependence between the L´evy parent returns generating processand the inverse-stable subordinator modelling the trades waiting time. We call it the modelwith dependent returns and trade duration (DRD).Apart from being natural outcomes of time-changed random walk tick-by-tick price models,these two models are strongly supported by the econometric analysis by Engle (2000) and Dufourand Engle (2000), and confirmed in numerous empirical studies later on. The evidence is thattrading activity is inversely correlated with price impact, i.e. the ‘volatility’ of the asset price:the fewer the trades (longer duration), the more sluggish the price innovations; conversely,intense trading (short duration) is associated with higher price excursions. Remarkably, thisprinciple is captured in our setting.We describe such equity models in a semimartingale dynamic setting, leading to no-arbitragepricing relations under appropriate equivalent risk neutral measures. Using the results of Ju-rlewicz et al. (2012) on the Fourier-Laplace transforms of anomalous diffusions, we furtherprovide familiar Parseval-Plancherel formulae for option prices in the spirit of Lewis (2001).Additionally, we study the moments and serial correlation properties of the model and showthat skewness and kurtosis of the asset returns in the DRD model converge for large times, anddo not vanish, contrary to L´evy models, leading in particular to profound differences on thelong-term volatility smile.Finally we characterize the large-maturity behaviour of Call options and find that the con-vergence rate is much slower than in standard L´evy or stochastic volatility regimes. We uncovera relationship according to which a declining implied volatility level implies a slowly decayingskew, at least compared to that of L´evy and exponentially affine models. But we find that a(slowly) vanishing volatility level is a defining feature of these models, due to long-maturityprices converging much slower than in standard models. Ultimately, for the DRD model weshow that the vanishing rate of the skew is slower than the usual 1 /T , in line with market data.As illustrated in the calibration in Section 8, the practical importance of anomalous diffusionmodel is that the ‘duration parameter’ β improves the cross-sectional fit to multiple maturitiescompared to a L´evy model, while having virtually no impact on the short-maturity calibration.This justifies the interpretation of β as a long term skew component.We believe the contribution of this work to be manifold. We establish an explicit structuralconnection between trade duration and skew persistence; we introduce an analytical model thataccounts for trades duration and dependence between trade waiting times and returns, consistentwith the Econometrics literature; we systematically unify the treatment of SL models underthe umbrella of a single time-changed representation and the corresponding analytic pricingformulae; finally we extend the analysis of the ‘Beta-time’ process in Meerschaert and Scheffler(2004) and Jurlewicz et al. (2012), providing its moments and statistical properties through itstime-changed representation.In Section 2 we introduce fundamental building blocks and some useful notations. In Sec-tion 3 we introduce the CTRWs components of the base tick-by-tick model and the convergence3heorem leading to their limiting continuous-time versions. The anomalous diffusions are intro-duced in Section 4, together with their analytical properties and time-changed semimartingalerepresentations, while their statistical properties are characterized in Section 5. In Section 6 weshow how to construct equivalent pricing measures, and provide an integral price representationfor European Call option prices. This allows us to study in Section 7 the structure of the corre-sponding implied volatility, with a particular emphasis on its large-maturity properties. Finallyin Section 8, we numerically highlight interesting features of the SL and DRD models, and showthat both models allow for a good fit to market data. We follow here (Kyprianou, 2014, Chapter 1). In a market filtration (Ω , F , ( F t ) t ≥ , P ), a L´evyprocess X is uniquely characterized by its L´evy exponent, namely the function ψ X : C → C defined via the relation E (cid:2) e − i zX t (cid:3) = exp ( − tψ X ( z )), and given explicitly by the L´evy-Khintchineformula ψ X ( z ) = i zµ + z σ − (cid:90) R ( e − i zx − zx | x | < ) ν ( dx ) , (2.1)where µ ∈ R , σ ≥
0, and ν is a measure concentrated on R \ { } such that (cid:82) R (1 ∧ x ) ν ( dx ) isfinite. In order to guarantee some minimal properties of the asset pricing model (existence ofthe first moment), we always assume that (cid:90) | x | > e x ν ( dx ) < ∞ . (2.2)A subordinator L is an almost surely non-decreasing L´evy process, with L´evy measure ν L supported on (0 , ∞ ), and the L´evy-Khintchine representation for its Laplace exponent definedvia the relation E (cid:2) e − sX t (cid:3) = exp ( − tφ L ( s )) simplifies to φ L ( s ) = sµ − (cid:90) ∞ ( e − su − ν L ( du ) , (2.3)for µ >
0, and where (cid:82) ∞ uν L ( du ) < ∞ . A bivariate L´evy process ( X, L ), with L a subordinator,has joint Fourier-Laplace transform E [ e − i zX t − sL t ] = exp ( − tψ X,L ( z, s )) of the form ψ X,L ( z, s ) = i zµ X + sµ L + z σ − (cid:90) R (cid:90) ∞ (cid:0) e − i zx − su − zx | x | < (cid:1) ν X,L ( dx, du ) , (2.4)with L´evy-Laplace triplet (( µ X , µ T ) , σ, ν X,T ).For a process Y , we denote Y t − denotes the random variable of the left limits, Y t + that ofthe right limits, and with Y − = ( Y t − ) t ≥ and Y + = ( Y t + ) t ≥ the corresponding processes. If Y is a L´evy process, stochastic continuity implies that Y = Y − = Y + up to a modification . A process Y is said to be a modification of a process X if ∀ t , P ( X t = Y t ) = 1. t, ∞ ) of L is the random variable H t := inf { s > | L s > t } , (2.5)which is F -adapted by the Debut Theorem (Dellacherie and Meyer, 1978) and has continuouspaths if and only if L is strictly increasing. The process H is called the inverse-subordinator of L . Of particular interest for us here is the case where L is an α -stable subordinator, i.e. ψ L ( s ) = s α , α ∈ (0 , time change is a non-decreasing, almost surely finiteprocess ( T t ) t ≥ diverging almost surely to infinity for large times. In particular, both L and H are time changes. If X is an F t -adapted semimartingale, then its time change by T is the F T t -adapted semimartingale X T := ( X T t ) t ≥ . Further, if X is almost surely constant on all sets[ T t − , T t ] we say that X is continuous with respect to T ; in this case many other properties arepreserved, and the semimartingale characteristics of X scale with T .A triangular array of random variables is a collection of random variables ( Y ci , J ci ) i ∈ N ,c> indexed by a scale parameter c such that each ( Y ci ) i ∈ N and ( J ci ) i ∈ N is an i.i.d. sequence, but notnecessarily independent from each other. For fixed c the variable Y ci retains the interpretationof the i -th log-return, and J ci the time elapsed between two consecutive price moves. We cancanonically associate to ( Y ci , J ci ) two families of continuous-time random walks (CTRWs): R ct := [ t ] (cid:88) i =0 Y ci and T ct := [ t ] (cid:88) i =0 J ci , (2.6)and associate to T c the counting process N ct := max { n : T cn ≤ t } . The notation (cid:98) · indicates theFourier transform of probability measures, and the Laplace transform in the time variable isdenoted by L ( · , s ), where s is the new transformed variable. At a microscopic level, we postulate that the time series of returns and trade times, at the timescale c , are determined by a triangular array of random variables ( Y ci , J ci ) i ∈ N , where Y ci ∈ R determines the size of the returns implied by the equity price variation conditional to observinga price revision, and J ci > N c corresponds then to the total number of price movements at t , and the tick-by-tickreturns process Σ c is thus given by time-changing R c with N c :Σ ct := N ct (cid:88) i =0 Y ci . (3.1)At time t the price will have moved by a quantity (cid:80) ni Y ci if the n -th arrival time is recordedbefore t . Or, conditional to n price moves occurred by time s , the price will move again by Y ci t > s if the waiting time variable J cn +1 realizes at a value lesser than t − s . Weassume that there exists a constant risk-free market rate r > c, ∗ t := rt + N ct (cid:88) i =0 Y ci . (3.2)The reasons for this modification shall be explained further on. For the moment, we remarkthat this physical tick-by-tick model must be understood in the sense that only the price inno-vations correspond to market observations. Hence, the linear drift introduced in the randomwalk Σ c, ∗ between two price movements does not give rise to a traded value, and impacts theprice only at revisions time. However, further deterministic trends in the price dynamics, suchas risk premia, are still possible and can be captured by an appropriate choice of Y c . The continuous-time pricing model we describe here is based on a scaling limit of the CTRW Σ c, ∗ for an appropriately selected triangular array ( Y ci , J ci ). This setup encompasses classical math-ematical finance models: when ( Y ci ) i ∈ N are centered with finite variance and J ci = 1 for all i ,then the Central Limit Theorem yields a Brownian motion. If the Y ci have infinite variance andare in the domain of attraction of a stable process X , then their scaling limit yields exactly X .Considering random waiting times for J ic with finite expectation does not improve here the gen-erality of the setting since by the Renewal Theorem N ct ∼ t/ E [ J c ] in probability for large t .Therefore, in order to build processes in which the trade time duration information has impacton the distribution of the scaling limit of Σ c , one has to consider infinite-mean waiting times.Under this choice, taking the limit leads to an anomalous diffusion model for the asset pricedynamics. The following result is central to the entire anomalous diffusions theory. Theorem 3.1 (Becker-Kern, Henry, Jurlewicz, Kern, Meerschaert, Scheffler, Straka) . Assumethat ( Y ci , J ci ) i ∈ N ,c> forms a triangular array of random variables and set R c , T c and Σ c asin (2.6) - (3.1) . If there exists a bivariate L´evy process ( X, L ) , where L is a subordinator withinverse process H as in (2.5) , such that lim c ↑∞ ( R cc , T cc ) = ( X, L ) , (3.3) in the J -topology on the Skorokhod space D ( R × R + ) , then lim c ↑∞ Σ c = (( X − ) H ) + , (3.4) in the J -topology on D ( R ) , where (( X − ) H ) + is the right-continuous version of the processobtained by time-changing by H the left limits process of X . X H that was claimed there to be the limit, the latter can be shown tocoincide with ( X − ) H under such assumptions. This was remarked by Straka and Henry (2011),who also gave a version of the theorem which allows dependence between X and L , but excludesthe possibility of either X or L being a compound Poisson process (CPP). Another proof isobtained by combining (Jurlewicz et al., 2012, Theorem 3.1 and Remark 3.5), which finallyextends the result of Straka and Henry to CPPs. Remark . Unless the J ci are constant or exponentially distributed, the CTRW limit is notMarkovian. Example . For a sequence ( Y i ) i ∈ N of i.i.d. centered random variables with unit variance, let Y ci := c − / Y i ; consider further the i.i.d. sequence ( J ci ) i ∈ N distributed as Exp( λ ), for some λ > c to W λ for some Brownian motion W . Example . Assume that ( Y i ) i ∈ N and ( J i ) i ∈ N are independent sequences of iid random variablesbelonging to the domain of attraction of respectively an α -stable law X with α ∈ (1 , β -stable law L with β ∈ (0 , B n ) n ∈ N and( b n ) n ∈ N , with respective indices − /α and − /β such that B n (cid:80) ni =1 Y i and b n (cid:80) ni =1 J i convergerespectively to X and L almost surely. Then letting Y ci := B ( c ) Y i and J ci := b ( c ) J i , with B ( c ) := B [ c ] and b ( c ) := b [ c ] yields an explicit triangular array, and the theorem to the stable processescanonically associated with X and L . In this case, Theorem (3.1) collapses to (Meerschaert andScheffler, 2004, Theorem 4.2). Example . An explicit representation of the CGMY process as a CTRW limit can been ob-tained by appropriately tempering variables in the domain of attraction of a stable law, asexplained in Chakrabarty and Meerschaert (2011). Combining this with Example 3.2 providesanother explicit CTRW limit representation of (3.4) for a CGMY process X and a stable sub-ordinator L . It is remarkable that the CTRW limit in Theorem 3.1 enjoys a very high degree of analyticaltractability. For example the probability density of an inverse L´evy subordinator H is knownin terms of the L´evy measure of the original process L . Similarly, the law of X H t − can berecovered by integral transforms involving ν X,L and the other Fourier-Laplace characteristics,as explained in Meerschaert and Scheffler (2008) and Jurlewicz et al. (2012). We recall thefollowing from (Jurlewicz et al., 2012, Proposition 4.2):7 roposition 3.2.
Let X H t − be as in the CTRW limit in (3.4) , with law P t . Then L (cid:16) (cid:98) P t ( dz ) , s (cid:17) = 1 s φ L ( s ) ψ X,L ( z, s ) . (3.5)The formula of the Laplace transform of X H t − is particularly simple. Having at hand aspecification for X H t − in terms of the involved characteristic exponents, by virtue of (3.5) weare only one Laplace inversion away from the characteristic function, and we shall see that thisinversion can be computed explicitly in our cases. From a theoretical perspective, the Fourier-Laplace transform of the process provides an interesting connection between the stochasticrepresentation of anomalous diffusions via CTRW limits and the classical characterization oftheir laws as weak solutions of fractional abstract Cauchy problems. For details we refer thereader to Baumer et al. (2005); Meerschaert et al. (2013); Jurlewicz et al. (2012); Meerschaertand Scheffler (2008), and references therein. We introduce here the two anomalous diffusions to establish the connection between tradesduration and the implied volatility surface.
Definition 4.1.
Let X be a L´evy process, L an independent β -stable subordinator, and ( Y ci , J ci ) i ∈ N ,c> a triangular array satisfying (3.3). We define the underlying price S as S t = S exp( rt + Y t ) , S > , (4.1)with Y t := X H t − given by (3.4), and shall consider the following two cases: (SL) The purely subdiffusive L´evy model is such that ( Y ci , J ci ) i ∈ N satisfy the assumptions ofTheorem 3.1 with ( X, L ) in the right-hand side of (3.3); (DRD)
The model with dependent returns and trade durations is such that ( Y ci , J ci ) i ∈ N satisfy theassumptions of Theorem 3.1 with ( X L , L ) in the right-hand side of (3.3).The two models look very similar, the only difference being that the second requires con-vergence of the return innovations to the subordinated L´evy process X L instead of X . Yet,this difference is critical since this subordination is precisely what introduces coupling in theDRD model. We shall denote the CTRW limits Y SL and Y DRD and, correspondingly, theprice processes S SL and S DRD . The underlying standard L´evy model is S = ( S t ) t ≥ with S t = S exp( rt + X t ). Remark . For the SL model, since X is stochastically continuous and independent of H , then X H t − = X H t in law for each t >
0. 8 emark . As β tends to 1, L t tends to t in probability and almost surely. Therefore theusual conditional independence argument shows that S t tends in law to S t . So in the limitingcase, the L´evy models are recovered, and β can be interpreted as a parameter regulating thedivergence from L´evy, and therefore quantifies the degree of ‘anomaly’ of the diffusion. Example . When X is a Brownian motion and L an independent β -stable subordinator theresulting SL model is the subdiffusive Black-Scholes first introduced in Magdziarz (2009). Example . Cartea and Meyer-Brandis (2010) introduce a CTRW model with independenttrade duration and returns, where the conditional waiting time is modelled through a hazardfunction. They in particular consider the latter to be of Mittag-Leffler type P ( T n > t ) = E β ( − t β ) (see also (6.1) below), and the price innovations follow an arbitrary infinitely divisibledistribution. The resulting driving CTRW is a Fractional Poisson process (FPP) as in Laskin(2003); Mainardi et al. (2004) with parameter β . Since an FPP can be represented as a CPP,time-changed by an independent inverse β -stable subordinator (as proved by Meerschaert et al.(2011)), the FPP model by Cartea and Meyer-Brandis (2010) is included in our framework. Example . The original model in Scalas et al. (2000) and Mainardi et al. (2000) also admitsan FPP representation, where the returns innovations follow a stable distribution, and can bewritten in terms of a triangular array limit (Meerschaert and Scalas 2004).
Example . A comprehensive treatment of subdiffusive asset models obtained as fractionalcounterparts of popular L´evy models is provided in Cartea and del-Castillo-Negrete (2007),who tackle the option pricing problem by numerically solving the fractional partial differentialequations characterizing their transition probabilities. In view of the results of (Meerschaertand Scheffler, 2008), all such models admit a time-changed representation of SL type.We recall that a stable subordinator has no drift: therefore the sample paths of Y SL and Y DRD are Lebesgue almost everywhere constant (Bertoin, 1997, Chapter 2), and thusconveniently capture the idea of tick-by-tick trading and persistence of trade duration at alltime scales. This also implicates that all equivalent measures for Y are mutually singular withrespect to the usual diffusion processes. However, the discounted asset value necessarily containsa Lebesgue absolutely continuous part, orthogonal to all equivalent martingale measures for Y ,coming from discounting by the market numeraire (the bank account). Therefore, in order forthe Fundamental Theorem of Option Pricing Delbaen and Schachermayer 1994 to apply, weneed to cancel such part. This clarifies the choice (3.2) of modelling the interest rate effectsexternally to Y . Of course, nothing prevents that the physical dynamics Y itself have a driftin the component X . In Figure 1 we show sample paths of H and Y SL when X is a standardBrownian motion, for two different values of β . As β increases, reversion to respectively thelinear time and a standard Brownian return model with no trades duration effects is observed.The non-Markovian structure of the two processes captures the possible memory effects inprice formation when observing random waiting times between trades. As we shall see later,both the value of the process at time t and the time elapsed since the last price revision influence9he price evolution. Dependence between trade times and price returns is a widely acknowledgedfact, as pointed out in Engle and Russell (1998) and confirmed in several empirical studies. Thismakes the DRD model more realistic compared to the SL one, although the cost/benefit impactin terms of performance of embedding this feature remains to be assessed. For now, observethat the two models have the same number of parameters, so that modelling price/durationdependence does not add any dimension in the calibration and estimation.It would be useful to find a DRD model representation in terms of an independent timechange similar to the one for the SL model. Consider first the special case X = L in Theorem 3.1.Since L − is F t adapted then L H := (( L − ) H ) + (4.2)is an F H t -adapted time change. We can then consider the process X time changed by L H andsee in which relation are the process X L H and the limiting process (( X L ) − H ) + coming from Theo-rem 3.1, where we recall that X L is just a L´evy subordinated process. Although these processesbear some resemblance their construction is different: in particular now X is independent of thetime change. However, ( X − L ) H = X ( L − ) H (4.3)because X has left limits and L is increasing. But then(( X − L ) H ) + = ( X ( L ) − H ) + = X (( L − ) H ) + = X L H (4.4)since X is right-continuous. There are then two ways of looking at the DRD returns process.The CTRW limit definition gives us a dependent representation using a continuous time change.The equalities above give instead an independent representation employing a discontinuous timechange. Both will be useful in the sequel.Finally, let us briefly comment on the nature of the process L H . It is easy to show that, forany t ≥ L H t − = sup { s < t : s = L u , for some u ≥ } . (4.5)In light of this identification, the process ( L H t − ) t ≥ is sometimes called the last passage process (Bertoin, 1997, Chapter 1) and plays an important role in potential theory for L´evy processes:it can be seen as the discontinuous increasing process which represents the past time at which H has started resting at its current (time- t ) location. Now the times of discontinuity of ( L H t − ) t ≥ coincide with the points in the image of L isolated on their right, which on account of L beingdriftless, by (Bertoin, 1997, Chapter 1, Proposition 1.9) is a set of Lebesgue measure zero. Fromthis we conclude that L H is a right-continuous modification of the first passage time. Moreoverthe post-jump value of L H is exactly t , and in any case L Ht ≤ t almost surely. This ties in withthe interpretation of L H as a delayed calendar time.Notably, the distribution of L Ht is known. We have the following result:10 roposition 4.2. Denote by B a,b the Beta distribution with parameters a and b . For any t > , L Ht is distributed as t B β, − β .Proof. See (Becker-Kern et al., 2004, Example 5.5) or (Jurlewicz et al., 2012, Example 5.2).This underpins the greater analytic tractability of the DRD model with respect to the SLmodel: somewhat paradoxically, the more realistic model is also the more explicit. Proposi-tion 4.2 clarifies how the DRD model captures the paradigm of Engle (2000) and Dufour andEngle (2000). The DRD time-changed evolution obeys a form of delayed calendar time whosemass in [0 , t ] concentrates more around 0 or t depending on whether β is close to zero or one(Figure 3). This mass represents the quantity of delay one has to apply to X to obtain thecurrent price value. When L has a low β , that is when duration of trade is higher, the price evo-lution is stickier, since t B β, − β is much smaller than t with high probability. This is associatedwith a reduced impact of the individual trades on the price process because the informationalcontent of sporadic trading is low. Conversely, as t B β, − β is close to t with high probability(namely when β is close to one) we observe a higher trading activity, typically associated withthe presence of informed traders. In such a case the contribution of each single trade to theprocess of price formation is greater, and the impact of trading on price higher. A similar rea-soning applies to the SL model. Here combining subordination with independence ‘delays’ theevolution of X for the time necessary to the next price revision to happen, but the resultingmove retains the variance of an earlier point-in-time position of the process X . Therefore, again,the lower the β , the stickier the price dynamics. We derive some statistical properties of the SL and DRD models and provide some initial insighton the structure of the volatility surface they generate, anticipating the full analysis in Section 7.We begin with the moments of the DRD model, whose analytic tractability plays a major role.The following proposition extends (Leonenko et al., 2014, Theorem 2.1) to higher cumulants.In this section, X is a given L´evy process, T an independent time change, and we let κ i and τ i denote their respective i -th cumulants, which we assume to exist for i = 1 , . . . , Proposition 5.1.
The process Y := X T has moments up to order four, and its cumulants read κ Y = τ κ , κ Y = τ κ + κ τ ,κ Y = τ κ + 3 κ κ τ + κ τ , κ Y = (3 κ + 4 κ κ ) τ + 6 κ κ τ + κ τ + κ τ . (5.1) Proof.
In our notation κ n = − (cid:16) i n ψ ( n ) X (0) (cid:17) . We proceed as in (Leonenko et al., 2014, Theorem2.1), where the usual conditioning argument yields E [ Y t ] = i ddz E (cid:2) e − i zY t (cid:3) (cid:12)(cid:12)(cid:12) z =0 = i ddz E (cid:104) e − ψ X ( z ) T t (cid:105) (cid:12)(cid:12)(cid:12) z =0 = − i ψ (cid:48) X (0) E [ T t ] , (5.2)11hich gives κ Y . Next E [ Y t ] = − d d z E (cid:2) e − i zY t (cid:3) (cid:12)(cid:12)(cid:12) z =0 = ψ (cid:48)(cid:48) X (0) E [ T t ] − ψ (cid:48) X (0) E (cid:2) T t (cid:3) , (5.3)Subtracting from (5.3) the square of (5.2) reconstructs τ and yields κ Y . Similarly, E [ Y t ] = − i d d z E (cid:2) e − i zY t (cid:3) (cid:12)(cid:12)(cid:12) z =0 = − ψ (cid:48)(cid:48)(cid:48) X (0) E [ T t ] + 3 E (cid:2) T t (cid:3) ψ (cid:48) X (0) ψ (cid:48)(cid:48) X (0) + i ψ (cid:48) X (0) E (cid:2) T t (cid:3) ; (5.4)calculating E [ Y t ] − E [ Y t ] E [ Y t ] + 2 E [ Y t ] and factoring the τ i as necessary we obtain κ Y . Thelast term κ Y is obtained analogously.The above proposition confirms the well-known fact that a L´evy model X subordinated by aL´evy process L creates non-zero skewness and kurtosis even in the presence of a mesokurtic andsymmetric parent process X such as a Brownian motion. Our situation here is identical, andcarries the message that trade duration alone can be a determinant of departure from normalityof returns (thus, in an option pricing perspective, creating volatility smile). However, the termstructure analysis of the moments is completely different. The key fact is that the moment timedispersion of a time-changed L´evy process only depends on the moments of the time change, andnot on the moments of X . In the usual L´evy subordination case, that is when T is a L´evy process,one then sees that the moments are linear in t , consistently with the fact that the subordinatedprocess is itself L´evy. As a consequence the skewness and kurtosis of the returns vanish withtime. In contrast, our framework produces a nonlinear time evolution of the moments, whichwe analyse in detail for the DRD model, where such evolution is polynomial. Proposition 5.2.
For any t ≥ , the first four cumulants of Y DRDt are κ Y = βκ t,κ Y = βκ t + κ − β ) βt ,κ Y = βκ t + 3 κ κ − β ) βt − κ − β ) β (2 β − t ,κ Y = βκ t + 4 κ κ + 3 κ β (1 − β ) t − − β ) β (2 β − κ κ t + κ (1 − β ) β (2 − β (1 − β )) t , (5.5) and the following asymptotic relations hold: lim t ↑∞ Skew( Y t ) = 2 √
23 1 − β (cid:112) (1 − β ) β sgn( κ ) , lim t ↑∞ Kurt( Y t ) = 1 β (1 − β ) − , lim t ↓ √ t Skew( Y t ) = κ (cid:112) βκ , lim t ↓ t Kurt( Y t ) = κ βκ . (5.6) Proof.
By explicitly integrating the Beta probability density function we have the central mo-12ents of T t : µ T = E (cid:2) L Ht (cid:3) = βt = τ , (5.7) µ T = V (cid:2) L Ht (cid:3) = 12 (1 − β ) βt = τ , (5.8) µ T = E (cid:2) ( L Ht − τ ) (cid:3) = −
13 (1 − β ) β (2 β − t = τ , (5.9) µ T = E (cid:2) ( L Ht − τ ) (cid:3) = β − β ) (2 − − β ) β )) t = τ + 3 τ . (5.10)Since in the DRD model X t and L Ht are independent, we can solve the above equations for τ i andsubstitute in (5.1) obtaining (5.5). Calculating further the normalized cumulants Skew( Y t ) = κ Y / ( κ Y ) / and Kurt( Y t ) = κ Y / ( κ Y ) and taking respectively the limits for large t and theleading order around t = 0 imply the limits in the proposition.In the DRD model, as the time scale gets larger, higher moments do not vanish, but convergeto a level that only depends on β , and not on the value of the L´evy cumulants (the sign of κ dictates the sign of the skewness). As frequently noted, leptokurtosis and negative skewness ofreturns are important drivers of implied volatility smiles. It thus makes sense to deduce thatnon-zero time limits of skewness and excess kurtosis determine persistence of the volatility smileover time. In contrast, for t close to zero, moment explosions are observed, as in the L´evycase; the rate of this explosion is exactly that of exponential L´evy models, including–up to anormalization by β –the constant factor. This suggests that the short-term smile/skew behaviourof the DRD implied volatility should be identical to that of the underlying L´evy model. We willverify these intuitions and make the matters more precise in Section 7.The analysis of the returns series properties stems from the observation that the models weare studying, although not Markovian with respect to their own filtration, admit a Markovianembedding. Remarkably, the marginal distributions of this embedding are known for the DRDprocess. For any t ≥
0, we define the backward renewal time V t := t − L Ht , (5.11)which represents the time elapsed from the current instant t to the previous price move. Knowingthe price at t and the time since the last price move is enough to fully describe the law of thefuture asset evolution. Proposition 5.3.
The following properties hold:(i) the pairs ( Y SL , V ) and ( Y DRD , V ) are time-homogeneous Markov processes;(ii) the process Y SL has correlated increments, whereas Y DRD has uncorrelated increments;(iii) the increments of Y SL are non-stationary, whereas the increments of Y DRD are weaklystationary; roof. Item (i) is proved in (Meerschaert and Straka, 2014, Theorem 4.1). For the SL model,statement (ii) can be deduced from (Leonenko et al., 2014, Example 3.2, Equation 9), since inour case E [ X ] (cid:54) = 0. In the case of the DRD model, for s ≤ t , we can write (we drop the modelsuperscript for convenience) E [ X t X s ] = E [( X t − X s ) X s ] + E (cid:2) X s (cid:3) = ( t − s ) s E [ X ] + s V [ X ] + s E [ X ] = ts E [ X ] + s V [ X ] , so that by independence and conditioningCov( Y t , Y s ) = E (cid:2) L Ht L Hs (cid:3) E [ X ] + E (cid:2) L Hs (cid:3) V [ X ] − E (cid:2) L Ht (cid:3) E (cid:2) L Hs (cid:3) E [ X ] = Cov (cid:0) L Ht , L Hs (cid:1) E [ X ] + E (cid:2) L Hs (cid:3) V [ X ] . (5.12)Thus, considering increments and using the above, together with Proposition 5.1,Cov( Y t − Y s , Y s ) = Cov( Y t , Y s ) − V [ Y s ] = E [ X ] (cid:0) Cov( L Ht , L Hs ) − V [ L Hs ] (cid:1) = E [ X ] Cov( L Ht − L Hs , L Hs ) , (5.13)so absence of returns autocorrelation is equivalently checked on L Ht . Now (Meerschaert andStraka, 2014, Example 5.4) give the conditional transition probabilities p t ( y , v , dy, dv ) := P ( L Ht ∈ dy, V t ∈ dv | y , v ) of the Markov process ( Y t , V t ) as: p t ( y , , dy, dv ) = v − β Γ(1 − β ) ( t − v ) β − Γ( β ) δ y + t − v ( dy ) dv { Let S be of SL or DRD type under P , and Q ∼ P an equivalent measuresuch that ( e − rt S t ) t ≥ is a L´evy exponential martingale under Q . Define the Radon-Nikodymderivative Z := d Q /d P , consider its time change Z H , for any H of the form (2.5) , and introduce he measure (cid:101) Q via d (cid:101) Q /d P := Z H . Then (cid:0) e − rt S DRDt (cid:1) t ≥ and (cid:0) e − rt S SLt (cid:1) t ≥ are martingalesrespectively under Q and (cid:101) Q .Proof. For the S SL models the statement follows from (Torricelli, 2020, Theorem 2) by taking L t to be a β -stable subordinator, S t = e − rt S SLt , X t = Z t and H t = 1, which corresponds to nochange of measure in the subordinator. For the DRD model it suffices to observe that L H is abounded family of stopping times and thus e − rt S DRD = exp (cid:0) X ∗ L H (cid:1) is a martingale under Q byDoob’s Optimal Sampling Theorem.Again we emphasize that this is a subset of all the possible equivalent martingale measuresand that for technical reasons we ignore a market price of duration risk. A model in whichthis risk can be priced can be obtained for example by considering for L the wider class oftempered stable subordinators, which is closed under the Esscher transform. This class, alongwith related questions of market completeness, is studied in Torricelli (2020); see also Fries andTorricelli (2020) for the situation when trade duration is caused by market suspensions. Having established that the risk-neutral specification comes in the form of a time-changed mar-tingale exponential, Proposition 3.5 can be combined with standard integral price representa-tions to yield semi-closed-form valuation formulae. Remarkably, the characteristic functions ofthe log-price in the SL and DRD models admit a very simple representation in terms of theone-parameter Mittag-Leffler function E a ( z ) := ∞ (cid:88) k =0 z k Γ( ak + 1) , (6.1)where Γ is the usual Gamma function, and of the confluent Hypergeometric function F ( a, b ; z ) := ∞ (cid:88) k =0 ( a ) k ( b ) k z k k ! . (6.2) Theorem 6.2. Let Y be either process in Definition 4.1, and F ( · ) a contingent claim on S maturing at T . Assume that x (cid:55)→ f ( x ) := F ( e x ) is Fourier-integrable and let S f be the domainof holomorphy of its Fourier transform (cid:98) f . Let Φ t ( z ) := E [ e − i zY t ] be the characteristic functionof Y t taken according to the relevant measure as described in Proposition 6.1, denote by S Y itsholomorphy domain, and assume S f ∩ S Y (cid:54) = ∅ . The price P of the derivative paying F ( S T ) attime T is given by P = E (cid:2) e − rT F ( S T ) (cid:3) = e − rT π (cid:90) i γ + ∞ i γ −∞ Φ T ( z ) (cid:98) f ( z )( S e rT ) i z dz. (6.3) The value γ ∈ R is chosen such that the integration line lies in S f ∩ S Y and Φ t ( z ) = (cid:40) E β (cid:16) − ψ X ( z ) t β (cid:17) , if Y = Y SL , F ( β, , − tψ X ( z )) , if Y = Y DRD . (6.4)16 roof. Under the given assumptions, the Plancherel representation (6.3) is standard (see Lewis(2001) for example), and we only need to prove (6.4). In the SL model, by independence of X and L we have ψ ( s, z ) = φ L ( s ) + ψ X ( z ) = s β + ψ X ( z ), and Proposition 3.2 then yields L (Φ t ( z ) , s ) = s β − s β + ψ X ( z ) . (6.5)Inverting the right hand-side, as in (Haubold et al., 2011), one obtains (6.4). In the DRD modelafter conditioning and applying Proposition 4.2, we obtainΦ t ( z ) = E (cid:2) exp (cid:0) − ψ X ( z ) L Ht (cid:1)(cid:3) = E [exp ( − tψ X ( z ) B β, − β )] , (6.6)and the statement follows from the characteristic function of B β, − β . Remark . Fast computational routines for the Mittag-Leffler and the confluent hypergeomet-ric functions are available in most software packages. Also, the two functions can be unified ina single software implementation by observing that the three-parameter Mittag-Leffler function E a,b,c ( z ) = ∞ (cid:88) k =0 ( c ) k z k Γ( ak + b ) (6.7)is such that E a, , ( z ) = E a ( z ) and E , ,c ( z ) = F ( c, , z ). Furthermore if a = b = c = 1,then (6.7) reverts to the standard exponential, which is consistent with the fact that S SL and S DRD revert to the exponential L´evy model S . Remark . The function E β is entire and F ( β, , − tψ X ( · )) is regular in the complex planewithout the negative real axis; hence S f ∩ S Y (cid:54) = ∅ depends on the domain of ψ X and (cid:98) f only. Remark . If X H has an FPP structure, then (6.3) coincides with the formula given by (Carteaand Meyer-Brandis, 2010, Theorem 3), when the jump sizes have infinitely divisible distribution.One sees that the pricing formulae are formally obtained from the standard L´evy case byreplacing the exponential function with two different kinds of ‘stretched exponentials’. Theparameter β relaxes the shape of the characteristic function, in particular in the tails, therebygenerating large-maturity prices very different from the base case. This overcomes the ‘curseof exponentiality’ of the standard models (both L´evy and exponentially-affine), for which thelong-maturity option prices follow Laplace-type asymptotics of leading order exp( − T ) / √ T . Wewill detail this better, together with its implications on the volatility surface, in Section 7 below.Note that the two functions (6.1) and (6.2) have very different behaviours. In Figure 2, wecan see for example that (6.1) has a cross-over region where its decay transitions from super tosub-exponential, whereas in (6.2), the integrand always dominates the exponential. This has aclear impact on the shape of the volatility surface, as illustrated numerically in Section 8.17 Time asymptotics of the volatility surface Bearing in mind the discussion so far, we naturally expect implications of trade duration (at leastin the form we chose to model it) on the volatility surface. The anomalous diffusions processes weconstructed are subdiffusions, and as such have a slower distributional dispersion rate than thebenchmark L´evy models, hence a slower option price convergence for large maturity. That said,since Black-Scholes is a L´evy model, inversion of the Black-Scholes formula using subdiffusiveoption prices should generate a vanishing implied volatility term structure in order to match theslower price time evolution.Less intuitive is to find a reason why the long-term skew should decline slower than standardL´evy and stochastic volatility models. A first answer is provided by Section 5: skewness andkurtosis in our models do not tend to zero as time grows but converge to some strictly positivelevel. Therefore Gaussian temporal returns aggregation is precluded, and time reversion to aflat volatility might be pushed further away in time . However, as we shall show, an exhaustiveanswer is provided by the fact that skew and level of the implied volatility are connected, andthe property of a vanishing asymptotic implied volatility is sufficient to hamper the skew timedecay.In this section we generically indicate with Q any of the two risk-neutral measures of Propo-sition 6.1. Without loss of generality, we assume here r = 0 and S = 1 and, denote C ( K, T )the Call option price with strike K and maturity T . In the Black-Scholes model dS t = σS t dW t ,with σ > 0, the price of such a Call option is given by C BS ( K, T, σ ) = S N (cid:16) d (cid:16) σ √ T (cid:17)(cid:17) − K N (cid:16) d (cid:16) σ √ T (cid:17) − σ √ T (cid:17) , (7.1)where d ( z ) := − log( K ) z + z , N denotes the standard Gaussian cumulative distribution function,and n its derivative, the Gaussian density function. For K, T ≥ 0, the implied volatility σ ( K, T )is the unique non-negative solution to C ( K, T ) = C BS ( K, T, σ ( K, T )), and the implied volatilityskew is defined as S ( K, T ) := ∂σ∂K ( K, T ) . (7.2)It is known by Rogers and Tehranchi (2010) that S ( K, · ) converges to zero as the maturityincreases, for each K . We begin with the following model-free lemma which, under some mildassumptions on the underlying distribution, connects the time decay of the skew with its level. Lemma 7.1. Let ( S t ) t ≥ be a martingale such that the law of S t is absolutely continuous foreach t and converges to zero in distribution as t tends to infinity.(i) For any K ≥ , if lim T ↑∞ √ T σ ( K, T ) = ∞ then, as T tends to infinity, S ( K, T ) = 2 T σ ( K, T ) (cid:20) K ) − T σ ( K, T ) + O (cid:18) T − σ ( K, T ) (cid:19)(cid:21) − Q ( S T ≥ K ) √ T n ( d ( σ ( K, T ) √ T )) ; (7.3) Gaussian aggregation is by no means responsible of the smile flattening, as shown by Rogers and Tehranchi(2010). ii) as T tends to zero, S (1 , T ) = (cid:114) πT (cid:34) − Q ( S t ≥ − σ (1 , T ) √ T √ π + O (cid:0) σ (1 , T ) T (cid:1)(cid:35) . (7.4) Proof. We only prove the first statement, as the second one is proved in (Gerhold et al., 2016,Lemma 2). Since S t has an absolutely continuous law, then by (Figueroa-L´opez et al., 2011,Lemma C.1), S in (7.2) exists, ∂ K C ( K, T ) = − Q ( S T ≥ K ), and the chain rule yields S ( K, T ) = − ∂ K C BS ( K, T, σ ( K, T )) + Q ( S T ≥ K ) ∂ σ C BS ( K, T, σ ( K, T )) . (7.5)Set z = √ T σ ( K, T ). Using the formulae for the Black-Scholes Delta and Vega: S ( K, T ) = N ( − d ( z )) − Q ( S T ≥ K ) √ T n ( d ( z )) . (7.6)where we recall that N ( · ) is the standard Gaussian cumulative distribution function. Since, as x tends to infinity N ( − x ) = n ( x ) x (cid:18) − x + O (cid:0) x − (cid:1)(cid:19) and 1 d ( x ) = 2 x + 4 log( K ) x + O (cid:0) x − (cid:1) , (7.7)then N ( − d ( z )) √ T n ( d ( z )) = 1 √ T d ( z ) (cid:18) − d ( z ) + O (cid:0) d ( z ) − (cid:1)(cid:19) = 2 √ T z (cid:18) K ) − z + O (cid:0) z − (cid:1)(cid:19) , (7.8)and (7.3) follows by substituting z and combining the above with (7.6). Remark . If S = exp( X ) is a martingale for some L´evy process X , from the proof ofProposition 6.1, our models can be written as S T t for some time change T t , so that S T t convergesto zero almost surely as t tends to infinity, provided we know this to hold for S t . Such a propertyfor exponential L´evy models can be proved using fluctuations identities, since the assumption E [ X ] < X t diverges to −∞ . Because of Jenseninequality, a negative first moment is always the case for X t when S is a martingale. Regardingthe absolute continuity of the price process, this follows from the fact that the law of the involvedprocesses are weak solutions of fractional Cauchy problems. These can be found using argumentsanalogous to (Jurlewicz et al., 2012, Examples 5.2-5.4).Part (i) of this lemma implies that the level and skew of the implied volatility are entangled:one cannot modify the leading order 1 /T of the skew decrease without postulating a zero ordiverging asymptotic implied volatility level. In turn, a declining implied volatility can only beattained through a convergence rate of option prices distributions to the spot price slower thanGaussian, which is precisely the distinguishing feature of anomalous diffusions-based models.Part (ii) is an already known fact, originally observed in (Gerhold et al., 2016, Lemma 2), whichhighlights a very stringent relationship between the prices of Digital options and the small-timeat-the-money skew. It will be used later in Corollary 7.5.19 heorem 7.2. As T tends to infinity, we have the following asymptotic expansions for the Callprice C ( K, T ) , for any K ≥ :(i) in the DRD model with β ∈ (0 , , with the interpretation that L t = t when β = 1 , thereexist C β and c β > such that C ( K, T ) = 1 − { β (cid:54) =1 } C β Γ(1 − β ) 1 T β (cid:20) O (cid:18) T (cid:19)(cid:21) − c β Γ( β ) e − T ψ X (i / T / − β (cid:20) O (cid:18) T (cid:19)(cid:21) ; (7.9) (ii) in the SL model with β ∈ (0 , , there exists C β > such that C ( K, T ) = 1 − C β Γ(1 − β ) 1 T β (cid:18) O (cid:18) T (cid:19)(cid:19) . (7.10) Proof. Since we are under the assumptions of Theorem 6.2, we can consider the price represen-tation for a Call option C ( K, T ) = 1 − π (cid:90) ∞−∞ e ( i u + ) log( K ) u + 1 / T (cid:18) u + i2 (cid:19) du (7.11)which can be obtained from (6.3) by moving the integration contour inside the strip (cid:61) ( z ) = 1 / T tends toinfinity of C ( K, T ) under the integral sign. So once we determined the asymptotic expansionof Φ T we can integrate the resulting expression to get the asymptotic equation of interest.Assume β < does not occur. The asymptoticexpansion of F ( a, b, z ), for large | z | is (Luke, 2012, Chapter 4): F ( a, b, z ) ∼ Γ( b )Γ( b − a ) z − a e i δπa F (cid:18) a, a − b, − z (cid:19) + Γ( b )Γ( a ) z a − b e z F (cid:18) b − a, − a, z (cid:19) (7.12)with δ = 1 if (cid:61) ( z ) > δ = − (cid:61) ( z ) = 1 / 2, since (cid:60) ( ψ X ( z )) > 0, inequation (6.4) for large T | ψ X ( z ) | we have the well-defined behaviour F ( β, , − T ψ X ( z )) ∼ ( T ψ X ( z )) − β Γ(1 − β ) F (cid:0) β, β, ( T ψ X ( z )) − (cid:1) + e − T ψ X ( z ) ( − T ψ X ( z )) β − Γ( β ) F (cid:0) − β, − β, − ( T ψ X ( z )) − (cid:1) ∼ ( T ψ X ( z )) − β Γ(1 − β ) + e − T ψ X ( z ) ( − T ψ X ( z )) β − Γ( β ) (7.13) The asymptotic behaviour of complex-valued functions can be different in different regions of the complexplane, which is normally referred to as Stokes phenomenon. A complex-valued function has a limit along adirection if it eventually takes values in only one of such areas. x → F ( a, b ; x ) = 1, for all a, b . In order to substitutein (7.11) the above expression for large T , we need a uniformity argument in u . Notice first thatas | z | tends to infinity, | ψ X ( z ) | also tends to infinity because the risk-neutral drift of X mustbe nonzero by (Bertoin, 1997, Corollary 1.1.3). This implicates that along any line (cid:61) ( z ) = c , | ψ X ( z ) | is strictly increasing. Also ψ X is even in its real part and odd in its imaginary part, sothat | ψ X ( · + i / | on such sets must be an even function. We conclude that | ψ X ( · + i / | , has apositive minimum at the origin. Therefore so long as T is much larger than 1 / | ψ X (i / | we canreplace Φ T in (7.11) with (7.13). Integrating the first term of the resulting expression producesthe first term in (7.9) with C β = 12 π (cid:90) ∞−∞ e (i u +1 / 2) log( K ) ( u + 1 / ψ X ( u + i / β du. (7.14)Regarding the exponential sub-leading terms we have to analyse I β ( T ) := (cid:90) ∞−∞ e (i u +1 / 2) log( K ) e − T ψ X ( u +i / ( u + 1 / ψ X ( u + i / − β du, (7.15)which can be treated using the saddle point method as in Andersen and Lipton (2012). Fromthe previous discussion, ψ X ( · + i / 2) has a stationary point in 0 and further by ψ (cid:48)(cid:48) X (i / > T I β ( T ) ∼ √ π √ Kψ X (i / − β (cid:112) ψ (cid:48)(cid:48) X (i / T , (7.16)which yields the second term in (7.9) with c β = 4 √ Kψ X (i / − β (cid:112) πψ (cid:48)(cid:48) X (i / . (7.17)When β = 1 the whole proof collapses to the well-known steepest descent argument (Ander-sen and Lipton, 2012, Section 7) for the L´evy models price representation integral.In the SL model we have, for any given β < 1, that so long as πβ/ < θ < min { π, πβ } theasymptotic series for E β is given by (Haubold et al., 2011, Equation 6.5) E β ( z ) = e z /β β n − (cid:88) k =1 − βk ) 1 z k + O (cid:0) z − n (cid:1) , for | arg( z ) | < θ, n − (cid:88) k =1 − βk ) 1 z k + O (cid:0) z − n (cid:1) , for θ < | arg( z ) | ≤ π. (7.18)Since (cid:60) (Ψ X ( u + i / > 0, for all α in the line (cid:61) ( z ) = 1 / T big enough such that πβ < | arg( − ψ X ( u + i / T β ) | , so that for T > T the Stokes lines are not crossed. The correctexpression is thus the second line in (7.18), and we can repeat what argued in the DRD case.21 emark . The second term in (7.9) is clearly negligible for large T compared to the leadingorder, when β is smaller than 1. However for fixed T , as β approaches one its contributioncannot be neglected. This term has been included to clarify the convergence to the L´evy model.Such correction is not present in the SL model, and as β = 1 the price approximation simplybreaks down (however, by dominated convergence we still have convergence of prices).Proposition 7.2 clarifies the aforementioned slower convergence of Call prices compared toL´evy (or exponentially-affine stochastic volatility) models. As already remarked, it can bethought of as a direct consequence of the slow, subdiffusive time spread of the asset returns.More specifically, the nature of the distribution implies that the pricing integral does not obeythe Laplace decay rate, since the integrand is not of the form exp( − T f ( x )) g ( x ). One insteadobtains a vanishing long-term volatility, and hence by Lemma 7.1 a persistent long-term skew,as we illustrate below: Corollary 7.3. For β ∈ (0 , , the leading-order asymptotic for large T of the implied volatilityin both the DRD and SL model satisfies σ β ( K, T ) ∼ (cid:115) T W (cid:18) K T β Γ(1 − β ) πC β (cid:19) , (7.19) where W is the Lambert function and C β > . Furthermore, for all K , α > / , lim T →∞ T − α S β ( K, T ) = 0 and lim T →∞ S β ( K, T ) √ T = 0 . (7.20) Proof. Using k = log K , the first-order expansion of the Black-Scholes price is simply C BS ( K, T, σ ) = 1 − √ K exp (cid:16) − σ T (cid:17) σ √ πT (cid:18) O (cid:18) T (cid:19)(cid:19) . (7.21)By definition of implied volatility, to determine the leading order of σ in the DRD (respec-tively SL) model, we need to equate the above respectively to the Call option price expan-sions (7.9) (resp. (7.10)) and then solve for σ . We haveexp (cid:18) − σ T (cid:19) √ Kσ √ T π = C β Γ(1 − β ) T − β , (7.22)with C β = C βi , i = 1 , z = σ T / M = √ K Γ(1 − β ) / ( C β √ π ), w = M T β , then the equality (7.22) reads e z z = w . Since w > z can be performed along the real axis so that W is well-defined, and (7.19)follows. Since W ( T ) ∼ log( T ) as T tends to infinity, then σ β ( K, T ) ∼ (cid:114) log( M T β ) T , (7.23)therefore T α σ β ( K, T ) converges to zero for all α < / 2, which means that the first term of (7.3)tends to zero slower than T − α , for all α > / 2, but faster than T − α .22tudying the asymptotics of the last term in (7.3), similar arguments to those of Propo-sition 7.2 imply that the long-term price decay for the Digital option I { S T ≥ K } is identical tothat of the Call option, namely c/T β for some c > 0. Then substituting (7.23) together with d ( x ) ∼ x/ 2, in the second term of (7.9), the proof follows from the asymptotic equivalence Q ( S T ≥ K ) √ T n ( d ( √ T σ β ( K, T ))) ∼ c exp (cid:16) T σ β ( K,T ) (cid:17) T β +1 / = c exp (cid:16) log( M T β )2 (cid:17) T β +1 / = cM √ T . (7.24)In light of the corollary above, persistence of the skew is to be interpreted as follows: the skewdeclines slower than any power of T − bigger than 1 / /T ) butalways faster than T − / . It is then natural to ask if these structural differences in the impliedvolatility of anomalous diffusion models also affect the small-maturity limit. This turns out notto be the case, at least for the DRD model, and the underlying L´evy model asymptotics aremaintained. More precisely, we have the following for Digital option prices: Proposition 7.4. If the underlying L´evy process X is such that Q ( S t ≥ 1) = c + c ε t ε + o ( t ε ) , (7.25) for some c , c ε , as t tends to zero, with < ε ≤ , then, with c β,ε := Γ( β + ε )Γ( β )Γ(1+ ε ) , Q (cid:0) S DRDt ≥ (cid:1) = c + c β,ε c ε t ε + o ( t ε ) . (7.26) Proof. Proposition 4.2 allows us to write Q (cid:0) Y DRDt ≥ (cid:1) = (cid:90) t Q ( X s ≥ s β − ( t − s ) − β Γ( β )Γ(1 − β ) ds. (7.27)Now, notice that E (cid:104) B aβ, − β (cid:105) = Γ( β + a )Γ( β )Γ(1+ a ) for all a > 0, and that for t sufficiently small, Q ( X t ≥ 0) = c + c ε t ε + f ( t ) , (7.28)where f ( t ) = o ( t ε ) is a bounded function in a neighbourhood of the origin. The zero- andfirst-order terms of (7.26) are then clear, and by dominated convergence,lim t → E [ f ( t B β, − β )] t − ε = (cid:90) lim t → f ( ts ) t ε s β − (1 − s ) − β Γ( β )Γ(1 − β ) ds = 0 , (7.29)which yields the small-o order ε of the remainder.By combining this result with Lemma 7.1 we have the desired analogy of the short-termATM skew of the DRD and L´evy volatility asset pricing models. Corollary 7.5. Under the assumptions of Proposition 7.4 the DRD model and its underlyingL´evy model have the same short-term ATM skew rate. roof. Substituting respectively (7.25) and (7.26) in (7.4) yields the claim.As extensively discussed in Gerhold et al. (2016), Equation (7.25) essentially encompasses allthe popular L´evy models and features very different behaviours: for example c = 1 if the processhas finite variation, whereas c = 1 / c ε = 1 / d/ ( σ √ π ), ε = 1 / σ and risk-neutral drift d (Gerhold et al., 2016, Theorem 1 and Lemma 2). Asit emerges, the critical value for which higher order terms are needed is c = . In the DRDmodel, introducing c β,ε does not change the asymptotic analysis, as c remains the same. Corollary 7.6. If X satisfies (7.25) , then the DRD model and the underlying exponential L´evymodel S have the same short-maturity at-the-money skews. In the next section we bring together all these results and see how they lead to modelcalibration improvements when a persistent implied volatility skew is observed. We visualize the volatility surfaces extracted from the DRD and SL models in Figures 4 to 7.For X , we use a Brownian motion (Figures 4 and 5) and a CGMY process with parameterstaken from Carr. et al. (2001) (Figures 6 and 7), and consider moneynesses ± 40% ATM andmaturities up to two years. In each figure, the smile of the anomalous diffusion is compared tothat of its underlying L´evy model S .First and foremost the slow decay of the volatility skew in anomalous diffusion models pre-dicted by Proposition 7.2 and Corollary 7.3 is apparent in all cases. Even though our resultonly predict an asymptotic rate of skew vanishing, our numerical tests, at least in the DRDcase, indicate that such a slower rate manifests itself already very early on. More research isnecessary to see whether and how Proposition 7.2 can be improved.In Figures 4 and 5 the volatility smile and skew of the anomalous diffusion model are presenteven if the L´evy generating returns process is a Brownian motion. In other words, this confirmsthat introducing infinite-mean trade durations in a standard CTRW is alone sufficient to generatea smile, consistently with Proposition 5.2. The smile appears rather symmetric, in line with theintuition that trade duration should have little skew impact, as it does not influence out-of-the-money prices any differently than in-the-money ones. This already suggests some orthogonalitybetween β and the L´evy parameters. In the Brownian motion case, β is thus ‘overloaded’, beingresponsible for both the smile convexity and its decay rate. This is relaxed in 6 and Figures 7by endowing X with a proper L´evy structure (CGMY); there a short-term skew arises while theskew term structure maintains its slower flattening rate, dictated by β .In Figures 4 and 6 we observe the repercussion on the implied volatilities of the ‘cross-over’phenomenon (Figure 2) generated by the Mittag-Leffler and exponential types of the character-24stics functions of the SL and pure L´evy models. The level of the SL surfaces transitions from ashort-term regime where the implied volatilities are higher to a long-term one in which they arelower than those of the underlying L´evy models (eventually tending to zero). Such a transitionseems to be very sharp.The time sections from Figures 6 and 7 are shown in Figures 8 and 9 and further highlightthe remarks above. Figures 10 to 13 highlight the convergence of the time sections to those ofthe underlying L´evy model as β approaches one. For the SL model this convergence is fromabove, while it is from below for the DRD model. Note also that the DRD model exhibits asharper ATM skew than the SL model. Corollary 7.3 and Proposition 7.4 suggest that, from a calibration viewpoint, the models shouldbehave as follows: the L´evy parameters have a short-time scale effect, unaffected when intro-ducing β , and they should hence absorb the short-time skew and smile. However, β is thevery component governing the long-term structure of the surface, where the L´evy structure isflat and has no impact, and should thus allow to pick up the long-term skew. To test this wegenerate 3-month and 6-month volatility skews from a given L´evy model S , which representour baseline synthetic market data. In order to generate two scenarios of persistent volatilityskew, while keeping the 3-month fixed, we shift the 6-month skew forward to make it coinciderespectively with the 1-year and 18-month skews. We then cross-sectionally calibrate S , S SL and S DRD to the 3-month and 6-month skews in the baseline scenario, and the 3-month and1-year (respectively 18-month) sections in the first and second scenarios.The calibration problem has been set-up as follows. Let C ( K, T ; β, Γ ) be the theoreticalCall prices from the SL or DRD models, where Γ denotes the L´evy parameters for X , and by C ( K, T ) the synthetic market prices obtained with the procedure described above. The goalis to minimize the root mean squared error (RMSE) on a set of quoted prices of n maturitiesand m strikes. The calibration parameters are then:arg min β, Γ (cid:118)(cid:117)(cid:117)(cid:116) m × n m,n (cid:88) i,j =1 | C ( K i , T j ; β, Γ ) − C ( K i , T j ) | . (8.1)We used m = 8, corresponding to ± 20% moneyness from S = 100 and a step size of ∆ K = 5.Further n = 2 with T j chosen according to the scenarios illustrated above explained above.To solve the optimization problem above we use the Differential Evolution global minimiza-tion algorithm. Its core MATLAB implementatio, from Gilli and Schumann (2011), is freelyavailable at the GitLab repository https://gitlab.com/NMOF/NMOF2-Code . This code usesthe Heston Heston (1993) characteristic function and its integral form for the Call price to beused in the objective function (8.1). The numerical quadrature method is Gauss-Legendre. Wethus only need to replace it with our characteristic function Φ T , which consists of the functions25 a (for the SL model) and F (for the DRD model), whose .m -files are available from the Mathworks website, and the characteristic exponent ψ X of the underlying driving returns L´evyprocess X . We make two distinct choices for X : a Normal Inverse Gaussian (NIG; Barndorff-Nielsen 1997) process, with the parametrization from Cont and Tankov (2004); and a VarianceGamma (VG; Madan et al. 1998) process. We therefore calibrate a total of four anomalous dif-fusions models, namely: SL-VG, SL-NIG, SL-VG and DRD-VG, plus the two underlying L´evymodels VG and NIG (using the obvious modification of (8.1)).The parameters of the Differential Evolution algorithm have been chosen as follows (fordetails see (Gilli and Schumann, 2011)). A total of n G = 100 number of generations (thehalting condition) is produced, each with a population of n P = 10 d individuals spanning thesearch space, where d is the number of parameters to be estimated. Therefore n P = 40 forboth the DRD and SL models, and n P = 30 for the L´evy models. We used a mutation factor F = 0 . C R = 0 . 95 of a mutation surviving before the selection test isapplied. Further, every 10 generations we perform a Direct Search for solution improvement,each consisting of two searchers running a maximum of 100 iterations with a tolerance level of0.001. This choice of parameters follows commonly applied rules of thumb which seem to workwell for all models. The integral truncation for the Gauss-Legendre method is at ± 200 and theintegrand is evaluated at 50 nodes.The results are shown in Tables 1 to 3. In Table 1 we represent the baseline scenario: allthree models perfectly fit the synthetic L´evy market data. As expected, the β parameter in theSL and DRD models calibrates to one, and produces no improvement on the S calibration. Inthe scenarios with a persistent long-term skew, the total error for the L´evy model is greater thanthat of both the SL and DRD models with β < 1. Comparing among them the two scenarios weobserve as expected that for both models, β is smaller in the second scenario than in the firstone, owing to a steeper long-term skew in the latter case. This can be interpreted as an assetwith a more prolonged trade duration.Comparing errors across the models, the SL model shows a better fit in all cases. Howeverthis should not necessarily be interpreted as an overall superiority: the better calibration mightbe only due to the synthetic market data generated by a L´evy model, and the SL distributionsbeing closer to L´evy. We have proposed the use of anomalous diffusion processes in the context of option pricing,which allows to naturally incorporate trade durations between price moves. Using limits ofCTRWs whose inter-arrival times distribution obeys a power law to model asset returns, weanalysed the impact on the term structure of the returns distribution and on the correspondingimplied volatility. More specifically, the observed volatility skew persistence on the market can26e explained by a non-negligible impact of trade time randomness even in the long-term priceevolution.We analysed both cases when the price innovations are either dependent or independent fromthe waiting times between trades. Both models are consistent with the econometric observationthat shorter duration generates sharper variations in the price revisions. Finally, we remarkedthat even though the two models lead to similar large-maturity implied volatility properties,their different distributional properties produce rather different shapes of volatility surfaces.Numerical experiments confirm that for option pricing anomalous diffusions models have thepotential to capture the slow decay of the volatility skew while retaining the short-term goodproperties of pure L´evy models. References Andersen, L. and Lipton, A. (2012). Asymptotic exponential L´evy processes and their volatilitysmile: survey and new results. International Journal of Theoretical and Applied Finance ,16(1).Barndorff-Nielsen, O. E. (1997). Normal inverse Gaussian distributions and stochastic volatilitymodelling. Scandinavian Journal of Statistics , 24:1–13.Baumer, B., Meerschaert, M. M., and Scheffler, H. P. (2005). Space-time fractional derivativeoperator. Proceedings of the American Mathematical Society , 133:2273–2282.Becker-Kern, P., Meerschaert, M. M., and Scheffler, H. P. (2004). Limit theorems for coupledcontinuous time random walks. The Annals of Probability , 32:730–756.Bertoin, J. (1996). L´evy Processes . Cambridge University Press.Bertoin, J. (1997). Subordinators: examples and applications. In Lectures on Probability Theoryand Statistics , volume 1717, chapter 1, pages 1–91. Springer.Carr., P., Geman, H., Madan, D., and Yor, M. (2001). The fine structure of asset returns: anempirical investigation. The Journal of Business , 75:305–332.Cartea, A. and del-Castillo-Negrete, D. (2007). Fractional diffusion models of option prices inmarkets with jumps. Physica A , 374:749–763.Cartea, A. and Meyer-Brandis, T. (2010). How duration between trades of underlying securitiesaffects option prices. The Finance Review , 14:749–785.Chakrabarty, A. and Meerschaert, M. M. (2011). Tempered stable laws as random walks limits. Statistics and Probability Letters , 81:989–997.Cont, R. and Tankov, P. (2004). Financial Modelling with Jump Processes . Chapman & Hall.27elbaen, F. and Schachermayer, W. (1994). A general version of the fundamental theorem ofasset pricing. Mathematische Annalen , pages 463–520.Dellacherie, C. and Meyer, P. A. (1978). Probability and Potential , volume 1. Elsevier.Dufour, A. and Engle, R. (2000). Time and the price impact of a trade. The Journal of Finance ,55:2467–2498.Engle, R. (2000). The econometrics of ultra-high frequency data. Econometrica , 68:1–22.Engle, R. and Russell, J. (1998). Autoregressive conditional duration: a new model for irregularlyspaced transaction data. Econometrica , 66:1127–1162.Figueroa-L´opez, J. E., Forde, M., and Jacquier, A. (2011). The large-time smile and skew forexponential L´evy models. Working paper .Fries, C. and Torricelli, L. (2020). An analytical valuation framework for financial assets withtrading suspensions. To appear in SIFN .Gerhold, S., G¨ul¨um, I. C., and Pinter, A. (2016). Small-maturity asymptotics for at-the-moneyimplied volatility slope in L´evy models. Applied Mathematical Finance , 23:135–157.Gilli, M. and Schumann, E. (2011). Numerical Methods and Optimization in Finance . AcademicPress.Haubold, H. J., Mathai, A. M., and Saxena, R. K. (2011). Mittag-leffler functions and theirapplications. Journal of Applied Mathematics .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:327–343.Jacod, J. (1979). Calcul Stochastique et Probl`emes de Martingales . Lecture notes in Mathemat-ics. Springer.Jurlewicz, A., Kern, P., Meerschaert, M. M., and Scheffler, H. P. (2012). Fractional governingequations for coupled random walks. Computers and Mathematics with Applications , 64:3021–3036.Kyprianou (2014). Fluctuations of L´evy processes with Applications . Springer-Verlag BerlinHeidelberg.Laskin, N. (2003). The fractional Poisson process. Communication in Nonlinear Science , 8:201–213.Leonenko, N. N., Meerschaert, M. M., Schilling, R. L., and Sikorskii, A. (2014). Correlationstructure of time-changed L´evy processes. Communications in Applied and Industrial Math-ematics , 6. 28ewis, A. (2001). A simple option formula for general jump-diffusion and other exponentialL´evy processes. OptionCity.net .Luke, Y. L. (2012). The Special Functions and Their Approximations . Academy Press.Madan, D., Carr, P., and Chang, E. (1998). The variance gamma process and option pricing. European Finance Review , 2:79–105.Magdziarz, M. (2009). Black-Scholes formula in subdiffusive regime. Journal of StatisticalPhysics , 136:553–564.Mainardi, F., Gorenflo, R., and Scalas, E. (2004). A fractional generalization of the Poissonprocess. Vietnam Journal of Mathematics , 32:53–64.Mainardi, F., Raberto, M., Gorenflo, R., and Scalas, E. (2000). Fractional calculus andcontinuous-time finance II: the waiting time distribution. Physica A , 287:468–481.Meerschaert, M. M., Nane, E., and Vellaisamy, P. (2011). The fractional Poisson process andthe inverse stable subordinator. Electronic Journal in Probability , 16:1600–1620.Meerschaert, M. M., Nane, E., and Vellaisamy, P. (2013). Transient anomalous sub-diffusionson bounded domains. Proceedings of the American Mathematical Society , 141:699–710.Meerschaert, M. M. and Scalas, E. (2004). Coupled continuous time random walks in finance. Physica A , 41:114–118.Meerschaert, M. M. and Scheffler, H. P. (2004). Limit theorems for continuous-time randomwalks with infinite mean waiting times. Journal of Applied Probability , 41:623–638.Meerschaert, M. M. and Scheffler, H. P. (2008). Triangular array limits for continuous randomwalks. Stochastic Processes and their Applications , 118:1606–1633.Meerschaert, M. M. and Scheffler, H. P. (2010). Erratum to: ‘triangular array limits for contin-uous random walks’. Stochastic Processes and their Applications , 120:2520–2521.Meerschaert, M. M. and Straka, P. (2014). Semi-Markov approach to continuous time randomwalk limit processes. The Annals of Probability , 42:1699–1723.Rogers, L. and Tehranchi, M. (2010). Can the implied volatiliy move by parallel shifts? Financeand Stochastics , 14:235–248.Scalas, E., Gorenflo, R., and Mainardi, F. (2000). Fractional calculus and continuous-timefinance. Physica A , 284:376–384.Straka, P. and Henry, B. I. (2011). Lagging and leading coupled continuous time random walks,renewal times and their joint limits. Stochastic processes and their applications , 121:324–336.Torricelli, L. (2020). Trade duration risk in subdiffusive financial models. Physica A , 541.29 Figure 1: Paths of X H (blue) and H (green) in the SL model. β = 0 . β = 0 . X is a driftless Brownian motion with diffusion parameter σ = 0 . t for the SL and DRD model with the exponential,with β = 0 . t = 0 . 5. We used the compensated geometric Brownian motion characteristicexponent φ X ( z ) = σ ( z − i z ) / (cid:61) ( z ) = 1 / L Ht ∼ t B β, − β . For each t the total integral at some value x has the interpretation of the probability that the time for the background L´evy process X ranat most up to x .Figure 4: SL implied volatility surface based ongeometric Brownian motion; σ = 0 . β = 0 . 7. Figure 5: DRD implied volatility surface fromgeometric Brownian motion; σ = 0 . β = 0 . κ σ θ -0.2983 -0.1039 -0.2984 -0.1036 -0.2984 -0.1038 β - - 1.0000 0.9977 0.9999 0.9999RMSE 0.0084 0.0191 0.0084 0.0191 0.0084 0.0191Table 1: Calibration to the 1-month L´evy smile generated by the base model S . The parameters( κ, σ, θ ) = (0 . , . , − . 3) for the VG model and (0 . , . , − . 1) for the NIG model.31igure 6: SL implied volatility surface basedon a CGMY L´evy model, with C = 6 . , G =18 . , M = 32 . , Y = 0 . , β = 0 . 7. Figure 7: DRD implied volatility surface basedon a CGMY L´evy model, with C = 6 . , G =18 . , M = 32 . , Y = 0 . , β = 0 . β tends to one, with T = 0 . 25. Figure 11: Convergence of the DRD skew to theCGMY one as β tends to one, with T = 0 . β tends to one, for T = 0 . 75. Figure 13: Convergence of the DRD skew to theBS volatility as β tends to one, for T = 0 . κ σ θ -0.1696 -0.0556 -0.1739 -0.0546 -0.2824 -0.0995 β - - 0.8669 0.8837 0.7224 0.6271RMSE 0.3681 0.2061 0.2651 0.1729 0.2952 0.1791Table 2: Calibration to the 1-month and 1-year shifted L´evy smile generated by the basemodel S . The parameters ( κ, σ, θ ) are (0 . , . , − . 3) for VG and (0 . , . , − . 1) for NIG.Parameter L´evy SL DRDVG NIG VG NIG VG NIG κ σ θ -0.1354 -0.0785 -0.1571 -0.0711 -0.2566 -0.1104 β - - 0.8305 0.8634 0.6546 0.5704RMSE 0.4857 0.2705 0.3612 0.2307 0.4157 0.2442Table 3: Calibration to 1-month and 18-month shifted L´evy smiles generated by S . Theparameters ( κ, σ, θ ) are (0 . , . , − . 3) for the VG model and (0 .. 3) for the VG model and (0 .. , .. 3) for the VG model and (0 .. , .. , − ..