FRACTIONAL DIFFERENTIAL EQUATIONS DRIVEN BY LÉVY NOISE

This paper considers a general class of fractional differential equations driven by Lévy noise. The singularity spectrum for these equations is obtained. This result allows to determine the conditions under which the solution is a semimartingale. The prediction formula and a numerical scheme for approximating the sample paths of these equations are given. Almost sure, uniform convergence of the scheme and some numerical results are also provided.

Barndorff-Nielsen [7,6] used discrete or continuous-type superposition of Ornstein-Uhlenbeck processes with Lévy motion input to obtain a class of random processes with LRD and infinitely divisible marginal distributions, while Igloi and Terdik [25] and Oppenheim and Viano [36] obtained long memory by aggregating continuous-time short-memory Gaussian processes with random coefficients.
In a continuous-time framework, it is known that LRD can be obtained by replacing ordinary differential operators by fractional differential operators in differential equations driven by white noise (see, for example, Gay and Heyde [23], Inoue [26], Viano et al. [42], Chambers [14], Anh et al. [1], Anh and Leonenko [3]).Following this approach, Anh et al. [2] introduced a class of general fractional differential equations (FDE) driven by Lévy noise, whose solutions are obtained as convolutions of the Green functions of the corresponding deterministic FDEs with Lévy noise or stochastic path integrals with respect to Lévy processes.The main advantage of this approach is that LRD can be effected via the Green function of the fractional operator involved, while the noise term can be used to represent the effects of non-Gaussianity or multifractality.Anh et al. [2] obtained explicit results on the Green functions, correlation functions, spectra and higher-order spectra of particular forms of these FDEs.In particular, they demonstrated that these equations can be used to model the stochastic volatility of log price processes and macroeconomic processes with long memory.
In this paper, we continue this line of investigation.Some elements of the above class of FDEs, as presented in Anh et al. [2], will be recalled in the next section.Their multifractality, which is described by the corresponding singularity spectrum, is obtained in Section 3.This result provides not only some sample path properties of the FDEs but also an important device to determine the conditions under which the solution is a semimartingale.This semimartingale representation is derived in Section 4. Section 5 gives a closed-form prediction formula, while Section 6 develops a numerical scheme for approximating the sample paths of these equations.Almost sure, uniform convergence of the scheme and some numerical results will also be given.

Fractional Differential Equations Driven by Lévy Noise
Let us first recall the definitions of fractional derivative and fractional integral (see Samko et al. [40], Djrbashian [18], Podlubny [37] among others).Assuming reasonable behaviour for f (t) , its Riemann-Liouville fractional derivative is defined as .., and its Riemann-Liouville fractional integral is defined as (2.2) In this paper, we consider fractional differential equations of the general form ) where L is Lévy noise.Note that Lévy noise L has the following properties: (1) it is infinitely divisible; (2) its probability distribution is translation invariant and (3) L (t) and L (s) are independent if t = s (see Mueller [33], for example).
As defined in Podlubny [37], p.150, the function G (t − τ ) satisfying the following conditions a) is called the Green function of Eq. (2.3).Its solution in terms of the Green function can be written as where formally L (t) = t 0 L (s) ds, and the stochastic integral is known to exist (Anh et al. [2]).Alternatively, the integral (2.5) may be interpreted pathwise.Let v p (f ) be the p−variation of f ; then If G is of bounded q−variation and L is of bounded p−variation with q −1 + p −1 > 1 then the integral may be interpreted • in the Riemann-Stieltjes sense whenever G and paths of L have no discontinuities at the same points; • in the Moore-Pollard-Stieltjes sense whenever G and paths of L have no one-sided discontinuities at the same points; • always in the sense defined by Young [43] (see Dudley and Norvaiša [19] for details).
When the integral exists, the Love-Young inequality (Theorem 2.1 of Dudley and Norvaiša [19]) holds, where We should note that the Green function G (t) , t ≥ 0, is usually defined by its Laplace transform (2.7) Consider the fractional integral equation We now use the Laplace transform method to obtain an explicit representation of the solution to (2.8).Let L T (t) = L (t ∧ T ).The Laplace transform of L T is well defined and will be denoted by L T (p).As L T (t) is bounded for all t ∈ [0, ∞), X (t) is also bounded on this interval and hence its Laplace transform X (p) is well defined.Applying the Laplace transform to (2.8) we have Inverting the Laplace transform gives (2.12) The integration by parts formula (see Dudley and Norvaiša [19], p.67) shows the equivalence of (2.5) and (2.12) for t < T .We will use the representation (2.8) to study the path properties of the solution of fractional differential equations when β n ≥ 1.For β n < 1 we will interpret the solution to (2.3) as , where Y (t) solves We next recall some basic definitions on the local scaling properties of the paths of a process X (t) on some interval [0, T ]; for futher details see, for example, Jaffard [27] or Riedi [38] .A typical feature of a fractal process X (t) is that it has a non-integer degree of differentiability, characterised by its local Hölder exponent h (t) defined by for t sufficiently close to t and P t (.) being the Taylor polynomial of X at t.The sets which form a decomposition of the support of X according to its singularity exponents, can be highly interwoven and dense on [0, T ].The singularity spectrum of X is then defined as d (h) = dim (E h ) where dim is the Hausdorff dimension.A process X is said to be multifractal if the support of its singularity spectrum has a non-empty interior.A classical example of a multifractal process is the multiplicative cascade on the interval [0, 1] (Mandelbrot [32]).However, such multiplicative cascades are not suitable for the stochastic models considered here as they do not possess stationary increments and are only defined on some finite interval.An example of a stochastic process with stationary increments and defined on [0, ∞) which is also a multifractal is Lévy motion.Jaffard [28] showed that all Lévy motions are multifractal with the exception of Brownian motion, compound Poisson processes, deterministic motion and their convolutions.The singularity spectrum of a Lévy motion without Brownian component was shown to be where γ is given by and is called the upper index of the Lévy measure.

