Monte carlo estimation of the solution of fractional partial differential equations

Remark 1.1

Given a Levy measure ν on ℝ₊ satisfying

the operator −_tD_a is defined by

When {Xsx}s⩾0 is Brownian motion, then A_x would be 12Δ, where Δ = ∑i=1d∂∂xi2. If {Ysa,t} is the deterministic drift, i.e. −tDa=−ddt and g = 0, then (1.1) becomes

the heat equation that we are more familiar with.

We assume {Xsx}s⩾0 is isotropic β-stable.(What ‘isotropic’ means is explained in Section 2, after Lemma 2.1) In this paper we shall investigate some properties of the representation (1.2) and its Monte-Carlo estimator, i.e.

where h > 0 is the step length, Ttk are iid samples of T_t, and tik = (i − 1)h. Note that we can sample the stopping time T_t (see Lemma 2.1 below), then sample the isotropic β-stable process {Xsx} and finally simulate the estimator (1.4).

In Section 2 we mainly focus on the situation when g = 0, i.e. the estimator now is

To make central limit theorem and Berry-Esseen bound hold, we only need to estimate the tail of the stable process at some stopping time, i.e. P[|XTtx|>s] for large s. And we begin with showing that the order of the tail of multidimentional stable distribution has the same order of the tail of each component of itself. In Section 3 we study the property of the Monte-Carlo estimator when the forcing term g ≠ 0. We estimate the upper bound of the second moment of the estimator and then, the L² error between the estimator and the solution. Besides, we use there properties to show that the central limit theorem holds using the triangular arrays. In Section 4 we give numerical examples, demonstrating the performance of our simulation algorithm.

Theorem 2.1

1. For all continuous functionϕ : ℝ^d → ℝ,

   uNt,x→a.s.ut,x,asN→∞.(2.2)

2. Let S_N(t, x) = N (u_N(t, x) − u(t, x))/σ(t, x) andWbe the standard normal distribution. Ifϕ(x) satisfiesϕx=O(|x|β2+δ),whereδ > 0, then the central limit theorem holds, i.e. for all bounded uniformly continuous funtionψ,

   EψSNt,x→EψWasN→∞.

3. LetYt,x:=ϕ(Ttx)−E[ϕ(XTtx)],denote 𝔼[Y(t, x)²] = σ(t, x)², 𝔼[∣Y(t, x)∣³] = ρ(t, x). Ifϕ(x) satisfiesϕx=O(|x|β3+δ),whereδ > 0, then for allC³functionsψ : ℝ → ℝ,

   |EψSNt,x−EψW|⩽0.433||ψ‴||∞ρt,xNσt,x3.

   HereC³means the space of functions with bounded third derivatives.

In other words, the central limit theorem can be written using convergence in distribution:

Since the estimator is unbiased, Theorem 2.1(i) holds because of the strong law of large numbers. For (ii), it is the standard central limit theorem and we only need to show that EϕXTtx2<∞. For (iii), it is a version of the Berry-Esseen bound and we need to show that E[|ϕ(XTtx)|3]<∞. These facts are evident if ϕ(x) is bounded. To deal with unbounded ϕ(x), let us recall the following fact: for any random variable U,

It is finite if ℙ[∣U∣ > t] = O(t^−(2+δ)), where δ is a positive constant. Now let us look back at our problems. Once we know the tail of XTtx and the growth rate of ϕ(x), the tail of ϕ(XTtx) would be clear as well as the finiteness of the moments of ϕ(XTtx).

Luckily, we have the following result:

Proposition 2.1

Assume that {X_s}_s⩾0is aβstable process, thenP[|XTtx|>u]=O(u−β).

To prove Proposition 2.1, we need a little lemma telling us that the distribution of T_t is analytically accessible:

Lemma 2.1

Denoteā := t − a, thenTt=da¯ηαwhereη ∼ S_α(1,1,0).

Proof

Note that Ysa,t=dt−s1/αη.{Tt>s}={Ysa,t>a}, since τ has monotone paths. Hence

□

Together with the facts that Xsx is β stable and Lemma 2.1, we have

Also we need Lemma 2.2 and Lemma 2.3 given below. Before that let us explain what ‘isotropic’ means in our assumption of {X_s}_s⩾0.

For d-dim β-stable random variable U = (U₍₁₎, …, U_(d)) on ℝ^d, there are a finte measure λ on sphere S and γ in ℝ^d such that the characteristic function of U satisfies

for β ≠ 1 and vice versa. Hence each component of U is 1-dim stable random variable and the stability index is still β. Besides, for 1-dim β-stable random variable V whose characteristic function has form

we use the notation V ∼ S_β(σ,ρ,μ). We say a d-dim stable random variable U is isotropic if its coordinates have the same distribution, i.e. U_(i) ∼ S_β(σ,ρ,μ) i = 1, …, d. We say a process {X_s}_s⩾0 is isotropic stable if X₁ is an isotropic stable random variable.

Lemma 2.2

LetU = (U₍₁₎, …, U_(d)) be an isotropicd-dimβ-stable random variable, andU_(i) ∼ S_β(σ,ρ,μ), then ℙ[∣U∣ > s] ∼ s^−βass → ∞.

