Lévy interest rate models with a long memory

Hainaut, Donatien Article Lévy interest rate models with a long memory Risks Provided in Cooperation with: MDPI – Multidisciplinary Digital Publishing Institute, Basel Suggested Citation: Hainaut, Donatien (2021) : Lévy interest rate models with a long memory, Risks, ISSN 2227-9091, MDPI, Basel, Vol. 10, Iss. 1, pp. 1-28, https://doi.org/10.3390/risks10010002 This Version is available at: https://hdl.handle.net/10419/258313 Standard-Nutzungsbedingungen: Die Dokumente auf EconStor dürfen zu eigenen wissenschaftlichen Zwecken und zum Privatgebrauch gespeichert und kopiert werden. Sie dürfen die Dokumente nicht für öffentliche oder kommerzielle Zwecke vervielfältigen, öffentlich ausstellen, öffentlich zugänglich machen, vertreiben oder anderweitig nutzen. Sofern die Verfasser die Dokumente unter Open-Content-Lizenzen (insbesondere CC-Lizenzen) zur Verfügung gestellt haben sollten, gelten abweichend von diesen Nutzungsbedingungen die in der dort genannten Lizenz gewährten Nutzungsrechte. Terms of use: Documents in EconStor may be saved and copied for your personal and scholarly purposes. You are not to copy documents for public or commercial purposes, to exhibit the documents publicly, to make them publicly available on the internet, or to distribute or otherwise use the documents in public. If the documents have been made available under an Open Content Licence (especially Creative Commons Licences), you may exercise further usage rights as specified in the indicated licence. https://creativecommons.org/licenses/by/4.0/   Citation: Hainaut, Donatien. 2022. Lévy Interest Rate Models with a Long Memory. Risks 10: 2. https:// doi.org/10.3390/risks10010002 Academic Editor: Mogens Steffensen Received: 9 November 2021 Accepted: 20 December 2021 Published: 23 December 2021 Publisher’s Note: MDPI stays neutral with regard to jurisdictional claims in published maps and institutional affiliations. Copyright: © 2021 by the author. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https:// creativecommons.org/licenses/by/ 4.0/). risks Article Lévy Interest Rate Models with a Long Memory Donatien Hainaut † UCLouvain, LIDAM, Louvain-La-Neueve, 1348 Ottignies-Louvain-la-Neuve, Belgium; [email protected] † Current address: 20 Voie du Roman Pays, Louvain-La-Neueve, 1348 Ottignies-Louvain-la-Neuve, Belgium. Abstract: This article proposes an interest rate model ruled by mean reverting Lévy processes with a sub-exponential memory of their sample path. This feature is achieved by considering an Ornstein–Uhlenbeck process in which the exponential decaying kernel is replaced by a Mittag–Leffler function. Based on a representation in term of an infinite dimensional Markov processes, we present the main characteristics of bonds and short-term rates in this setting. Their dynamics under risk neutral and forward measures are studied. Finally, bond options are valued with a discretization scheme and a discrete Fourier’s transform. Keywords: interest rate; Lévy process; Mittag–Leffler function; mean reverting process 1. Introduction From the 1980s to current years, many interest rate models were proposed in the literature. Their common aim is to explain changes in bond or swap quotes and to replicate risks within the interest rates market. Three dominating frameworks coexist: short-term rate, forward rate and the Libor market models. In this last approach proposed by Brace et al. (1997), interest rates are driven by geometric processes. In the forward rate model pioneered by Heath et al. (1992), the term structure of rates is specified through instantaneous forward rates. Mercurio and Moraleda (2000) and Falini (2010) proposed a forward model with a humped structure of volatilities. Li et al. (2020) developed a forward rate model unifying most existing Gaussian models for interest rates. Short-term rate models such as those promoted by Hull and White (1990) specify a mean reverting dynamic for the instantaneous risk-free rate. The framework developed in this article belongs to this third category. The family of short-term rate models gathers multiple frameworks and we refer the reader, e.g., to Boero and Torricelli (1996) or to Schmidt (2011) for a review. On the other hand, various processes were proposed to explain the interest rate risk. The literature is too vast to be exhaustive. We restrict our attention to contributions related to this work and do not attempt at providing a general overview, referring instead to Brigo and Mercurio (2006) for detailed accounts on the topic. For instance, Eberlein and Raible (1999) were among the first to propose a short-term rate model driven by Lévy processes. In their setting, the short-term rate is ruled by Ornstein–Uhlenbeck processes reverting in an exponential manner to a mean level. Within this framework, Eberlein and Kluge (2005) derived analytical formulae for the prices of caps and floors using bilateral Laplace transforms. In a similar setting, Hainaut and MacGilchrist (2010) proposed a pentanomial tree for pricing derivatives. Hainaut (2013) studied the properties of a Gaussian shortterm rate with a Markov Switching Multifractal volatility. Moreno and Platania (2015) proposed a square-root model replicating economic cycles in the dynamics of interest rates. Hainaut (2016) and Njike and Hainaut (2020) developed unicurve and multicurve models in which the short-term rate is exposed to self-excitating jumps. Fontana et al. (2020) studied a modelling framework exploiting the the self-exciting behavior of continuous-state branching processes with immigration (CBI). A wide majority of short-term rate models are driven by Markov mean reverting processes. In this setting, the interest rate depends on its path through the integral of past Risks 2022,10, 2. https://doi.org/10.3390/risks10010002 https://www.mdpi.com/journal/risks Risks 2022,10, 2 2 of 28 occurrences weighted by an exponential decaying function, called the memory kernel. In this framework, the influence of previous variations of interest rates decreases exponentially over time. In this article, we consider an alternative that consists in replacing this memory kernel by a Mittag–Leffler decaying function. This function is sub-exponential and may be seen as a generalization of the exponential. In this setting, the interest rate remembers its previous occurrences for a longer period than in the exponential framework. The motivations for studying such a model are multiple. Firstly, the features replicate the persistence empirically observed in short term rates, as illustrated in Section 3, or in nominal yields and inflation, as shown by Golinski and Zaffaroni (2016). Secondly, the proposed model has a sufficient level of analytical tractability for most of financial applications such as option pricing. Finally, it also offers a novel alternative to existing approaches based on the fractional Brownian motion for introducing long memory, as for instances studied in Cheridito et al. (2003)orMaejima and Yamamoto (2003). A consequence of this change of memory kernel is that the interest rate process is not Markovian anymore. Bond prices depend then on the sample path of past rates. Nevertheless, we take advantage of the representation of the memory kernel as a Laplace Stieltjes integral to rewrite the short-term rate as an infinite dimensional Markov process. The last part of the article studies the properties of the model under a forward measure and proposes a pricing method for bond options. 2. A Lévy Model with a Mittag Leffler Kernel We consider that the short-term rate is driven by d∈N independent Lévy processes, denoted by (L(j) t)t≥0 for j= 1, . . . , d . These processes are defined on a probability space Ω , endowed with the natural filtration (Ft)t≥0 and a risk neutral measure, Q . Each process (L(j) t)t≥0 has independent and stationary increments. It is fully described by a triplet µj,σ2 j,νj where µj∈R , σj∈R+ and νj( . ) is a Lévy measure. The moment generating function (mgf) is denoted by φ(j) t(ω)for j=1, . . . , dand is equal to the following: φ(j) t(ω):=Eexpω(L(j) t−L(j) 0)| F0 =exptµjω+1 2ω2σ2 j+ZReωz−1−ωz1|z|<eνj(dz) (1) =exptψj(ω), where e> 0. Notice that if the Lévy measure, νj( . ) , has no singurality at z= 0, we can set e to zero. The function ψj(ω) is the characteristic exponent of this Lévy process that is always defined for pure imaginary numbers. Without loss of generality, we assume that L(j) 0= 0. According to the Lévy Itô decomposition, each L(j) t is the sum of three components: a deterministic drift µjt , a Brownian motion with variance σ2 j and a jump process, Jj(t , z) , of intensity νj(dz) well defined on [−∞ , 0 )∪( 0, +∞] . This Lévy measure is such that the probability of observing k jumps between [τ1 , τ2] of a size included in a set B⊂R\{0} is given by the following: PJj([τ1,τ2]×B) = k=e−Rτ2 τ1RBνj(dz)dt Rτ2 τ1RBνj(dz)dtk k!, (2) for j= 1, . . . , d . If W(j) tt≥0 for j= 1, . . . , d , are independent Brownian motions on (Ω,F,Q),L(j) tmay be split as the sum of a drift, a Brownian motion and a jump process: dL(j) t=µjdt +σjdW(j) t+ZRz˜ Jj(dt ,dz), (3) Risks 2022,10, 2 3 of 28 where the following is the case. ˜ Jj(dt ,dz) = Jj(dt ,dz)−1|z|<eνj(dz)dt . Without loss of generality, we assume that EL(j) t= 0. This constraint implies that the drift compensates jumps larger than ein absolute value: µj=−Z|z|≥ezνj(dz),j=1, . . . , d. From Cont and Tankov (2003, Lemma 15.1, p. 482), for any integrable function f:R→C, the following relation holds. EeωRt sf(u)dL(j) u|Fs=expZt sψj(ωf(u))du,j=1, . . . , d. (4) This can be proved by approaching f( . ) with a stepwise function and by using the property of independence of increments. This property will be useful in later developments. We postulate that the risk-free rate, (rt)t≥0 , is the sum of a deterministic function ϕ(t): R→Rand of dprocesses: X(d) tt≥0: rt=ϕ(t) + d ∑ j=1 X(j) t, (5) where X(j) tis such that the following is the case. X(j) t=gj(t)X(j) 0+Zt 0gj(t−u)dL(j) u. (6) The function gj( . ):R+→R is continuously decreasing with an initial value gj( 0 ) = 1. As gj( . ) determines the influence of past realizations of L(j) t on the current value of the process, we call it “memory kernel” or “kernel function” in the remainder of this article. When the function gj(t) decreases exponentially, the processes X(j) tj=1,...,d are Markov and mean-reverting as recalled in Appendix A. Such a model was, e.g., developed in Hainaut and MacGilchrist (2010) for pricing interest rate derivatives and will serve as benchmark in numerical illustrations. However, X(j) tt≥0 and (rt)t≥0 are in general not Markov for any other non-exponential kernels. This feature makes difficult the evaluation of bond prices at a given time t> 0 as their value depends on the entire sample path of interest rates up to t . Nevertheless, for a few kernel functions that admit a representation as a Laplace–Stieltjes integral, we will show that the interest rate may be represented as an infinite dimensional Markov process. In particular, we consider two types of kernel functions, both based on the Mittag–Leffler function, denoted by Eα(.)where α∈[0, 1]. gj(t) = Eαj(−βjt)or gj(t) = Eαj(−βjtαj). (7) In order to understand the motivation for working with such functions, we need to review the main properties of the Mittag–Leffler function. The Mittag–Leffler function of order α> 0 plays a fundamental role in the fractional calculus and can be considered as an extension of the exponential function. It is defined as an infinite sum. Eα(t) = ∞ ∑ n=0 tn Γ(nα+1). Risks 2022,10, 2 4 of 28 In this article, we assume that α∈[ 0, 1 ] . For this range of values, Eα(−t) is an intermediary between the power and the exponential decreasing functions. E0(−t)=1 1+t∀|t|<1 and E1(−t) = e−t. We refer the reader to the book of Gorenflo et al. (2020) for a detailed presentation of this function. For a complete presentation of fractional calculus, we recommend the book of Baleanu et al. (2012) and the book of Kochubei and Luchko (2019). Let us recall that a function f:( 0, ∞)→∞ is called completely monotonic if it possesses derivatives f(n)(t) of any order n=0, 1, 2 . . . and the derivatives are alternating in sign. (−1)nf(n)(t)≥0∀t∈((0, ∞). The above property is equivalent to the existence of a representation of the function f in the form of a Laplace–Stieltjes integral with non-decreasing density and non-negative measure dγ(.)such that the following is the case. f(t) = Z∞ 0e−utdγ(u). The Mittag–Leffler function of negative argument Eα(−βt) is completely monotonic for all 0 ≤α≤ 1 and β∈R+ . A proof of this result is available in Gorenflo et al. (2020) on page 47. The authors show the following: Eα(−βt) = Z∞ 0e−tudγα,β(u), (8) where the derivative of γα,β(u)is equal to the following. dγα,β(u) du =1 πα ∞ ∑ k=1 (−1)k−1 k!sin(παk)Γ(αk+1)uk−1 βk. (9) Given that Eα( 0 ) = 1, we have R∞ 0dγα,β(u) = 1 and γα,β( . ) is then a measure of probability on R+ . In the remainder of this article, Eα(−βt) is called the decreasing Mittag– Leffler (ML) kernel. The DML behaves at short-term as follows: Eα(−βt) = 1−βt Γ(1+α)+· · · ∼ exp−βt Γ(1+α), when t→ 0. From Haubold et al. (2011), we known that when t→∞ , DML converges to the following. Eα(−βt)∼(βt)−1 Γ(1−α). The function Eα(−βt) interpolates for intermediate times t between the decreasing exponential and the inverse power law. The exponential models fast decay for small time t , whereas the asymptotic inverse power law entails a slow decrease at long term. This point is illustrated in the left plot of Figure 1, which compares the ML kernel to the exponential and inverse power functions. Risks 2022,10, 2 5 of 28 0 5 10 15 20 25 30 0.0 0.2 0.4 0.6 0.8 1.0 Mittag Leffler Tv y ML Exp. Inverse 0 5 10 15 20 25 30 0.0 0.2 0.4 0.6 0.8 1.0 Power Mittag Leffler Tv y PML Stretched exp. Power Figure 1. Left and right plots: Comparison of the Mittag–Leffler (ML), Eα(−βt) , and power Mittag– Leffler (PML), Eα(−βtα) kernels with exponential and power decreasing functions ( α= 0.7 and β=0.5). The Mittag–Leffler function is also related to the fractional calculus. In order to explain this link, we recall that Caputo’s fractional derivative of order α∈( 0, 1 ) for a function h(t) : R+→R,C1with respect to tis defined by the following. ∂α ∂tαh(t) = 1 Γ(1−α)Zt 0(t−s)−α∂ ∂sh(s)ds (10) When α= 1, this derivative corresponds to the first order derivative. The solution of the fractional differential equation of the following: ∂α ∂tαy(t) = −βy(t)0<α<1 with the initial condition y( 0 ) = b0 is y(t) = Eα(−βtα) . We abusively call this function the power decreasing Mittag–Leffler kernel (PML). The Laplace transform of Eα(−βtα) is given by sα−1 sα+β where Re(s)>|β|1 α . By inverting this transform, we can prove as in Gorenflo and Mainardi (1997) that the PML is the Laplace transform of the following: Eα(−βtα) = Z∞ 0e−ut dγp α,β(u), (11) where the derivative of γp α,β(u)is given by the following. dγp α,β(u) du =1 π βuα−1sin(απ) u2α+2βuαcos(απ) + β2. (12) Since Eα( 0 ) = 1, we have R∞ 0dγp α,β(u) = 1 and γp α,β(u) is a measure of probability on R+. As underlined by Mainardi (2020), the PML behaves at short-term as follows: Eα(−βtα) = 1−βtα Γ(1+α)+· · · ∼ exp−βtα Γ(1+α). when t→ 0. From Erdélyi et al. (1955), we known that when t→∞ , PML converges to the following. Eα(−βtα)∼β1/αt−α Γ(1−α). Risks 2022,10, 2 6 of 28 As a consequence, the function Eα(−βtα) interpolates for intermediate times t between the stretched exponential and the negative power law. The stretched exponential models the very fast decay at short-term whereas the asymptotic power law is due to the very slow decay for large time t. The right plot of Figure 1illustrates this convergence. Figure 2presents simulated sample paths of the short-term rate with a one dimensional ( d= 1) model ruled by a Brownian motion. These paths are computed for various α with the same random occurrences in order to make a comparison feasible. For the ML kernel, the trajectories are nearly similar. An analysis of figures reveals that the sample path is smoother for α= 0.50 than for α= 0.90. This trend becomes visible if we choose the highest value for β . For PML sample paths, the difference is clearly visible. Decreasing α reduces the volatility of rates and smooths the sample path. 0.0 0.1 0.2 0.3 0.4 0.5 0.6 −0.010 0.000 0.010 M.L. Time r(t) alpha=0.90 alpha=0.50 alpha=0.10 0.0 0.1 0.2 0.3 0.4 0.5 0.6 −0.010 0.000 0.010 Power M.L. Time r(t) alpha=0.90 alpha=0.50 alpha=0.10 Figure 2. Upper and lower plots: Comparison of sample paths of (rt)t≥0 with Mittag–Leffler (ML) and power Mittag—Leffler (PML) kernels. d= 1, β= 3, X(1) 0= 0, ϕ(t) = 0 and Lt are Brownian motions with σ1=0.05. To conclude this section, we present the first two moments of the short-term rate and its autocovariance function. The ML and PML kernels defined in Equation (7) are continuously decreasing and integrable functions on any bounded interval of R . From Equation (4), the mgf of X(j) tfor j=1, . . . , dis then equal to the following. EeωX(j) t|F0=expωgj(t)X(j) 0+Zt 0ψjωgj(t−u)du. Deriving mgf allows us to find the first moments of X(j) t conditionally to the initial information. On the other hand, from the following representation: rt=ϕ(t) + d ∑ j=1 gj(t)X(j) 0+ d ∑ j=1Zs 0gj(t−u)dL(j) u+ d ∑ j=1Zt sgj(t−u)dL(j) u, (13) of rt , we infer that the conditional expectation of the short-term rate is given by the following. E(rt|Fs)=ϕ(t) + d ∑ j=1 gj(t)X(j) 0+ d ∑ j=1Zs 0gj(t−u)dL(j) u, (14) Risks 2022,10, 2 7 of 28 On the other hand, the conditional variance of rt is equal to the sum of variances of Lévy processes times the integral of the squared kernels. V(rt|Fs)= d ∑ k=1σ2 j+ZRz2νj(dz)Zt sgj(t−u)2du . (15) Unfortunately the integral of gj( . )2 does not admit a closed form expression for the ML and PML kernels. Nevertheless, they can be numerically estimated. A direct calculation allows us to infer that autocovariance between rtand rufor u≤tis given by: C(rtru|Fs)= d ∑ k=1 EZt sgj(t−v)dL(j) vZu sgj(u−v)dL(j) v(16) = d ∑ k=1σ2 j+ZRz2νj(dz)Zu sgj(t−v)gj(u−v)dv . Contrary to the exponential case, this covariance function between rt and rt−∆ does not admit any closed form expression when t→∞. 3. Empirical Motivation This short section provides some empirical arguments motivating the developments performed in this article, particularly the choice of a ML or PML memory kernel. We fit univariate Gaussian models ( d= 1) with exponential ML and PML kernels to the Eonia time series from 13 March 2016 to 1 August 2019. The dataset counts nobs = 866 daily observations. We select this time window mainly because Eonia was relatively stable during this period, as illustrated in the right plot of Figure 3. This allows us to assume that the trend function, ϕ(t) , is constant. Considering a larger dataset would raise the question of modelling a non-linear trend in view of the left plot of Figure 3. 2000 2005 2010 2015 2020 012345 13/3/2000 to 1/2/2021 Date Eonia % 2017 2018 2019 −0.38 −0.34 −0.30 −0.26 13/3/2016−1/8/2019 Date Eonia % Figure 3. ( Left plot): Eonia time series from 13 March 2000 to 1 February 2021. ( Right plot): dataset used to estimate parameters. The dates of observations are denoted by t0 , t1 , . . . , tnobs−1 , whereas the step of time between two successive observations is ∆t . We first consider a one dimensional exponential kernel model as detailed in Appendix A. Under the assumption that ϕ(tk) = ϕ(tk+∆t) and X0=0, we find the following: ∆Lk+1=Ltk+1−Ltk≈rtk+1−rtk− k ∑ j=0e−κ(tk+1−tj)−e−κ(tk−tj)∆Lj, for k= 1, . . . , nobs − 1. Under the assumption that Lt is a Brownian motion without drift and of variance σ2 , parameters are estimated by log-likelihood maximization. We next Risks 2022,10, 2 8 of 28 consider models with ML and PML kernel functions. In this cases, the variations of the underlying Lévy process are approached by the following: ∆Lk+1=Ltk+1−Ltk≈rtk+1−rtk− k ∑ j=0g(tk+1−tj)−g(tk−tj)∆Lj, for k= 1, . . . , nobs − 1 where g( .) is either the ML or the PML function. Under the assumption of normality, parameters are estimated by log-likelihood maximization whereas the goodness of fit is measured with the Akaike Information Criterion (AIC). The results of this procedure are reported in Table 1. In terms of log-likelihood, exponential and ML models achieve the same goodness of fit but the exponential model is more parsimonious and preferred. They share the same estimated volatilities but the parameters of reversion κ and β slightly differ. The best fit, based on the AIC, is obtained with a power Mittag–Leffler kernel. The parameter α that drives the decay of the memory is small. This means that the process has a longer-term memory than the exponential model. On the other hand, the mean reversion speed is significantly higher than the ones of other models. We recall to the reader that estimated parameters reflect the dynamic of rt under the real measure P and not the risk neutral one. Table 1. Parameter estimates by log-likelihood maximization of one dimensional models with exponential, ML and PML kernels. Exponential ML PML σ1.4763 ×10−3σ1.4767 ×10−31.3927×10−3 κ0.6926 α0.9954 0.1924 β0.7218 1.2171 Log-likelihood 6809.75 Log-likelihood 6809.93 6861.15 AIC −13,615.50 −13,613.86 −13,716.30 4. Alternative Formulation If the memory is not exponential decaying, the interest rate process (rt)t≥0 is not Markovian. This feature makes difficult the evaluation of bond prices or interest rate derivatives because their values depend on the entire history of the rate process. Nevertheless, if the memory kernel admits a representation as a Laplace–Stieltjes integral, we can convert the short-term rate into an infinite dimensional Markov process. Approximating this process allows us to infer the dynamic of bond prices. As developments are similar in both cases, we denote by γj( . ) the measures γ(j) αj,βj( . ) or γp(j) αj,βj( . ) of Eαj(−βjt) and Eαj(−βjtαj) . Let us recall that the function gj(t) = R∞ 0e−ξtdγj(ξ) . Using the Fubini’s theorem, the processes X(j) tare then rewritten as follows: X(j) t=X(j) 0Z∞ 0e−tξdγj(ξ) + Zt 0Z∞ 0e−ξ(t−u)dγj(ξ)dL(j) u =Z∞ 0e−tξY(j,ξ) 0+Zt 0e−ξ(t−u)dL(j) u | {z } Y(j,ξ) t dγj(ξ) for j= 1, . . . , d . The Y(j,ξ) t are Ornstein–Uhlenbeck processes for j= 1, . . . , d and for all ξ∈R+ . Their initial values are X(j) 0=Y(j,ξ) 0 , and they obey to the following stochastic differential equation (SDE). dY(j,ξ) t=−ξY(j,ξ) tdt +dL(j) t. (17) Risks 2022,10, 2 15 of 28 where ϕ(n)(t)is the discretized version of ϕ(t). ϕ(n)(t) = −∂ ∂tln(P(0, t))− d ∑ j=1 ψj − n ∑ k=1 B(j,k)(0, t)m(j) k!. (37) As the number of processes in the discretized model is finite, we can apply Itô’s lemma in order to establish the dynamic of P(n)(s,t). Proposition 6. The zero coupon bond price P(n)(s , t) is a geometric Lévy processes solution of the SDE: dP(n)(s,t) P(n)(s,t)=r(n) sds − d ∑ j=1 σj n ∑ k=1 B(j,k)(s,t)m(j) kdW(j) s(38) + d ∑ j=1ZRe−z∑n k=1m(j) kB(j,k)(s,t)−1Jj(ds ,dz)−νj(dz)ds. with the terminal condition P(n)(t,t) = 1. Proof. If we remember that dY(j,k) s=dL(j) s−b(j) kY(j,k) sds , the Itô’s lemma provides us with the following differential for P(n)(s,t). dP(n)(s,t) = ∂P(n)(s,t) ∂sds − d ∑ j=1 n ∑ k=1 ∂P(n)(s,t) ∂Y(j,k)b(j) kY(j,k) sds (39) + d ∑ j=1 n ∑ k=1 ∂P(n)(s,t) ∂Y(j,k)µjds + d ∑ j=1 n ∑ k=1 ∂P(n)(s,t) ∂Y(j,k)ZRz˜ Jj(dt ,dz) + d ∑ j=1 n ∑ k=1 ∂P(n)(s,t) ∂Y(j,k)σjdW(j) s+1 2 d ∑ j=1 n ∑ k=1 n ∑ l=1 ∂2P(n)(s,t) ∂Y(j,k)∂Y(j,l)σ2 jds + d ∑ j=1ZRP(n)(s,t)e−z∑n k=1m(j) kB(j,k)(s,t)−1−z n ∑ k=1 ∂P(n)(s,t) ∂Y(j,k)Jj(ds ,dz). The first and second order partial derivatives of P(n)(s,t)are equal to the following. ∂P(n)(s,t) ∂Y(j,k)=−P(n)(s,t)B(j,k)(s,t)m(j) k, (40) ∂2P(n)(s,t) ∂Y(j,k)∂Y(i,l)=P(n)(s,t)B(j,k)(s,t)B(i,l)(s,t)m(j) km(i) l. (41) On the other hand, its partial derivative with respect to time is given by the following: ∂P(n)(s,t) ∂s=P(n)(s,t) ϕ(n)(s) + d ∑ j=1 n ∑ k=1 Y(j,k) se−b(j) k(t−s)m(j) k(42) − d ∑ j=1 ψj − n ∑ k=1 B(j,k)(s,t)m(j) k!!, Risks 2022,10, 2 16 of 28 where the characteristic exponent is developed as the following sum. ψj − n ∑ k=1 B(j,k)(s,t)m(j) k!=−µj n ∑ k=1 B(j,k)(s,t)m(j) k+1 2 n ∑ k=1 B(j,k)(s,t)m(j) k!2 σ2 j(43) ZR e−∑n k=1B(j,k)(s,t)m(j) kz−1+ n ∑ k=1 B(j,k)(s,t)m(j) k!z1|z|<e!νj(dz). Combining Equations (39)–(43) allows us to rewrite the dynamic of bond prices as follows. dP(n)(s,t) = P(n)(s,t) ϕ(n)(s) + d ∑ j=1 n ∑ k=1 Y(j,k) se−b(j) k(t−s)m(j) k!ds +P(n)(s,t) d ∑ j=1 n ∑ k=1 Y(j,k) sB(j,k)(s,t)m(j) kb(j) kds −P(n)(s,t) d ∑ j=1 σj n ∑ k=1 B(k)(s,t)m(n) kdW(j) s +P(n)(s,t) d ∑ j=1ZRe−z∑n k=1m(j) kB(j,k)(s,t)−1Jj(ds ,dz)−νj(dz)ds. We infer the result if we remember Equation (35) and notice that the following is the case. n ∑ k=1 Y(j,k) se−b(j) k(t−s)m(j) k+ n ∑ k=1 Y(j,k) sB(k)(s,t)m(j) kb(j) k= n ∑ k=1 Y(j,k) sm(j) k. As by construction limn→∞γ(n) j(z) = γj(z) , we immediately infer the dynamic of the bond price in the non-discretized model by considering the limit of Equation (38). Corollary 1. If the short-term rate is driven by Equation (19), the bond price is solution of the following SDE: dP(s,t) P(s,t)=rsds − d ∑ j=1 σjhj(s,t)dW(j) s(44) + d ∑ j=1ZRe−zhj(s,t)−1Jj(ds ,dz)−νj(dz)ds. where hj(s,t) = R∞ 0B(ξ)(s,t)dγj(ξ)and with the terminal condition P(t,t) = 1. In order to understand the role played by the memory kernel in the evolution of the term structure of interest rates, we analyze the expectation and variance of bond yields. The bond yield of a zero-coupon bond in the discretized model is defined as follows. y(n)(s,t) = −ln P(n)(s,t) t−s. Risks 2022,10, 2 17 of 28 From Equation (36), we infer that this yield is proportional to a weighted sum of mean reverting processes Y(j,k) s . Since EY(j,k) s|F0= 0 by construction, the expected yield conditionally to the initial filtration is equal to the following. Ey(n)(s,t)| F0=1 t−s Zt sϕ(n)(u)du − d ∑ j=1Zt sψj − n ∑ k=1 B(j,k)(v,t)m(j) k!dv!. (45) The variance of the yield is the sum of d variances of a linear combination of processes Y(j,k) s. If we remember that the covariance Y(j,k) sand Y(l,k) sis equal to the following: CY(j,k) sY(j,l) s|F0=σ2 j+ZRz2νj(dz)Zs 0e−b(j) k+b(j) l(s−u)du (46) =σ2 j+RRz2νj(dz) b(j) k+b(j) l1−e−b(j) k+b(j) ls, the variance of the yield, conditionally to F0, is the following triple sum. Vy(n)(s,t)| F0= d ∑ j=1 V  n ∑ k=1 B(j,k)(s,t)m(j) k t−sY(j,k) s|F0 (47) = d ∑ j=1  n ∑ k=1 n ∑ l=1 B(j,k)(s,t)m(j) k t−s B(j,l)(s,t)m(j) l t−sCY(j,k) sY(j,l) s|F0 . To illustrate these results, we consider an univariate jump-diffusion model ( d= 1) in which a single Lévy process is ruled by the following SDE. dL(1) t=−λ1η1dt +σ1dW(1) t+η1dN(1) t. (48) N(1) t is a Poisson process with constant intensity λ1 . The characteristic exponent is in this case equal to the following: ψ1(ω) = −λ1η1ω+1 2ω2σ2 1+λ1(eωη1−1), (49) whereas the Lévy measure is ν1(z) = λ1δη1(z) where δη1(z) is the Dirac measure located at z=η1. The instantaneous variance of dL(1) tis equal to the following. σ2 1+ZRz2ν1(dz) = σ2 1+λ1η2 1. Of course, we can consider other types of Lévy processes such as the variance gamma and the Normal inverse gaussian, but this does not fundamentally modify the conclusions drawn in this section. We first fit the curve ϕ(t) to the term structure of zero-coupon bond yields, bootstrapped from the ICE swap rates on the 26/2/21 at noon. Since this function is proportional to instantaneous forward rates, we have to interpolate bond yields. For this purpose, we use a Nelson–Siegel (NS) model. Appendix Brecalls this model and reports the swap rates curve on the 26/2/21 and parameter estimates. There exist more advanced interpolation methods such as, e.g., splines or kriging techniques as detailed in Cousin et al. (2016) but the NS model is sufficiently accurate for our purposes to understand the dynamic of yields over time. We calculate the term structure of expected yields in 5, 10 and 15 years with n=40 atoms . Increasing n does not significantly modify the results. A sensitivity analysis with respect to the number of atoms is proposed in Section 7. For the parameters, we set β= 1.5, λ1= 0.5, σ1= 0.01 and η1=− 0.0002, which are possible realistic values in Risks 2022,10, 2 18 of 28 view of results from Section 3. We remind the reader that α is the parameter tuning the memory of the model. If α= 1, the memory kernel is exponential and the model forgets in an exponential manner past fluctuations of interest rates. On the other hand, for lower values of α , the model forgets these variations according to a power decaying function. Figures 5and 6 show expected yields computed with the ML and PML models for α= 0.5 and α= 0.9. We use a discretization scheme with 40 atoms that cover ML and PML density up to their 90% percentiles. For both models, the range between short and long-term yields is narrowing with the time horizon. In the ML model, the lower α is, the longer the memory and the quicker the convergence to a bumped (but nearly flat) yield curve. In the PML model, expected yield curves with a low α dominate those with a high α and have a more pronounced curvature. 0 10 20 30 40 0.002 0.003 0.004 0.005 0.006 E(yield) in 5 years Maturity yield alpha=0.9 alpha=0.5 0 10 20 30 40 0.0054 0.0056 0.0058 0.0060 0.0062 0.0064 E(yield) in 10 years Maturity yield alpha=0.9 alpha=0.5 0 10 20 30 40 0.00625 0.00630 0.00635 0.00640 0.00645 E(yield) in 15 years Maturity yield alpha=0.9 alpha=0.5 Figure 5. Expected yield curves in 5, 10 and 15 years; α= 0.5 and 0.9; β= 1.5; λ1 = 0.5; η1=− 0.002; σ1=0.01. ML kernel with n=40 atoms. 0 10 20 30 40 0.002 0.003 0.004 0.005 0.006 E(yield) in 5 years Maturity yield alpha=0.9 alpha=0.5 0 10 20 30 40 0.0056 0.0058 0.0060 0.0062 0.0064 E(yield) in 10 years Maturity yield alpha=0.9 alpha=0.5 0 10 20 30 40 0.00635 0.00640 0.00645 0.00650 E(yield) in 15 years Maturity yield alpha=0.9 alpha=0.5 Figure 6. Expected yield curves in 5, 10 and 15 years; α= 0.5 and 0.9; β= 1.5; λ1 = 0.5; η1=− 0.002; σ1=0.01. PML kernel with n=40 atoms. Figures 7and 8report the expected standard deviations of future yields calculated with ML and PML models. In both cases, the term structures of expected standard deviations do not significantly evolve with the time horizon and are decreasing functions of the maturity. We also notice that standard deviations are, respectively, inversely and directly proportional to αin ML and PML models. Risks 2022,10, 2 19 of 28 0 10 20 30 40 0.001 0.002 0.003 0.004 Std(yield) in 5 years Maturity yield alpha=0.9 alpha=0.5 0 10 20 30 40 0.001 0.002 0.003 0.004 Std(yield) in 10 years Maturity yield alpha=0.9 alpha=0.5 0 10 20 30 40 0.001 0.002 0.003 0.004 Std(yield) in 15 years Maturity yield alpha=0.9 alpha=0.5 Figure 7. Expected standard deviations of yields in 5, 10 and 15 years; α= 0.5 and 0.9; β= 1.5; λ1= 0.5; η1=−0.002; σ1=0.01. ML kernel with n=40 atoms. 0 10 20 30 40 0.0000 0.0005 0.0010 0.0015 0.0020 0.0025 0.0030 Std(yield) in 5 years Maturity yield alpha=0.9 alpha=0.5 0 10 20 30 40 0.0000 0.0005 0.0010 0.0015 0.0020 0.0025 0.0030 Std(yield) in 10 years Maturity yield alpha=0.9 alpha=0.5 0 10 20 30 40 0.0000 0.0005 0.0010 0.0015 0.0020 0.0025 0.0030 Std(yield) in 15 years Maturity yield alpha=0.9 alpha=0.5 Figure 8. Expected standard deviations of yields in 5, 10 and 15 years; α= 0.5 and 0.9; β= 1.5; λ1= 0.5; η1=−0.002; σ1=0.01. PML kernel with n=40 atoms. 7. Pricing of Bond Options The purpose of this section is to illustrate how long-term memory models can be used for the pricing of derivatives. We investigate this for European options written on a zero coupon bond. Our choice is motivated by the fact that these are the building blocks of many structured products. The underlying bond has a maturity noted t2 while the option expires at time t1≤t2 with a strike price K . Without loss of generality we focus on a call option but all developments remain valid for any other type of European options. The value at time s≤t1≤t2 of the call option is denoted by C(s) and is equal to the expected discounted payoff under the risk neutral measure Q. C(s) = Ee−Rt1 srudu(P(t1,t2)−K)+| Fs. This price does not admit any closed form expression. Nevertheless, we can apply the standard technique of changes of numeraire combined with a discrete Fourier’s transform (DFT) to evaluate this call option. The first stage consists to define a forward measure, denoted by F that has P(s , t1) as numeraire, i.e., every assets discounted by P(s , t1) is a martingale. If S(0) t=eRt 0rudu is the cash account (numeraire of Q ), the Radon–Nykodym derivative defining this change of measure is dF dQt1 =S(0) 0 S(0) t1P(0,t1) . Using standard arguments, the call price is equal to the product of a discount bond and of the expected payoff under this forward measure. C(s) = P(s,t1)EF(P(t1,t2)−K)+| Fs. (50) Risks 2022,10, 2 20 of 28 The next proposition details the dynamics of underlying Lévy processes under the forward measure. Proposition 7. Under the forward measure F with P(s , t1) as numeraire, the process W(j)F ss≥0 is the following: dW(j)F s=dW(j) t+hj(u,t1)du . (51) for j =1, . . . d is a Brownian motion. Furthermore, the random measure is the following: ˜ JF j(du,dz) = ˜ Jj(du,dz)−1|z|<ee−hj(u,t1)z−1νj(dz)du , (52) for j= 1, . . . d defines a jump process under ˜ P with the time-varying measure νF j(dt , dz) = e−hj(u,t1)zνj(dz). Proof. From Applebaum (2004, chap. 5, secs. 2 and 4), the forward change of the following measure: dF dQt1 =S(0) 0 S(0) t1P(0, t1)(53) =exp d ∑ j=1Zt1 0−hj(u,t1)σjdW(j) u−1 2hj(u,t1)2σ2 jdu − d ∑ j=1Zt1 0ZRe−hj(u,t1)z−1+hj(u,t1)z1|z|<eνj(dz)du + d ∑ j=1Zt1 0Z|z|≥e−hj(u,t1)z˜ Jj(du ,dz)! is a martingale under the measure Q , which defines an equivalent measure F under which W(j)F ss≥0 and ˜ JF j(du , dz) for j= 1, . . . d are, respectively, Brownian motions and random jump measures satisfying Equations (51) and (52). Combining Corollary 1and Proposition 7allows us to infer that a zero-coupon bond of maturity t≥t1is ruled under the forward measure by the SDE. dP(s,t) P(s,t)= rs+ d ∑ j=1 σjhj(s,t1)2+ d ∑ j=1ZRe−zhj(s,t1)−12νj(dz)!ds (54) + d ∑ j=1ZRe−zhj(s,t1)−1JF j(ds ,dz)−e−hj(s,t1)zνj(dz)ds − d ∑ j=1 σjhj(s,t1)dW(j)F s. This equation emphasizes that the average bond return under the forward measure is the risk-free rate plus a positive time varying premium that converges toward zero when s→t1 . The next proposition reveals that the processes Y(j,ξ) t do not revert anymore toward zero under F. Proposition 8. Let us define the following function. θ(j,ξ)(u) = (−hj(u,t1)σj+RRz1−e−hj(u,t1)z1|z|<eνj(dz)ξ>0 0ξ=0. (55) Risks 2022,10, 2 21 of 28 Under the forward measure F , Y(j,ξ) t is a process reverting to θ(j,ξ)(t)/ξ at a speed ξ such that the following is the case: Y(j,ξ) t=e−ξ(t−s)Y(j) s+Zt se−ξ(t−u)θ(j,ξ)(u)du (56) +σjZt se−ξ(t−u)dW(j)F u+Zt se−ξ(t−u)ZRz˜ JF j(du,dz), for all j =1, . . . , d and ξ∈R+. Proof. By definition of Y(j,ξ) s and from Corollary 1, we immediately infer its dynamic under the forward measure. dY(j,ξ) s=−ξY(j,ξ) sds +σjdW(j)F s+ZRz˜ JF j(ds,dz)(57) +µj−hj(s,t1)σj+Z|z|<eze−hj(s,t1)z−1νj(dz)ds . Since we impose that µj=−RRz−z1|z|<eνj(dz) for j= 1, . . . , d in order to ensure that EL(j) t=0, we deduce that Y(j,ξ) sreverts to θ(j,ξ)(t)/ξ: dY(j,ξ) s=ξ1 ξθ(j,ξ)(s)−Y(j,ξ) sds +σjdW(j)F s+ZRz˜ JF j(ds,dz), and the solution of this SDE is Equation (56). Equations (53) and (57) are useful for simulating sample paths of interest rates and bond prices under the forward measure, e.g., for pricing exotic options. One solution to price a call option on a zero-coupon bond consists in numerically inverting the Fourier’s transform of the density of its bond log-return under the forward measure. This step is performed with a standard discrete Fourier’s transform algorithm (DFT) that approximates the probability density function of the log-return under the forward measure F . The expectation of the payoff under the forward measure is next computed with this distribution. The moment generating function of the bond log-return, presented in the next proposition, admits a closed form formula and is a key result to implementing the DFT procedure. Proposition 9. The moment generating function (mgf) of ln P(t1 , t2) under measure F , conditionally to the filtration Fs, is given by the following. EFeωln P(t1,t2)|Fs=(58) exp − d ∑ j=1Zt1 sψj−hj(u,t1)du −ωZt2 t1 ϕ(u)du +ω d ∑ j=1Zt2 t1 ψj−hj(u,t2)du! ×exp −ω d ∑ j=1Z∞ 0Y(j,ξ) se−ξ(t1−s)B(ξ)(t1,t2)dγj(ξ)! exp d ∑ j=1Zt1 sψj−Z∞ 0 1 ξ1−(1+ω)e−ξ(t1−u)+ωe−ξ(t2−u)dγj(ξ)du!. Risks 2022,10, 2 22 of 28 Proof. The mgf can be rewritten as an expectation under the risk neutral measure using Bayes’ rule: EFeωln P(t1,t2)|Fs= EdF dQt1 eωln P(t1,t2)|Fs EdF dQt1 |Fs, (59) where dF dQt1 is detailed in Equation (53). From Equation (24) and due to the independence between the processes Y(j,ξ) t1and Y(k,ξ) t1for j6=k, expectation (59) becomes the following. EFeωln P(t1,t2)|Fs= exp − d ∑ j=1Zt1 sψj−hj(u,t1)du −ωZt2 t1 ϕ(u)du +ω d ∑ j=1Zt2 t1 ψj−hj(u,t2)du! × d ∏ j=1 Eexp−Zt1 shj(u,t1)dL(j) u−ωZ∞ 0Y(j,ξ) t1B(ξ)(t1,t2)dγj(ξ)| Fs. From Equation (18) and after a change of integration order, the integral of processes Y(j,ξ) t1is the following. Z∞ 0Y(j,ξ) t1B(ξ)(t1,t2)dγj(ξ) =Z∞ 0Y(j,ξ) se−ξ(t1−s)B(ξ)(t1,t2)dγj(ξ) + Z∞ 0Zt1 se−ξ(t1−u)dL(j) uB(ξ)(t1,t2)dγj(ξ) =Z∞ 0Y(j,ξ) se−ξ(t1−s)B(ξ)(t1,t2)dγj(ξ) + Zt1 sZ∞ 0e−ξ(t1−u)B(ξ)(t1,t2)dγj(ξ)dL(j) u. This allows us to rewrite EFeωln P(t1,t2)|Fsas follows, EFeωln P(t1,t2)|Fs= exp − d ∑ j=1Zt1 sψj−hj(u,t1)du −ωZt2 t1 ϕ(u)du +ω d ∑ j=1Zt2 t1 ψj−hj(u,t2)du! ×exp −ω d ∑ j=1Z∞ 0Y(j,ξ) se−ξ(t1−s)B(ξ)(t1,t2)dγj(ξ)! × d ∏ j=1 Eexp−Zt1 shj(u,t1)−ωZ∞ 0e−ξ(t1−u)B(ξ)(t1,t2)dγj(ξ)dL(j) u| Fs. By definition of hj(u,t):=R∞ 0B(ξ)(u,t)dγj(ξ), we have the following: hj(u,t1)−ωZ∞ 0e−ξ(t1−u)B(ξ)(t1,t2)dγj(ξ) = Z∞ 0B(ξ)(u,t1)−ωe−ξ(t1−u)B(ξ)(t1,t2)dγj(ξ), wherein the integrand is equal to the following. B(ξ)(u,t1)−ωe−ξ(t1−u)B(ξ)(t1,t2) =1 ξ1−(1+ω)e−ξ(t1−u)+ωe−ξ(t2−u). Risks 2022,10, 2 23 of 28 Combining these three last equations leads to the result. The expression (58) of the mgf involves integrals of Y(j,ξ) s with respect to ξ∈R+ . These integrals do not admit analytical formulas and are in practice approached with the discretization scheme introduced in Section 6. Let n be the number of atoms of the discrete approximation of γj for j= 1, . . . , d . To lighten future developments, we adopt the following notation: h(n) j(u,t):= n ∑ k=1 m(j) k b(j) k1−e−b(j) k(t−u), g(n) j(ω,u,t1,t2):= n ∑ k=1 m(j) k b(j) k1−(1+ω)e−b(j) k(t1−u)+ωe−b(j) k(t2−u), for j= 1, . . . , d . The discretized version of the mgf of the bond log-return under the forward measure is as follows: EFeωln P(n)(t1,t2)|Fs= exp d ∑ j=1Zt1 sψj−g(n) j(ω,u,t1,t2)−ψj−h(n) j(u,t1)du! ×exp Zt2 t1 ω d ∑ j=1 ψj−h(n) j(u,t2)−ω ϕ(n)(u)du! ×exp − d ∑ j=1 n ∑ k=1 m(j) kY(j,k) se−b(j) k(t1−s)B(j,k)(t1,t2)!, where ϕ(n)(t) , the discretized version of ϕ(t) is detailed in Equation (37). The integrals with respect to times are numerically computed, e.g., with a Simpson’s rule. Let us denote by f(n)(x)the probability density of ln P(n)(t1,t2)conditionally to Fsand the following. Υ(n)(iω) = EFeiωln P(n)(t1,t2)| Fs=Zeiωxf(n)(x)dx , The above is its Fourier’s transform. The probability density function can, therefore, be expressed as the real part of the inverse Fourier’s transform: f(n)(x) = 1 2πZ+∞ −∞ Υ(n)(iω)e−i x ωdω(60) =1 πReZ+∞ 0 Υ(n)(iω)e−i x ωdω, and is numerically computed with the algorithm detailed in the next proposition. Proposition 10. Let M be the number of steps used in the Discrete Fourier Transform (DFT) and ∆x=2xmax M−1be this step of discretization. Let us denote ∆ω=2π M∆xand the following: ωm= (m−1)∆ω, for m= 1 . . . M . The values of f(n)(x) and the pdf of ln P(n)(t1 , t2)|Fs at points xk=−M 2∆x+ (k−1)∆xare approached by the following sum: f(n)(xk)≈2 M∆x Re M ∑ m=1 $mΥ(n)(iωm)(−1)m−1e−i2π M(m−1)(k−1)!. (61) Risks 2022,10, 2 24 of 28 where $m=1 21{m=1}+1{m6=1}. The proof of this result is based on the trapezoidal approximation of the integral in Equation (60). Figure 9shows the probability density functions of P(n)( 5, 10 ) under the forward measure of numeraire P(n)(0, 5). We consider for this illustration a one factor Lévy model ( d= 1) in which the driving process is a jump-diffusion such as presented in Equation (48), with the same parameters as those of Section 6. The DFT parameters are xmax = 0.2 and M= 2 10 . The plots reveals that variances under F of bond prices display the same sensitivity to memory parameter α as under the risk of neutral measures. PML and ML bond variances are, respectively, directly and inversely proportional to α. 0.95 0.97 0.99 1.01 0 20 40 60 80 100 120 PML P(t1,t2) pdf alpha=0.9 alpha=0.5 0.95 0.97 0.99 1.01 0 20 40 60 80 ML P(t1,t2) pdf alpha=0.9 alpha=0.5 Figure 9. Probability density function of P(n)( 5, 10 ) under the forward measure of numeraire P(n)( 0, 5 ) , α= 0.5 and 0.9; β= 1.5; λ1 = 0.5; η1=− 0.002; σ1= 0.01. ML and PML kernels with n=40 atoms. After computation of the log-bond density, the call price is calculated by the following sum: C(s) = P(s,t1)EF(P(t1,t2)−K)+| Fs ≈P(s,t1) M ∑ k=1 1{xk≥ln K}(exk−K)f(n)(xk)∆x, where 1{xk≥ln K} is an indicator variable equal to one on [ln L , +∞) and zero otherwise. Figure 10 shows call prices on P(n)( 5, 10 ) for various strike prices in ML and PML kernels. In the PML model, increasing the memory parameter α drives up the option values with respect to whatever the strike price is. For the ML kernel, we observed the opposite trend. This is relevant with our previous conclusions about the variance that is the main driving factor of option prices. Figure 11 allows us to confirm the convergence of the call option when the number, n , of atoms in the discretization scheme grows. We recall that the partition ranges from ξ0= 0 up to ξn that is a percentile of γj(z) . Here, percentiles 90%, 95% and 97% are considered. For the ML kernel (right plot), the convergence is quick and the difference between prices computed with the 90% and 97% is negligible (around 2.6 10 −5 ). For the PML model (left plot), prices are higher (and then conservative) with a 90% percentile and converge with less atoms than with a partition covering 97% of γj(.).