Multifractality of Fractional Differential Equations
We first give an extension of the mapping property of fractional integrals on the space ) by the fractional integral of order α (see Samko et al. [40], Theorem 3.1).

Lemma 1 Let X (t) be a process with bounded sample paths. Denote the local Hölder exponents of X (t) by h X (t). Define Y (t) by
with α ∈ (0, 1).Then the following inequality holds: Proof For some t > t (the case of t < t can be dealt with in a similar fashion), Adding and subtracting the terms where ) We first show that J 3 is a local polynomial for t > 0. Now, For |t − t| sufficiently small and t > 0, J 3 is a local polynomial and so will not affect the local Hölder exponents.From the definition of local Hölder exponents for s in a neighbourhood of t, for any ε > 0. Applying this to J 1 , we get For J 2 , the interval of integration is divided into [2t − t , t] and [0, 2t − t ].Applying (3.9) to J 2 on the first interval of integration yields For the second inteval of integration, (3.18)where Applying Fubini's theorem, the order of integration and summation can be changed to yield As α + h X (t) > 1 we may choose an ε > 0 small enough so that α + h X (t) − ε > 1 and the first integral on the left hand side is finite.The second integral is of the order where The second integral is of the order The equations (3.21) and (3.27) yield where C is a constant independent of t .Combining (3.8), (3.11), (3.14) and (3.28) we have where P t (t ) is a local polynomial of Y (t) and for any ε > 0. Letting ε → 0 we see that Theorem 1 Let X (t) be the solution to the fractional integral equation (2.8) with β n ≥ 1 and let L (t) a process with bounded sample paths.The singularity spectrum of Proof From our previous discussion on the solution of (2.8) we can say that, if L (t) has bounded sample paths, then X (t) will also have bounded sample paths.Clearly the local Hölder exponents of I β n −1 L (t) and coincide by (2.8).As the local Hölder exponent of the sum of two functions is the infimum of the two, except perhaps when they are equal, in which case it may be larger, we may apply Lemma 1 to conclude that and Theorem 1 follows.Remark 1 A specific example of Theorem 1 is obtained when we let L (t) be a Lévy motion.As the sample paths of all Lévy motions are right continuous with left limits, almost surely, they must be bounded on any finite interval [0, T ].When β n = 1 and 1/γ < 2 − β n−1 , where γ is the upper index of the Lévy measure and L (t) has no Brownian component, then Theorem 1 implies that the singularity spectrum of X (t) is given by (2.16).
When β n > 1 the right hand side of (2.8) is the fractional integral of Lévy motion for which we do not know the singularity spectrum.The following lemma provides an upper bound on the singularity spectrum.

Lemma 2 Let L (t) be a Lévy motion with upper index γ ≥ 1 and no Brownian component. The singularity spectrum of the fractional integral of Lévy motion is bounded by
Proof The proof follows from Lemma 1 which gives h I β L (t) ≥ β + h L (t), the multifractal formalism (see Riedi [38], Section 7.1) and the form (2.16) of the singularity spectrum of Lévy motion.

