Search tips
Search criteria 


Logo of nihpaAbout Author manuscriptsSubmit a manuscriptHHS Public Access; Author Manuscript; Accepted for publication in peer reviewed journal;
J Comput Phys. Author manuscript; available in PMC 2016 July 15.
Published in final edited form as:
J Comput Phys. 2015 July 15; 293: 14–28.
doi:  10.1016/
PMCID: PMC4465221



Fractional derivatives and integrals are convolutions with a power law. Multiplying by an exponential factor leads to tempered fractional derivatives and integrals. Tempered fractional diffusion equations, where the usual second derivative in space is replaced by a tempered fractional derivative, govern the limits of random walk models with an exponentially tempered power law jump distribution. The limiting tempered stable probability densities exhibit semi-heavy tails, which are commonly observed in finance. Tempered power law waiting times lead to tempered fractional time derivatives, which have proven useful in geophysics. The tempered fractional derivative or integral of a Brownian motion, called a tempered fractional Brownian motion, can exhibit semi-long range dependence. The increments of this process, called tempered fractional Gaussian noise, provide a useful new stochastic model for wind speed data. A tempered difference forms the basis for numerical methods to solve tempered fractional diffusion equations, and it also provides a useful new correlation model in time series.

Keywords: Fractional calculus, anomalous diffusion, random walk

1. Introduction

Fractional derivatives were invented by Leibnitz soon after the more familiar integer order derivatives [39, 54], but have only recently become popular in applications. They are now used to model a wide variety of problems in physics [25, 36, 39, 52, 53, 57, 68], finance [23, 28, 44, 38, 63, 64], biology [4, 2, 22, 27, 37], and hydrology [1, 7, 8, 15, 17, 65]. Fractional derivatives can be most easily understood in terms of their connection to probability [47, 48]. Einstein [20] explained the connection between random walks, Brownian motion, and the diffusion equation tp(x,t)=x2p(x,t). Sokolov and Klafter [66] review the modern theory of anomalous diffusion, where the integer order derivatives in the diffusion equation are replaced by their fractional analogues: tβp(x,t)=xαp(x,t). A fractional space derivative of order α < 2 corresponds to heavy tailed power law particle jumps P[J>x]xα (the famous Lévy flight), while a fractional time derivative of order β < 1 models heavy tailed power law waiting times P[W>t]tβ between jumps. Hence fractional space derivatives model anomalous super-diffusion, where a plume of particles spreads faster than the traditional diffusion equation predicts, and fractional time derivatives model anomalous sub-diffusion.

The goal of this paper is to describe a new variation on the fractional calculus, where power laws are tempered by an exponential factor. This exponential tempering has both mathematical and practical advantages. Mantegna and Stanley [41] proposed a truncated Lévy flight to capture the natural cutoff in real physical systems. Koponen [30] introduced the tempered Lévy flight as a smoother alternative, without a sharp cutoff. Cartea and del-Castillo-Negrete [12] developed the tempered fractional diffusion equation that governs the probability densities of the tempered Lévy flight. Unlike the truncated model, tempered Lévy flights offer a complete set of statistical physics and numerical analysis tools. Random walks with exponentially tempered power law jumps converge to a tempered stable motion [13]. Probability densities of the tempered stable motion solve a tempered fractional diffusion equation that describes the particle plume shape [3], just like the original Einstein model for traditional diffusion. Tempered fractional derivatives are approximated by tempered fractional difference quotients, and this facilitates finite difference schemes for solving tempered fractional diffusion equations [3]. The tempered diffusion model has already proven useful in applications to geophysics [46, 73, 74] and finance [10, 11]. In finance, the tempered stable process models price fluctuations with semi-heavy tails, resembling a pure power law at moderate time scales, but converging to a Gaussian at long time scales [5]. Since the anomalous diffusion eventually relaxes into a traditional diffusion profile at late time, this model is also called transient anomalous diffusion [74].

Kolmogorov [29] invented a new stochastic model for turbulence in the inertial range. Mandelbrot and Van Ness [40] pointed out that this stochastic process is the fractional derivative of a Brownian motion, and coined the name fractional Brownian motion. Fractional Brownian motion can exhibit a very useful property called long range dependence, where correlations fall off like a power law with time lag. A new variation on this model, called tempered fractional Brownian motion, is the tempered fractional derivative of a Brownian motion [49]. Tempered fractional Brownian motion can exhibit semi-long range dependence, with correlations that fall off like a power law at moderate time scales, but then eventually become short-range dependent at long time scales. This extends the Kolmogorov model for turbulence to also include low frequencies, and in fact tempered fractional Brownian motion provides a time-domain stochastic process model for the famous Davenport spectrum of wind speed [6, 16, 26, 55], which is used to design electric power generation facilities.

2. Tempered fractional diffusion

We begin by recalling the connection between random walks, Brownian motion, and the diffusion equation (see [48] for complete details). Given a random walk of mean zero particle jumps S(n) = X1 + . . . + Xn, the Central Limit Theorem implies that n−1/2S(nt) [implies] B(t) in distribution. The probability density p(x, t) of the Brownian motion limit B(t) solves the diffusion equation tp(x,t)=Dx2p(x,t). This useful connection between Brownian motion, random walks, and the diffusion equation assumes finite variance particle jumps. Power law jumps with density f(x) = Cαx−α−11[C1/α,∞)(t) for 1 < α < 2 have a finite mean but an infinite variance. Subtract the mean, and apply the extended central limit theorem [48, Theorem 3.37] to get n−1/αS(nt) [implies] A(t). Now the probability density p(x, t) of the α-stable Lévy motion A(t) solves the fractional diffusion equation tp(x,t)=Dxαp(x,t).

Tempered fractional diffusion applies an exponential tempering factor to the particle jump density. Consider a random walk Sε(n) with particle jump density


using the incomplete gamma function, and define the Poisson jump rate


for any ε > 0. To ease notation, we begin with the case of positive jumps.

Theorem 2.1. Suppose 0 < α < 1. Given a random walk Sε(n)=X1ε++Xnε with independent jumps, each having probability density function (1), and an independent Poisson process Ntε with rate (2), as ε → 0 we have the convergence


where the limit is a tempered stable process whose probability density function p(x, t) has Fourier transform