Proof

Since {|U|=U12+⋯+Ud2>s}⊃{|U1|>s}, we have

Since {|U|>s}⊂{max1⩽i⩽d|Ui|>s/d}⊂∪i=1d{|Ui|>s/d}, we have

Now recall the well known result of the tail of 1-dim stable random variable: if V ∼ S_β(σ,ρ,μ), then

where Cβ=(∫0∞x−βsin⁡xdx)−1=1−βΓ2−βcosπβ/2 (see [15], Property 1.2.15).

Hence for any ϵ > 0, there exists some M, such that for all s > M and i = 1, …, d,

Hence for s>dM,

Therefore ℙ[∣X∣ > s] ∼ s^−β as s → ∞.□

Lemma 2.3

LetU, Vbe positive random variables such that

whereC₁ > C₂, then

Proof

Given a positive number M, there exsits T and ϵ > 0, such that for all t > T,

If we pick M big enough, we have ℙ[U − V > t] ⩾ Ct^−α for some constant C. On the other hand, for large t,

Therefore ℙ[U − V > t] ∼ t^−α as s → ∞.□

Lemma 2.2 tells us the order of tail of high dimensional stable process. Lemma 2.3 shows the order of the difference between certain random variables and we can apply it to the logarithm of (2.4), i.e. log∣X₁∣ + αβ logā − αβ logτ₁.

Now we can come back to the proof of Proposition 2.1.

Proof

Proof of Proposition 2.1. Now let us estimate the tail of XTtx. For large u > 0,

where A := log∣X₁∣, B := αβ log(η), r := log(u − ∣x∣) − αβ logā. (Note that for large u we have r > 0).

Let X₁ = (X₍₁₎, …, X_(d)) and X_(i) ∼ S_β(σ,ρ,μ), i = 1, …, d. By the Proof of Lemma 2.2, for any ϵ > 0, there exists some M, such that for all s > M and i = 1, …, d,

Hence for s>dM,

and for t>log⁡(dM),

Now let us discuss (2.7) in three conditions. For r>log⁡(dM):

1. When A > 0, B > 0, we have

   P[A−B>r,A>0,B>0]⩽P[A>r]=P[|X1|>es]⩽d1+β/2ϵ+Cβσβe−βr.(2.8)

2. When A > 0, B < 0, pick integer k = ⌊r/S⌋, and we divide the event {A + (−B) > r} into k parts:

   {A+−B>r}=⋃i=1k−1{A+−B>r,−B∈i−1kr,ikr}⋃{A+−B>r,−B>k−1kr}⊂⋃i=1k−1{A>k−ikr,−B∈i−1kr,ikr}⋃{−B>k−1kr}⊂⋃i=1k{A>k−ikr,−B>i−1kr}.(2.9)

   Hence,

   P[A+−B>r,A>0,B<0]⩽∑i=1kP[A>k−ikr,−B>i−1kr]=∑i=1kP[A>k−ikr]P[−B>i−1kr].(2.10)

   Recall that

   P[log⁡|X1|>k−ikr]⩽d1+β/2ϵ+Cβσβe−k−ikβr.(2.11)

   Using the result (3.7) that we shall mention later, we have

   E[η−2α]=(cos⁡(πα2))22Γ1+2α.(2.12)

   By the Markov inequality,

   P[αβlog⁡(η−1)⩾i−1kr]=P[η−1>eβαi−1kr]⩽E[τ1−2α](eβαi−1kr)2α⩽2e−2i−1kβr.(2.13)

   Combining (2.10), (2.11) and (2.13), we have

   P[A+−B>r,A>0,B<0]⩽∑i=1kd1+β/2ϵ+Cβσβe−k−ikβr2e−2i−1kβr=2d1+β/2ϵ+Cβσβ∑i=1ke−i−2kβre−βr⩽2d1+β/2ϵ+Cβσβeβr/k1−e−βr/ke−βr⩽2d1+β/2ϵ+Cβσβe2βSeβS−1e−βr.(2.14)

3. When A < 0, B < 0, then

   P[A−B>r,A<0,B<0]⩽P[A<0,−B>r]⩽P[−B>r]=P[αβlog⁡(η−1)>r]=P[η−1>eβαr]⩽E[η−α]eβαrα⩽e−βr.(2.15)

   Combining the three conditions above, we know that for large u,

   P[|XTtx|>u]⩽1+2e2βSeβS−1d1+β/2ϵ+Cβσβ+1e−βr=Oe−β(logu−|x|−αβlog⁡a¯)=Ou−β.(2.16)

   □

   Now let us finish the proof of our main result. First let us see this for Theorem 2.1(ii).

Proof

Proof of Theorem 2.1

and by (2.3)E[(ϕ(XTtx))2] is finite.□

For the proof of Theorem 2.1(iii), E[|ϕ(XTtx)|3] is finite because of the similar argument. For the rest of proof, see [14], page 356, Variant Berry-Esseen Theorem.

Remark 2.1