Semimartingale Representation
An interesting problem in the development of a stochastic calculus for fractional differential equations is to determine when the solution has the semimartingale representation.Fractional Brownian motion, which is characterized by a spectral density of the form is known to be a semimartingale if and only if γ = 1 (Lin [31]).Another dynamic model, fractional Riesz-Bessel motion, which is characterised by a spectral density of the form is known to be a semimartingale when α+γ ↓ 3/2, (Anh and Nguyen [4]).Using Lemma 1, appropriate conditions on a fractional differential equation for its solution to be a semimartingale can be given.Theorem 2 Let X (t) be the solution to (2.3) with L (t) being a Lévy motion whose Lévy measure has upper index γ and without Brownian component.Then the following results hold: is unbounded on any finite interval and thus not a semimartingale.
, where Y (t) is a process with a dense set of discontinuities, almost surely.It follows that X (t) is unbounded on any finite interval and hence X (t) is not a semimartingale.
(ii) For any sequence of partitions From Anh et al. [2] it is known that as ω → ∞.It follows that for any ε > 0, A function will be of bounded 1-variation if and only if it is differentiable almost everywhere.Therefore, if the local Hölder exponent of X (t) is greater than one almost everywhere, it will be of bounded 1-variation.From Lemma 1 the local Hölder exponent of X (t) will be greater than one almost everywhere if is of bounded 1-variation.Substituting X (t) = −B (t) + L (t) into (4.8) and applying Lemma 1, it is seen that this is equivalent to showing that I 1−β n−1 L (t) is of bounded 1-variation.Using the same argument as in (ii) we conclude that B (t) is of bounded 1-variation if β n−1 < 1/γ, and hence X (t) is a semimartingale.Remark 2 We can allow the Lévy motion in Theorem 2 to possess a Brownian component, in which case we set γ = 2 regardless of the Lévy measure.Parts (ii) and (iii) of the theorem hold with the same proof if a Brownian component is included.For part (i) we note that G (t) ∼ Ct β n −1 as t → 0 (see Anh et al. [2]).Applying Proposition 13 of Carmona et al. [13] shows that X (t) will have no quadratic variation if β n < 1 and hence X (t) is not a semimartingale for β n < 1.
Remark 3 As in the proof of Theorem 1, we do not make use of the assumption that L (t) is a Lévy motion, only that it is a semimartingale with Hölder exponent 1/γ almost everywhere.As all semimartingales have finite 2−variation, Lemma 4.3 of Dudley and Norvaiša [19] implies that semimartingales have a Hölder exponent greater than or equal to 1/2 almost everywhere.Following the arguments of the proof of Theorem 2 (iii) we see that X (t) will be a semimartingale if β n = 1, β n−1 < 1/2 and L (t) is a semimartingale.

Prediction Formula
Consider first the simple case of the fractional differential equation (5.1) which was studied in Anh et al. [2].It is assumed that we have observed the sample path of X (t) from its initial condition X (0) = 0 and so F t = σ (X (s) : 0 ≤ s ≤ t).
To compute the conditional expectation E (X (t + u) |F t ) we can take advantage of the form of the solution as follows: (5.3) A simple change of variable yields Theorem 1.2 of Djrbashian [18] gives the solution to the above intregral equation as (5.5) For the more general form of fractional differential equation (2.3) with β n = 1, the prediction formula is given by the solution to the integral equation where whose solution may be written as (5.7) G (t) being the Green function for the fractional differential equation (5.8)

Numerical Approximation of Sample Paths
This section will provide a numerical scheme to approximate the sample paths of equations of the form (2.3).We first show that their Green function is completely monotonic on [0, ∞) if β n = 1.Bernstein's theorem then implies that this Green function can be written in the form where µ (dλ) is a finite Borel measure on [0, ∞) (Feller [22]).In view of (2.5) and (6.1), the solution of (2.3) with β n = 1 is then given by 2) It should be noted that t 0 e −λ(t−s) dL (s) is the solution of the Orstein-Uhlenbeck-type equation dY (λ, t) = −λY (λ, t) dt + dL (t) , Y (λ, 0) = 0. (6.3) We will follow Carmona et al. [13] and Camona and Coutin [12] which give a fast and efficient algorithm for simulating Gaussian LRD processes based on the representation where µ (dλ) is a Borel measure satisfying and B (t) is a Brownian motion.They prove mean-square and almost-sure convergence in this case.An example of a process with representation (6.4) is a type II fractional Brownian motion.When (6.1) is square integrable on [0, ∞) the process (6.4) is not stationary, but does converge in mean square to a stationary Gaussian process.Furthermore, (6.4) can be made stationary by letting the initial condition in (6.3) be a Gaussian random function on [0, ∞) with mean zero and covariance function We will show that the solution to a certain class of FDEs has the representation (6.2).
In this case we extend the almost sure convergence result to include Lévy motion as a driving term.

Approximation algorithm
The following approximation scheme is adapted from Carmona et al. [13].Define a compact set K ⊂ [0, ∞) by K = [r −m , r n ] with m, n being positive integers and r > 1. Denote the geometric partition of K by π = {A i } with A i = r i , r i+1 , i = −m, . . ., n − 1.Consider the following approximation of X (t) : where Y r i , t is the solution to the Ornstein-Uhlenbeck-type equation subject to the initial condition Y (λ, 0) = 0. Finally, the Ornstein-Uhlenbeck-type process is approximated by for 0 < t ≤ ∆ and for (n − 1) ∆ < t ≤ n∆.The approximation of X (t) is then