for any t ≥ 0.

Proof. The Poisson random variable satisfies P(Ntε=n)=eλet(λεt)nn! for n ≥0, and then


by the law of total probability. Apply the Fourier-Stieltjes transform f^(k)=eikxF(dx) where f(x) = F′(x), nothing that f^ε(k)n is the Fourier transform of probability distribution of Sε(n), to get


using the Taylor series for the complex exponential function. Rewrite this in the form


A straightforward computation using integration by parts and the formula for the gamma probability density [48, Eq. (3.14)] shows that


Use this formula twice to compute


when 0 < α < 1. Then as ε → 0, we get


for any t ≥ 0.

It follows from Theorem 2.1 that [p with hat](k,t) = etD[(λ+ik)α–λα] is the Fourier transform of the tempered Lévy flight A(t), a stochastic process with tempered power law jumps. We define the tempered fractional derivative xα,λf(x) as the function with Fourier transform [(λ+ik)αλα]f^(k) when 0 < α < 1. Note that


and invert the Fourier transform to see that the probability density p(x, t) of the tempered Lévy flight A(t) is the point source solution to the tempered fractional diffusion equation


We will also define the negative tempered fractional derivative xα,λf(x) as the function with Fourier transform [(λ − ik)α λα]f(k) when 0 < α < 1. A combination positive and negative jumps with particle jump density


for p, q ≥ 0 with p + q = 1 leads to a tempered stable limit process A(t) having both positive and negative jumps. A simple extension of Theorem 2.1 shows that (3) holds, and now the limit density p(x, t) solves the tempered fractional diffusion equation


To write the tempered fractional derivative in real space, use the formula (7) to see that


and then use the shift property of the Fourier transform ∫ e−ikxf(xy) dx = e−ikyf(k) to see that


for 0 < α < 1. A similar argument shows that


Theorem 2.2. Suppose 1 < α < 2. Given the mean-centered random walk Sε(n)=X1ε++Xnε with independent jumps, each having probability density function (1) and mean με=αCε1α(α1), and an independent Poisson process Ntε with rate (2), as ε → 0 we have the convergence


where the limit is a mean zero tempered stable Lévy motion whose probability density function p(x, t) has Fourier transform


for any t ≥ 0.

Proof. An easy computation shows that the jump variable Xε with probability density (1) has mean με=xfε(x)dx=αCε1α(α1). Then it follows from (5) that the density of Sε(Ntε)tλεμε has Fourier transform


for any t ≥ 0. Rewrite this in the form


using the fact that Γ(2 − α) = (1 − α)Γ(1 α). Use (6) and another integration by parts [48, p. 58] to see that


Let C = α(α − 1)/Γ(2 − α) and compute


using (13). Then as ε → 0, we get


for any t ≥ 0.

It follows from Theorem 2.2 that [p with hat](k, t) = etD[(λ+ik)α−λαikαλα−1] is the Fourier transform of the tempered Lévy flight A(t). We define the tempered fractional derivative xα,λf(x) as the function with Fourier transform [(λ + ik)α − λαikαλα−1]f(k) when 1 < α < 2. Since [partial differential]t[p with hat](k, t) = D[(λ + ik)α − λα − ikαλα−1][p with hat](k, t), we can invert the Fourier transform to see that the probability density p(x, t) of A(t) is the point source solution to the tempered fractional diffusion equation


Figure 1 shows a typical particle path, illustrating the effect of tempering to cool the long particle jumps.

Figure 1
Tempered fractional diffusion particle paths with with α = 1.2, from [3], showing the effect of tempering on long particle jumps.

To write this tempered fractional derivative in real space, use the formula (14) to see that


and invert the Fourier transform to get


for 1 < α < 2.

The negative tempered fractional derivative xα,λf(x) is the function with Fourier transform [(λ − ik)α −λα +ikαλα−1]f(k) when 1 < α < 2. In real space, we can write


The particle jump density (8) leads to a tempered stable limit process A(t) having both positive and negative jumps. A simple extension of Theorem 2.2 shows that (11) holds, and now the limit density p(x, t) solves the tempered fractional diffusion equation


If p = q then this produces a symmetric profile, with semi-heavy tails. That is, the tails of the probability density p(x, t) resemble a power law p(x, t) |x|α−1 for large |x| and moderate t > 0, but as t → ∞ the tails relax to a Gaussian profile. These semi- heavy tails are typical in applications to finance [10, 11]. Figure 2 fits a symmetric tempered stable density function to macroeconomic data on inflation rates. The data shows a sharper peak and a heavier tail than the best fitting Gaussian curve.

Figure 2
Symmetric tempered stable model fit (solid line) with α = 1.1, λ = 12, D = 0.1 to one-step bilinear ARMA forecast errors for annual inflation rates. The best fitting normal density (dotted line) fails to capture the sharp data peak.

The tempered space-fractional diffusion equation models transient super-diffusion, where a particle plume initially spreads faster than the traditional diffusion equation predicts, but later relaxes to a typical diffusion profile. The corresponding model for transient sub-diffusion assumes tempered power law waiting times between particle jumps. Suppose that the n th particle jump S(n) = X1 + · · · + Xn occurs at a random time Tn = W1 + · · · + Wn where Wn are all independent with density d(t) = Bβt−β−1e−ηt1[B1/β,)(t) for some 0 < β < 1. Theorem 2.1 yields random walk convergence to another stable Lévy motion Dt with Laplace transform g(s, t) = e−tψD(s) where the Laplace symbol ψD(s) = [(η + s)βηβ]. The number of jumps by time t > 0 is given by the renewal process Nt = max{n ≥ 0 : Tn ≤ t}, the inverse of Tn (graph Tn versus n, then swap the axes, to get the graph of Nt versus t). Then the CTRW (continuous time random walk) defined by S(Nt) gives the particle location at time t > 0.

A general result from CTRW limit theory [45, Theorem 2.1] shows that the CTRW converges at late time to A(Et), a tempered Lévy flight with the time variable replaced by an independent inverse tempered stable subodinator Et = inf{u > 0 : Du > t}. Since P[Etu]=P[Dut], we can take derivatives and then Laplace transforms to see that the density h(u, t) of Et has Laplace transform