1. If C_α < C_βσ^β, by Lemma 2.3, we have

   P[A−B>r]⩾P[A−B>r,A>0,B>0]⩾Ct−β,(2.18)

   where C is a constant that can be chosen from the proof of Lemma 2.3. This result means the order t^−β is the best one.

2. In the proof of Proposition 2.1 we need r = log(u − ∣x∣) − αβ logā and r>log⁡(dM). Hence there exists some constant M₀ such that for u > M₀, (2.16) holds and M₀ has order d^1/2.

Besides, we can roughly give the upper bound of E[ϕ(XTtx)2].

Example 2.1

If ϕ(x) satisfies ϕx⩽|x|βδ+2, where δ > 0, then from Remark 2.1 we know that there exists some M₀ such that for all t > M₀,

where M1=1+e2βSeβS−1d1+β/2(ϵ+Cβσβ)+1,M2=2a¯−αβM1. Hence,

Note that M₀ has order d^1/2 and M⁽²⁾ has order d^1+β/2, This upper bound has order d^1+β/2.

Theorem 3.1

Assume|ϕ(x)|=O(|x|β2+δ)for ∣x∣ → ∞, where δ > 0.

1. E[uNht,x−ut,x2]→0asN → ∞, h → 0.

2. (CLT with a bias correction) LethN=N−2βγ,ut,x=EZt,xwhere

   Z(t,x)=ϕ(XTtx)+∫0TtgXsxds,

   andWbe the standard normal distribution, then for all bounded uniformly continuous functionψ,

   E[ψ(NuNhNt,x−ut,x/VarZt,x)]→E[ψW]asN→∞.

Let

be the approximation of Z(t, x). And let

where Yhk(t,x)=ϕXTtkx,k+∑i=1⌊Ttk/h⌋hgXtikx,k,k=1,…,N.Yhk(t,x) are the iid copies of Y_h(t, x). Note that for random variable U, let V be its approximation and V^k, k = 1, …, N be the iid copies of V. The L² error satisfies

Therefore, to estimate the L² error 𝔼[(u(t, x) − u_N(t, x))²], we only need to study varY_h(t, x) and 𝔼Z(t, x)-𝔼Y_h(t, x), and the following propositions answer these questions.

Proposition 3.1

There exists a constantMt,x1(depending ont, x) such thatVarY_h(t, x) ⩽ Mt,x1.

Proposition 3.2

There exists a constantMt,x2(depending ont, x) such thatE[|Zt,x−Yht,x|]⩽Mt,x2hγβ.

Proposition 3.3

There exists a constantMt,x3(depending ont, x) such thatE[|Zt,x−Yht,x|2]⩽Mt,x3h2γβ.

Sections 3.1 and 3.2 give proofs of these propositions. Section 3.3 is the proof of our CLT.

Remark 3.1

• Non-asymptotic confidence interval: Combining (3.1), Proposition 3.1 and Proposition 3.2 we have

  E[(ut,x−uNht,x)2]⩽1NMt,x1+(Mt,x2)2h2γβ,(3.2)

  where u(t, x) is the solution of problem (1.1). Now we can construct the confidence interval using the Markov inequality:

  P[|ut,x−uNht,x|>r]⩽E[(ut,x−uNht,x)2]/r2⩽1r21NMt,x1+(Mt,x2)2h2γβ.(3.3)

  Hence we can pick suitable N and h such that P[|ut,x−uNht,x|>r] < 1 − ϵ for some small ϵ.

• Asymptotic confidence interval: We can use CLT in Theorem 3.1 to get the asymptotic confidence interval. In other words, the central limit theorem can be written using convergence in distribution:

  NuNhNt,x−ut,x→dN0,VarZt,xas n→∞.(3.4)

  Once we have the upper bound M(t, x) of VarZt,x (e.g. see Example 2.1), it is easy to see that it yields a 100(1 − α)% asymptomic confidence interval uNhN±Mt,xNzα/2 for u(t, x), where z(t) satisfies Φ(z(t)) = 1 − t and Φ is the distribution function of the standard normal distribution. See Section 4.3 for a simple example.

Before the calculation, we need the following results (see [6], page 162):

1. For constants c > 0, η ∈ (−1,β) and a symmetric β-stable 1-dim process U_t with 𝔼[e^izU_t}] = e^{−tc∣z∣β}, we have

   E[|Ut|η]=tcη/β2ηΓ1+η2Γ1−ηβπΓ1−η2.(3.5)

   Recall that each component of X₁ − X₀, denoted by X_(j), is symmetric with 𝔼[e^izX_(j)] = e^{−c∣z∣β} and c > 0.

2. If 0 < α < 1 and {X_t} is a stable subordinator with 𝔼[e^−uX_t] = e^−tc′uα, where c′ is some constant, then for −∞ < η < α,

   E[Xtη]=tc′η/αΓ(1−ηα)Γ1−η.(3.6)

   Since η ∼ S_α(1,1,0), we have E[e−uη]=exp⁡{−1cosπα2uα} (see [15], Proposition 1.2.12). Hence,

   E[ηη]=1cosπα2η/αΓ(1−ηα)Γ1−η.(3.7)