Representation as sums of OU-type processes
For the approximation (6.11) to be useful in simulating sample paths of FDEs, we need to determine if the Green function has the representation (6.1).While it is clear that this cannot be the case for β n > 1 as G (0) = 0, the Green function will have the representation (6.1) if β n = 1.For this we need part (iii) of Theorem 5.10 from Inoue [26] (see also its Theorems 1.1 and 1.2).Theorem 3 Let R (t) be the covariance function of a stationary process satisfying where σ (dλ) is a finite Borel measure on (0, ∞).Then where and ρ (dλ) is a Borel measure satisfying

Theorem 4 The Green function of a fractional differential equation is completely monotonic on
Proof With β n = 1, the Green function of (2.3) has the Laplace transform (Podlubny [37], p. 158).Let Clearly, (6.18) satisfies (6.15) if β 0 > 0 and has Laplace transform Hence letting ζ = ip and R (0) = A −1 n in (6.13) we see that g (p) is the Laplace transform of a function which is completely monotonic on [0, ∞) for β 0 > 0. For β 0 = 0, the same approach is followed with the exception of using Theorem 8.5 of Okabe [34] instead of Theorem 5.10 of Inoue [26] .completes the proof.
It is not sufficient to know that we can write the Green function in the form (6.1).We need to know the measure µ (dλ).For this purpose, we rely on the following Theorem 1.3-5 of Djrbashian [18] for the two-parameter Mittag-Leffler function, which can be defined by the series expansion It is noted that Theorem 5 If ρ < 1 and µ ∈ (0, 1 + ρ) , then the following formula is true: We now consider some examples of fractional differential equations whose Green function has the representation (6.1).
Example 2 Following the description in Okabe [35], we consider a sphere of mass m and radius r moving in a fluid with viscosity η and density ρ.Denoting the velocity of the sphere by X (t), the random force by W (t) and the drag force by F (t) , Newton's second law reads Solving a linearised Navier-Stokes equation subject to incompressibility and sticky boundary conditions we have Newton's equation then becomes the Stokes-Boussinesq-Langevin equation where m * = m+ 2 3 πr 3 ρ is the effective mass of the sphere.If W (t) is white noise, that is if W (t) = dB (t)  dt (understood in the random distribution sense), where B (t) is Brownian motion, then the solution to the Stokes-Boussinesq-Langevin equation has the moving average representation where the kernel G (t) has the representation with

Almost sure convergence
We first prove the rate of convergence for the approximation to an Ornstein-Uhlenbecktype process.Lemma 3 Let q > 1 and let p be such that q −1 + p −1 > 1. Assume L (t) is of bounded p−variation, the difference between the Ornstein-Uhlenbeck-type process and its approximation defined by (6.9 -6.10) is uniformly bounded on [0, T ] by By elementary arguments the q−variation, q > 1, of f ∆ (s) is It follows that 1−1/q + λ∆ (6.31) for some finite constant C. Applying the Love-Young inequality to the difference of Y (n∆) and Y ∆ (n∆) we get for some finite constant C. When (n + 1) ∆ ≥ t > n∆ then from (6.10), Applying the (6.32) and the Love-Young inequality to (6.33) we get the difference Clearly this is bounded by (6.25) and proof is completed.Theorem 6 Suppose that K N is a sequence of compact sets growing to (0, ∞) , r N → 1 and ∆ N → 0 such that as N → ∞.Then, almost surely, Proof We first consider the error in approximating X (t) by where the inequality is derived using the Love-Young inequality.We now consider the error in approximating X K (t) by X π (t): For any λ 2 > λ 1 , where (6.43) is an application of the Love-Young inequality.Now,

Simulations
Sample paths of the fractional differential equation were generated in MATLAB using the algorithm described in Subsection 6.1.The simulations were carried out for three different types of Lévy motion: α−stable, inverse Gaussian and normal inverse Gaussian.
A common way to define an α−stable random variable is by its characteristic function:

ψ.Figure 1 :Figure 2 :Figure 3 :
Figure 1: A sample path of the two-term fractional differential equation driven by the generalized derivative of an α-stable process [13] Using the inequality 1 − e −r i N (r N −1)t ≤ r i N (r N − 1) t, (6.44) becomes|X K N (t) − X π (t)| ≤ 2C 1,p L (t) (p) (r N − 1)Using the triangle inequality we may combine equations (6.38), (6.46) and (6.49) to conclude that the almost sure convergence holds if the conditions of the theorem are satistfied.This completes the proof.Remark 4 In Example 1, if β n < 1, then K C N µ (dλ)0 and so a.s.convergence does not hold.However, mean-square convergence still holds by the results of Carmona et al.[13].