see [45, Theorem 3.1] for details. Since x = A(u) has density function p(x, u) and u = Et has density function h(u, t), the CTRW limit process A(Et) has density


Recall that the tempered Lévy density has Fourier transform [p with hat](k, u) = e−uψA(k) with Fourier symbol ψA(k) = pD[(λ + ik)α − λα − ikαλα−1] + qD[(λ − ik)α − λα + ikαλα1] for 1 < α < 2. Take Laplace and Fourier transforms to see that


then rearrange and invert to see that the limit density solves the space-time tempered fractional diffusion equation


where the forcing term f(t) = [var phi]D(t,∞) has Laplace transform f(s) = s−1ψD(s), see [45, Theorem 4.1] for complete details. If η = 0, the boundary term reduces to f(t) = tβ−1/Γ(1−β), as in the fractional kinetic equation of Zaslavsky [72]. Tempered jumps are the basis for the popular CGMY model in finance [10, 11]. Tempered waiting times were applied to problems in ground water hydrology [46], and tempering in both variables was used for sediment transport in rivers [74]. Tempered fractional diffusion was first proposed by Cartea and del-Castillo-Negrete [12], and developed further in [3, 46], see also [48, Sections 7.2 and 7.3]. A more general tempering scheme was developed in [13]. The theory of tempered stable processes was developed by Rosiński [60], and an exact simulation scheme for the sample paths of a tempered Lévy flight was presented in Cohen and Rosiński [14].

3. Tempered fractional Calculus

Given any λ > 0, we define the positive tempered fractional integral of a suitable function f(x) by


and the negative tempered fractional integral by


If λ = 0, these formulae reduce to the well-known Riemann-Liouville fractional integrals [48, 61]. Hence we will call these operators the Riemann-Liouville tempered fractional integrals. Since these operators are convolutions, their Fourier transforms are products. Since g(x) = xα−1e−λx1[0,)(x)/Γ(α) has Fourier transform ĝ(k) = (λ + ik)−α, it follows from the convolution property of the Fourier transform that I+α,λf(x) has Fourier transform (λ + ik)−αf(k). Likewise Iα,λf(x) has Fourier transform (λ − ik)−αf(k), and when λ = 0 we recover the well-known formulae for the Fourier transform of a Riemann-Liouville fractional integral [61].

For functions in the fractional Sobolev space


we define the Riemann-Liouville tempered fractional derivatives D±α,λf(x) to be the functions with Fourier transform (λ ± ik)αf(k). Then integration and differentiation are inverse operators: D±α,λI±α,λf(x)=f(x) for any fL2(R), and I±α,λD±α,λf(x)=f(x) for any fWα,2(R). We can relate back to the normalized tempered fractional derivative ±xα,λf(x) from Section 2 with Fourier transform [(λ ± ik)α − λα]f(k) when 0 < α < 1, so that we can use (9) and (10) to write


Since ±xα,λf(x) has Fourier transform [(λ ± ik)α − λα [minus-or-plus sign] ikαλα−1]f(k) for 1 < α 2, we can also use (16) and (17) to write


Using the shift property of the Fourier transform, we can also relate the Riemann-Liouville tempered fractional derivative to the (untempered) fractional derivative xαf(x) from Section 2. Since ∫ e−ikxeλxf(x) dx = f(k + iλ), and since xαf(x) has Fourier transform (ik)αf(k) it follows that xα[eλxf(x)] has Fourier transform (ik)αf(k + iλ). Another application of the shift property shows that eλxxα[eλxf(x)] has Fourier transform (i(kiλ))αf((kiλ) + iλ) = (ik + λ)αf(k) and then the uniqueness of the Fourier transform implies that Note that the formula (37) remains valid when α is a positive integer. One can also use (37) to connect xαf(x) with the normalized tempered fractional derivative xα,λf(x), see [48, Section 7.2] for details.

D+α,λf(x)=eλxxα[eλxf(x)]for anyα>0.

Using the inner product f,g=f(x)g(x)¯dx on L2(R), along with the Plancherel Theorem left angle bracketf, gright angle bracket = left angle bracketf, ĝright angle bracket, it follows that


so that the positive and negative Riemann-Liouville tempered fractional derivatives are adjoint operators in the Hilbert space L2(R). Essentially the same argument shows that the positive and negative Riemann-Liouville tempered fractional integrals are also L2(R) adjoints. The semigroup property


follows by another easy Fourier transform argument.

The space of rapidly decreasing functions S(R) consists of the infinitely differentiable functions g:RR such that


where n, m are non-negative integers, and g(m) is the derivative of order m. The space of tempered distributions S(R) consists of continuous linear functionals on S(R). The Fourier transform maps S(R) into itself. If f:RR is polynomial growth, so that ∫ |f(x)|(1 + |x|)−pdx < ∞ for some p > 0, then Tf(ϕ) = ∫ f(x)ϕ(x)dx := left angle bracketf, ϕright angle bracket is a tempered distribution, also called a generalized function, with Fourier transform T^f(φ)=f^,φ=f,φ^=Tf(φ^) for φS(R). See Yosida [70, Ch.VI] for more details. For any tempered distribution f, the tempered fractional integral I+α,λf(x) exist as a convolution with the tempered distribution g(x) = xα−1e−λx1[0,∞)(x)/Γ(α). The same holds for Riemann-Liouville fractional integrals (the case λ = 0), since g(x) = xα11[0,∞)(x)/Γ(α) is of polynomial growth. Since G(k) = (λ ± ik)α is of polynomial growth, it is also a tempered distribution, and hence it follows from the Fourier inversion formula on S(R) that G = ĝ for some tempered distribution gS(R). Then we may define the Riemann-Liouville tempered fractional derivative of a tempered distribution as the convolution with this function. If n −1 < α < n, we can also write (λ + ik)αf(k) = (λ + ik)n(λ + ik)α−nf(k), and inverting yields an alternative definition in terms of the Riemann-Liouville tempered fractional integral:


where D+n,λf(x)=eλx[eλxf(x)](n) is a tempered derivative of integer order, defined by (23). There are many technical issues behind the fractional calculus of tempered distributions [61, Chapter 8], and it would be interesting to extend that theory for tempered fractional calculus.

4. Tempered fractional Brownian motion

As noted in Section 2, a random walk S(n) = X1 + · · · + Xn with mean zero, finite variance particle jumps converges to a Brownian motion B(t). For this Gaussian stochastic process, B(0) = 0, B(t + s) − B(s) has the same density function as B(t) (independent increment property) and B(t + s) − B(s) is independent of B(s) (independent increment property). The stochastic integral ∫ f(x)B(dx) is defined for any fL2(R) as a Gaussian random variable with mean zero and variance left angle bracketf, fright angle bracket = ∫ f(x)2dx, such the the Itô isometry holds: the stochastic integrals ∫ f(x)B(dx) and ∫ g(x)B(dx) have covariance left angle bracketf, gright angle bracket = ∫ f(x)g(x)dx. See [62] for more details. Then for any H > 0 and λ > 0 we can define the tempered fractional Brownian motion


where (x)+ = 1[0,∞)(x), and 00 = 0. The Gaussian stochastic process (24) has mean zero, and its increments are stationary, but not independent. Using the fact that B(c dx) [similar, equals] c1/2B(dx) (same distribution), it is not hard to check that tempered fractional Brownian motion has a nice scaling property: BH,λ(ct) [similar, equals] cHBH,cλ(t) for every c > 0. If λ = 0 and 0 < H < 1 (required to keep the integrand in L2(R) when λ = 0), we get the fractional Brownian motion with Hurst index H. Kolmogorov [29] proposed this process with H = 1/3 as a model for turbulence. Mandelbrot and Van Ness [40] established the connection between fractional Brownian motion and fractional calculus. Next we explain that connection in the tempered case.

Stochastic integrals with respect to Brownian motion can also be constructed using white noise theory. Heuristically, one writes ∫ f(x)B(dx) = f(x)W (x)dx where the white noise W (x) = B′ (x) is the (weak) derivative of the Brownian motion. This integral does not exist path-wise, since the paths of a Brownian motion are not differentiable. However, we can define W (f) := ∫ f(x)W (x)dx = left angle bracketf, Wright angle bracket as a random distribution (generalized function) on L2(R) by specifying that this object has the same (Gaussian) distribution and covariance function as ∫ f(x)B(dx) [31]. Taking f(x) = 1[0,t](x), it follows that B(t)=0tW(x)dx, and in this sense, the white noise W (x) is the (distributional) derivative of B(x). Suppose that H > 1/2, and compute


for any t > 0. Using the fact that positive and negative Riemann-Liouville tempered fractional integrals are L2(R) adjoints, we can therefore use white noise theory to write


which shows that tempered fractional Brownian motion with H > 1/2 can be obtained by integrating a linear combination of tempered fractional integrals of a white noise. The same is true for 0 < H < 1/2, if we interpret the tempered fractional integral of negative order as a tempered fractional derivative, see [51] for more details. Now if we take λ = 0 and 0 < H < 1, we recover the fact that fractional Brownian motion is the fractional integral of a white noise [40, 56].

Fractional Brownian motion with 1/2 < H < 1 is a popular model for long range dependence. The increments Xn = BH(n) − BH(n − 1) form a stationary time series of Gaussian random variables with mean zero, and the covariance between Xn and Xn+j falls off like |j|2H−2 as j → ∞ [62, Proposition 7.2.10]. Tempered fractional Brownian motion with a small tempering parameter λ ≈ 0 exhibits semi-long range dependence. The increments Xn = BH,λ(n) − BH,λ(n − 1) form a stationary Gaussian time series, the covariance between X and Xn+j falls off like |j|2H −2 for moderate values of j, but then tempering makes the covariance fall off faster (short range dependence) as j goes to infinity. When 0 < H < 1/2, both processes exhibit negative dependence, with a negative covariance between Xn and Xn+j for all j sufficiently large. See [49] for further details.

5. Tempered fractional difference operator

Tempered fractional derivatives can also be defined as the limit of a tempered fractional difference quotient.

Theorem 5.1. If f and its derivatives up to order n > 1 + α exist and are absolutely integrable, then the Riemann-Liouville tempered fractional derivative




Proof. The proof is similar to [48, Proposition 2.1]. The fractional binomial formula


holds for any α > 0 and any complex |z| ≤ 1. Since e−ikjhf(k) is the Fourier transform of f(xjh), the Fourier transform of hαΔλαf(x) is


as h → 0, using a Taylor expansion of the complex exponential function. Since f and its derivatives up to order n > 1 + α exist and are absolutely integrable, we have f(k) ≤ C(1 + |k|n) for some C > 0 [48, Lemma 2.4], and so (λ + ik)αf(k)| ≤ C2 + k2)α/2/(1 + |k|n). Hence (λ + ik)αf(k) is absolutely integrable, and then the inversion formula for the Fourier transform [48, Theorem 1.4] implies that there exists a function D+α,λf(x) with this Fourier transform. Then the continuity theorem for the Fourier transform [48, Theorem 1.3] implies that (25) holds.

The tempered fractional difference quotient can be applied to develop simple, effective numerical schemes to solve the tempered fractional diffusion equation. For physical reasons [7, 8, 17, 65], the most important case in applications is 1 < α < 2. Consider the tempered fractional diffusion equation with drift


where 1 < α ≤ 2, v(x) ≥ 0, and D(x) ≥ 0 on a finite domain axb, with initial condition p(x, 0) = p0(x), boundary conditions p(a, t) = 0 and p(b, t) = B(t), and forcing function f(x, t) for all t ≥ 0. Using Theorem 5.1 along with (22) and the fact that h−α(1 − e−λh)αf(x) → λαf(x) as h → 0, it follows that


To write the details of our numerical code, set ti = iΔt, xj = a + jΔx, h = Δx, pij = p(xj, ti), Dj = D(xj), vj = v(xj + Djαλα−1, and fij = f(xj, ti). It is natural to consider the implicit Euler scheme


using the weights


However, this implicit scheme is unstable. The proof is very similar to [43, Proposition 2.3]. The explicit Euler scheme (replace i + 1 by i on the right-hand side) is also unstable. The proof is essentially identical to [43, Proposition 2.1].

In order to obtain a stable convergent numerical method, it is necessary to use a shifted finite difference approximation


where the exponentially tempered Grünwald weights are defined by


The proof that (56) holds is similar to Theorem 5.1, see [3, Proposition 3] for details. Now consider the implicit Euler scheme


Theorem 5.2. The implicit Euler scheme (29) for the tempered fractional diffusion equation (26) using the weights (28) is consistent and unconditionally stable.

Proof. The proof is similar to [43, Theorem 2.7]. The main difference is that the Grünwald weights wj are replaced by the tempered weights gj. It follows from the fractional binomial formula that


and hence this numerical method mass-preserving. Since 1 < α < 2 we have w1 = −α and wj > 0 for all j ≠ 1, so that g1 < 0 and gj > 0 for all j ≠ 1. Then


for any finite number K. The implicit Euler scheme can be written as a linear system of equations APi+1 = Pi + FiΔt in matrix form where Pi = [pi,0, . . . , pi,N]′, Fi = [fi0, . . . , fiN]′, N = (ba)/h, and the inequality j=0Kgj<0 implies that every eigenvalue ζ of the implicit iteration matrix A satisfies |ζ| ≥ 1. Hence the spectral radius of the inverse matrix ρ(A1) ≤ 1, and the result follows.

Note that the shifted finite difference approximation (27) is essentially a convolution with the exponentially tempered Grünwald weights. An application of Stirling's formula shows that wm ~ −αm−α−1/Γ(1 − α) as m → ∞ [48, Eq. (2.5)], and this connects the finite difference scheme (29) with the probability model from Section 2. In the probability model, the tempered fractional derivative codes power law jumps with density function f(x) ≈ x−α−1e−λx for x large. In the finite difference approximation of the tempered fractional derivative, the particle flux into location x from location xjh follows the same tempered power law model. Between time ti and time ti + Δt, the mass at location xj is reduced by an amount proportional to the current particle mass density p(xj, ti) at this location, and increased by an amount proportional to the mass at other locations xj − mΔx with a weight g = wme−λΔx that falls off like a tempered power law with distance mΔx.

Example 5.3. Consider the tempered fractional diffusion equation (26) with v(x) = 0, x [set membership] [0,1], α = 1.5, λ = 0.5, initial condition u(x, 0) = xβe−λx/Γ(β + 1) with β = 2.5, diffusion coefficient c(x = xα(1 + βα)/Γ(β + 1), forcing function


and boundary conditions u(0, t) = 0 and u(1, t) = B(t) := e−λ−t/Γ(β + 1). The exact solution is given by u(x, t) = xβe−λx−t/Γ(1 + β). A similar example was considered in [3, Example 7]. Set [F with macron]i = fi1, . . . , fi,N−2, fi,N−1 + DN−1h g0B(ti+1)]′, pi0 = 0, pi,N = B(ti), and letĀ denote the (N − 1) × (N − 1) matrix with entries


Then the vector of interior points Pi = [pi,1, . . . , pi,N−1]′ solves the iteration equation ĀPi+1 = Pi + [F with macron]iΔt. The additional term in the last entry of [F with macron]i accounts for cutting off the super-diagonal entry of the full iteration matrix A. The algorithm was implemented in MATLAB (code is available from the authors upon request), see Table 1 for a summary of results.

Table 1
Comparison of maximum errors and error rate for Example 5.3 at time t = 1.0.

The implicit Euler method in Theorem 5.2 has error Ox + Δt), consistent with Table 1. A Crank-Nicolson method for solving the tempered fractional diffusion equation (26) with 1 < α ≤ 2 is developed in [3]. The method is stable and consistent, and second order convergent in time. However, because the finite difference approximation (27) is only first order accurate, the method is Ot2 + Δx). An application of Richardson extrapolation yields Ot2 + Δx2), see [3] for more details. It is also possible to construct particle tracking solutions based on the Monte Carlo simulation scheme in [3, Section 4], see [46, 74]. A wide variety of effective numerical methods have been developed to solve fractional diffusion equations, see for example [18, 33, 34, 35, 21, 59, 67, 71]. It would certainly be interesting to extend these methods to tempered fractional diffusion equations. Even in the untempered case, many interesting open problems remain, including a rigorous proof that fractional boundary value problems are well-posed, see [19] for a nice discussion.

6. ARTFIMA model

Another useful application of the fractional difference operator arises in time series analysis. Given a sequence of independent Gaussian random variables Zn with mean zero and variance σ2 (white noise sequence), a fractional noise sequence can be constructed by setting


where typically −1/2 < α < 1/2. Since ΔαΔβ = Δα+β in general (this follows from the fractional binomial formula), one can also write Zn = ΔαXn. Since the fractional difference is the discrete analogue of the fractional derivative, a fractional noise is closely related to fractional Brownian motion. Indeed, the sum X1 + · · · +XtBH(t) when H + α = 1/2. More precisely, we have


for all t > 0, see [69, Theorem 4.6.1]. The fractional noise can also exhibit long range dependence. The covariance between Xj and Xj+h falls off like |j|2H2 as j → ∞. In applications, a fractional difference is applied to a time series data set, and the order α of the fractional difference is chosen to remove long range correlations from the data. See [50] for a recent application to ground water hydrology.

Fractional differencing can be combined with the standard ARMA(p, q) model for time series, to yield the ARFIMA(p, α, q) time series model. The ARMA(p, q) model is


and we say that Xt follows an ARFIMA(p, α, q) if ΔαXt follows an ARMA(p, q) model. The spectral density fX(k) of a stationary time series Xt is the discrete Fourier transform of its covariance function. If γ(h) is the covariance between Xt and Xt+h, the spectral density


The spectral density of the noise sequence Zn is a constant fZ(k) = σ2/(2π), which is the reason for the term “white noise,” as all frequencies have equal weight. Using the backward shift operator notation BXt = Xt−1, one can write the ARMA(p, q) model in the compact form Φ(B)Xt = Θ(B)Zt where Φ(z) = 1 − [var phi]1z − · · · − [var phi]pzp and Θ(z) = 1 + θ1z + · · · θzq. Then it follows from the theory of linear filters that the spectral density of an ARMA(p, q) time series is given by


using the complex absolute value. Indeed, for any time invariant linear filter Xt = ΣjψjZt−j a standard argument shows that fX(k) = |ψ(e−ik)|2 fZ(k) where ψ(z) = Σjψjzj, see for example [9]. For an ARMA(p, q) time series, the filter function is ψ(z) = Θ(z)/Φ(z). Since


in operator notation, a similar argument shows that the spectral density of an ARFIMA(p, α, q) time series is given by


As k → 0, the fractional term |1 − eik|−2α ~ |k|−2α blows up at a power law rate. This spectral density is the hallmark of long range dependence.

The tempered fractional difference


can provide a more flexible model for (semi-)long range dependence. We say that Xt follows an ARTFIMA(p, α, λ, q) if ΔλαXt follows an ARMA(p, q) model. Since


in operator notation, it follows that the spectral density of an ARTFIMA(p, α, λ, q) time series is given by


For small values of the tempering parameter λ ≈ 0, this spectral density grows like a power law as k decreases, but then remains bounded as k 0 because of tempering. More details on ARTFIMA(p, α, λ, q) time series will be included in a forthcoming paper by the authors.

7. Turbulence

Kolmogorov [29] invented fractional Brownian motion as a model for turbulence in the inertial range. The main idea is that the spectral density grows like a power law as the frequency decreases, and a scaling argument leads to H = 1/3. This anti-persistent model has been verified in many real world experiments, for a range of frequencies. Davenport [16] modified this model in his study of wind speed at low altitudes, recognizing that the observed spectral density cannot actually diverge at low frequencies, outside the inertial range. The Davenport model for wind turbulence is st = μ + Xt where μ=E[st] is the average wind speed, and the stochastic process Xt has the normalized spectral density


where V10 is the mean velocity (m/sec) at an altitude of 10 meters, D is the corresponding drag coefficient, and x = 1200k/V10 where k is the wave number, see also Li and Kareem [32].

The spectral density of a stochastic process is the Fourier transform of its covariance function. The increment process Xt = BH,λ(t + 1) − BH,λ(t) of a tempered fractional Brownian motion has spectral density


for all real k, where C > 0 is some constant [49, Remark 4.3]. Since e−ik − ~ −ik as k → 0, this leads to the low frequency approximation


which shows that the stochastic process Xt has the Davenport spectrum with H = 5/6 and λ = V10/1200. The process Xt is called a tempered fractional Gaussian noise.

Figure 4 plots the spectral density of this tempered fractional Gaussian noise, along with the corresponding fractional Gaussian noise (the case λ = 0). Here we take λ = 0.01 to represent an average wind speed of 12 m/sec. The low frequency range corresponds to large eddy production. In the high frequency range, energy is dissipated via shear stress. In the intermediate (inertial) range, energy is transferred from smaller to larger eddies, and this is where the Kolmogorov model pertains. In that model, the spectral density is proportional to |k|1−2H for moderate frequencies. For more details, see for example Beaupuits et al. [6].

In the inertial range, it can be seen from Figure 4 that the spectrum of tempered fractional Gaussian noise is almost indistinguishable from that of (untempered) fractional Gaussian noise. However, there is a significant difference at low frequencies. The spectral density of fractional Gaussian noise diverges to infinity, which is not physical. However, the spectral density of tempered fractional Gaussian noise remains bounded, providing a more reasonable model for turbulence, consistent with the Davenport spectrum.

The Davenport model for wind turbulence is widely used to design windmill farms, offshore oil platforms, and large antennae [6, 16, 26, 55], but to date the model provides only a spectrum for wind speed data. Tempered fractional Gaussian noise, the increments of a tempered fractional Brownian motion, provides a stochastic process model in the time domain, which can also be useful in applications. For example, it can be used to simulate realistic wind speed data. It may also be useful to apply tempered fractional Brownian motion with H ≈ 1/3 as a model for turbulence in a more general setting. In the inertial range, it agrees well with the original Kolmogorov model (fractional Brownian motion), but it would seem to provide a more reasonable physical model at low frequencies.

The ARTFIMA process from Section 6 is the discrete time analogue of a tempered fractional Gaussian noise. Hence it could also provide a useful model for turbulence, when applied for example to a time series of turbulent velocities. From (30) we see that the ARTFIMA(0, α, λ, 0) spectral density is proportional to |1 − e(λ+ik)|−2α and so, for small values of the tempering parameter λ ≈ 0, the spectral density is approximately |k|1−2H for small frequencies k, recalling that α = 1/2 − H. As k → 0, the spectral remains bounded, because of the tempering.

8. Summary

Fractional calculus is a very useful analytical toolbox, with many diverse applications to all areas of science and engineering. This paper is promoting a new variation called tempered fractional calculus, as a more flexible alternative with considerable promise for practical applications. A fractional derivative (or integral) is a (distributional) convolution with a power law. A tempered fractional derivative multiplies that power law kernel by an exponential factor. Fractional diffusion equations govern random walks with a power law probability density for the particle jumps. Multiplying that power law by an exponential tempering factor leads to a tempered fractional diffusion equation. Point source solutions to those equations are tempered stable probability densities, with semi-heavy tails that transition from power law to Gaussian. The emerging theory of tempered fractional calculus provides a sound mathematical basis for applications. Tempered fractional Brownian motion, the tempered fractional derivative or integral of a Brownian motion, extends the popular fractional Brownian motion, and provides a flexible model for semi-long range dependence that transitions from long to short range correlations. The tempered fractional difference operator can provide the basis for effective numerical schemes to solve tempered fractional differential equations. It is also useful as a model for time series analysis, where ARTFIMA(p, α, λ, q) models incorporate semi-long range dependence in a natural way. Tempered fractional Brownian motion provides a new stochastic process model for turbulence, consistent with the Davenport spectrum that extends the Kolmogorov theory of turbulence beyond the inertial range, and the ARTFIMA(p, α, λ, q) model provides a useful discrete time analogue.

Figure 3
Spectral density of tempered fractional Gaussian noise with H = 5/6 and λ = 0.01 (solid line) and of the corresponding fractional Gaussian noise with H = 5/6 (dotted line).


This work was partially supported by NSF grant DMS-1025486, a Fujian Governmental Scholarship, and the National Natural Science Foundation of China (Grant No. 11101344).


Publisher's Disclaimer: This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

Contributor Information

MARK M. MEERSCHAERT, Department of Statistics and Probability, Michigan State University, East Lansing MI 48823 ; ude.usm.tts@debucm URL:

FARZAD SABZIKAR, Department of Statistics and Probability, Michigan State University, East Lansing MI 48823 ; ude.usm.tts@2akizbas.

JINGHUA CHEN, School of Sciences, Jimei University, Xiamen, China, 361500 ; moc.361@zdzdhjc..


1. Baeumer B, Meerschaert MM, Benson DA, Wheatcraft SW. Subordinated advection-dispersion equation for contaminant transport. Water Resour. Res. 2001;37:1543–1550.
2. Baeumer B, Kovács M, Meerschaert MM. Fractional reproduction-dispersal equations and heavy tail dispersal kernels. Bull. Math. Biol. 2007;69:2281–2297. [PubMed]
3. Baeumer B, Meerschaert MM. Tempered stable Lévy motion and transient super-diffusion. J. Comput. Appl. Math. 2010;233:243–2448.
4. Barkai E, Garini Y, Metzler R. Strange kinetics of single molecules in living cells. Phys. Today. 2012;65:29.
5. Barndorff-Nielsen OE. Processes of normal inverse Gaussian type. Finance Stoch. 1998;2:41–68.
6. Pérez Beaupuits JP, Otárola A, Rantakyrö FT, Rivera RC, Radford SJE, Nyman L-Å. Analysis of wind data gathered at Chajnantor. ALMA Memo. 2004;497
7. Benson D, Wheatcraft S, Meerschaert M. Application of a fractional advection-dispersion equation. Water Resour. Res. 2000;36:1403–1412.
8. Benson D, Schumer R, Meerschaert M, Wheatcraft S. Fractional dispersion, Lévy motions, and the MADE tracer tests. Transport in Porous Media. 2001;42:211–240.
9. Brockwell PJ, Davis RA. Time Series: Theory and Methods. 2nd Ed. Springer-Verlag; New York: 1991.
10. Carr P, Geman H, Madan DB, Yor M. The fine structure of asset returns: An empirical investigation. J. Business. 2002;75:303–325.
11. Carr P, Geman H, Madan DB, Yor M. Stochastic volatility for Lévy processes. Math. Finance. 2003;13:345–382.
12. Cartea Á, del-Castillo-Negrete D. Fluid limit of the continuous-time random walk with general Lévy jump distribution functions. Phys. Rev. E. 2007;76:041105. [PubMed]
13. Chakrabarty A, Meerschaert MM. Tempered stable laws as random walk limits. Statist. Probab. Lett. 2011;81:989–997.
14. Cohen S, Rosiński J. Gaussian approximation of multivariate Lévy processes with applications to simulation of tempered stable processes. Bernoulli. 2007;13:195–210.
15. Cushman JH, Ginn TR. Fractional advection-dispersion equation: A classical mass balance with convolution-Fickian flux. Water Resour. Res. 2000;36:3763–3766.
16. Davenport AG. The spectrum of horizontal gustiness near the ground in high winds. Q. J. Royal Meteor. Soc. 1961;87:194–211.
17. Deng Z-Q, Bengtsson L, Singh VP. Parameter estimation for fractional dispersion model for rivers. Environ. Fluid Mech. 2006;6:451–475.
18. Diethelm K, Ford NJ, Freed AD, Luchko Yu. Algorithms for the fractional calculus: A selection of numerical methods. Comput. Methods Appl. Mech. Engrg. 2005;194:743–773.
19. Du Q, Gunzburger M, Lehoucq RB, Zhou K. Analysis and approximation of nonlocal diffusion problems with volume constraints. SIAM Review. 2012;54:667–696.
20. Einstein A. On the movement of small particles suspended in a stationary liquid demanded by the molecular kinetic theory of heat. Ann. Phys. 1905;17:549–560.
21. Ervin VJ, Roop JP. Variational solution of fractional advection dispersion equations on bounded domains in Rd. Numer. Meth. P.D.E. 2007;23:256–281.
22. Fedotov S, Iomin A. Migration and proliferation dichotomy in tumor-cell invasion. Phys Rev Lett. 2007;98:8101. [PubMed]
23. Gorenflo R, Mainardi F, Scalas E, Raberto M. Trends Math. Birkhäuser; Basel: 2001. Fractional calculus and continuous-time finance. III. The diffusion limit. Mathematical finance (Konstanz, 2000) pp. 171–180.
24. Herrmann R. Fractional Calculus an Introduction for Physicists. World Scientific; Singapore: 2011.
25. Hilfer R. Applications of Fractional Calculus in Physics. World Scientific; 2000.
26. Jang J-J, Jyh-Shinn G. Analysis of maximum wind force for offshore structure design. J. Marine Sci. Tech. 1999;7:43–51.
27. Jeon J-H, Martinez-Seara Monne H, Javanainen M, Metzler R. Anomalous diffusion of phospholipids and cholesterols in a lipid bilayer and its origins. Phys. Rev. Lett. 2012;109:188103. [PubMed]
28. Jurlewicz A, Wy lomańska A, Żebrowski P. Coupled continuous-time random walk approach to the Rachev-Rüschendorf model for financial data. Physica A. 2009;388:407–418.
29. Kolmogorov AN. Wiener spiral and some other interesting curves in Hilbert space. Dokl. Akad. Nauk SSSR. 1940;26:115–118.
30. Koponen I. Analytic approach to the problem of convergence of truncated Lévy flights towards the Gaussian stochastic process. Phys. Rev. E. 1995;52:1197–1199. [PubMed]
31. Kuo HH. White noise distribution theory. CRC Press; Boca Raton, Florida: 1996.
32. Li Y, Kareem A. ARMA systems in wind engineering. Probabilistic Engin. Mech. 1990;5:49–59.
33. Liu F, Ahn V, Turner I, Zhuang P. Numerical simulation for solute transport in fractal porous media. ANZIAM J. 2004;45:C461–C473.
34. Liu F, Ahn V, Turner I. Numerical solution of the space fractional Fokker-Planck equation. J. Comput. Appl. Math. 2004;166:209–219.
35. Lynch VE, Carreras BA, del-Castillo-Negrete D, Ferreira-Mejias KM, Hicks HR. Numerical methods for the solution of partial differential equations of fractional order. J. Comput. Phys. 2003;192:406–421.
36. Magdziarz M, Weron A, Weron K. Fractional Fokker-Planck dynamics: Stochastic representation and computer simulation. Phys. Rev. E. 2007;75:016708. [PubMed]
37. Magin RL. Fractional Calculus in Bioengineering. Begell House; 2006.
38. Mainardi F, Raberto M, Gorenflo R, Scalas E. Fractional Calculus and continuous-time finance II: the waiting-time distribution. Physica A. 2000;287:468–481.
39. Mainardi F. Fractional Calculus and Waves in Linear Viscoelasticity: An Introduction to Mathematical Models. World Scientific; 2010.
40. Mandelbrot B, Van Ness J. Fractional Brownian motion, fractional noises and applications. SIAM Review. 1968;10:422–437.
41. Mantegna RN, Stanley HE. Stochastic process with ultraslow convergence to a Gaussian: The truncated Lévy flight. Phys. Rev. Lett. 1994;73:2946–2949. [PubMed]
42. Meerschaert MM, Scheffler H-P. Limit theorems for continuous time random walks with infinite mean waiting times. J. Applied Probab. 2004;41:623–638.
43. Meerschaert MM, Tadjeran C. Finite difference approximations for fractional advection-dispersion flow equations. J. Comput. Appl. Math. 2004;172:65–77.
44. Meerschaert MM, Scalas E. Coupled continuous time random walks in finance. Physica A. 2006;370:114–118.
45. Meerschaert MM, Scheffler H-P. Triangular array limits for continuous time random walks. Stoch. Proc. Appl. 2008;118:1606–1633.
46. Meerschaert MM, Zhang Y, Baeumer B. Tempered anomalous diffusions in heterogeneous systems. Geophys. Res. Lett. 2008;35:L17403–L17407.
47. Meerschaert MM. In: Fractional calculus, anomalous diffusion, and probability, Fractional Dynamics. Metzler R, Klafter J, editors. World Scientific; Singapore: 2012. pp. 265–284.
48. Meerschaert MM, Sikorskii A. Stochastic Models for Fractional Calculus. De Gruyter; Berlin: 2012.
49. Meerschaert MM, Sabzikar F. Tempered fractional Brownian motion. Statist. Prob. Lett. 2013;83:2269–2275.
50. Meerschaert MM, Dogan M, Van Dam RL, Hyndman DW, Benson DA. Hydraulic conductivity fields: Gaussian or not? Water Resour. Res. 2013;49:4730–4737. [PMC free article] [PubMed]
51. Meerschaert MM, Sabzikar F. Stochastic integration for tempered fractional Brownian motion. Stoch. Proc. Appl. 2014;124:2363–2387. [PMC free article] [PubMed]
52. Metzler R, Klafter J. The random walk's guide to anomalous diffusion: A fractional dynamics approach. Phys. Rep. 2000;339:1–77.
53. Metzler R, Klafter J. The restaurant at the end of the random walk: recent developments in the description of anomalous transport by fractional dynamics. J. Physics A. 2004;37:R161–R208.
54. Miller K, Ross B. An Introduction to the Fractional Calculus and Fractional Differential Equations. Wiley and Sons; New York: 1993.
55. Norton DJ. Mobile offshore platform wind loads. Proc. 13th Offshore Techn. Conf., OTC 4123. 1981;4:77–88.
56. Pipiras V, Taqqu M. Integration questions related to fractional Brownian motion. Prob. Th. Rel. Fields. 2000;118:251–291.
57. Piryatinska A, Saichev AI, Woyczynski WA. Models of anomalous diffusion: Subdiffusive and superdiffusive cases. Physica A. 2005;349:375–424.
58. Ramos-Fernandez G, Mateos JL, Miramontes O, Cocho G, Larralde H, Ayala-Orozco B. Lévy walk patterns in the foraging movements of spider monkeys (Ateles geoffroyi). Behav. Ecol. Sociobiol. 2003;55:223–230.
59. Roop JP. Computational aspects of FEM approximation of fractional advection dispersion equations on bounded domains in R2. J. Comput. Appl. Math. 2006;193:243–268.
60. Rosiński J. Tempering stable processes. Stoch. Proc. Appl. 2007;117:677–707.
61. Samko SG, Kilbas AA, Marichev OI. Fractional Integrals and Derivatives. Gordon and Breach. 1993
62. Samorodnitsky G, Taqqu M. Stable non-Gaussian Random Processes. Chapman and Hall; New York: 1994.
63. Scalas E, Gorenflo R, Mainardi F. Fractional calculus and continuous-time finance. Physica A. 2000;284:376–384.
64. Scalas E. The Complex Networks of Economic Interactions. Springer; Berlin: 2006. Five years of continuous-time random walks in econophysics. pp. 3–16.
65. Schumer R, Benson DA, Meerschaert MM, Wheatcraft SW. Eulerian derivation of the fractional advection-dispersion equation. J. Contaminant Hydrol. 2001;48:69–88. [PubMed]
66. Sokolov IM, Klafter J. Anomalous diffusion spreads its wings. Physics World. 2005;18:29–32.
67. Tadjeran C, Meerschaert MM, Scheffler H-P. A second order accurate numerical approximation for the fractional diffusion equation. J. Comput. Phys. 2006;213:205–213.
68. Tarasov VE. Fractional vector calculus and fractional Maxwell's equations. Ann. Phys. 2008;323:2756–2778.
69. Whitt W. Stochastic-Process Limits. Springer; New York: 2002.
70. Yosida K. Functional Analysis. Sixth edition Springer; 1980.
71. Yuste SB, Acedo L. An explicit finite difference method and a new von Neumann type stability analysis for fractional diffusion equations. SIAM J. Numer. Anal. 2005;42:1862–1874.
72. Zaslavsky G. Fractional kinetic equation for Hamiltonian chaos. Physica D. 1994;76:110–122.
73. Zhang Y, Meerschaert MM. Gaussian setting time for solute transport in fluvial systems. Water Resour. Res. 2011;47:W08601.
74. Zhang Y, Meerschaert MM, Packman AI. Linking fluvial bed sediment transport across scales. Geophys. Res. Lett. 2012;39:L20404.