Iterated stochastic processes : simulation and relationship with high order partial differential equations

TL;DR

This paper introduces an Euler-based simulation algorithm for iterated stochastic processes, extends the Feynman-Kac formula to higher-order PDEs, and demonstrates their connection to complex stochastic models.

Contribution

It presents a natural, convergent simulation method for iterated processes and generalizes the Feynman-Kac formula to relate these processes to high-order PDEs.

Findings

The simulation algorithm converges almost surely with rate 1/4.

Extended Feynman-Kac formula to higher-order PDEs.

Simulated solutions to fourth order PDEs using iterated processes.

Abstract

In this paper, we consider the composition of two independent processes : one process corresponds to position and the other one to time. Such processes will be called iterated processes. We first propose an algorithm based on the Euler scheme to simulate the trajectories of the corresponding iterated processes on a fixed time interval. This algorithm is natural and can be implemented easily. We show that it converges almost surely, uniformly in time, with a rate of convergence of order 1/4 and propose an estimation of the error. We then extend the well known Feynman- Kac formula which gives a probabilistic representation of partial differential equations (PDEs), to its higher order version using iterated processes. In particular we consider general position processes which are not necessarily Markovian or are indexed by the real line but real valued. We also weaken some assumptions from previous works. We show that intertwining diffusions are related to transformations of high order PDEs. Combining our numerical scheme with the Feynman-Kac formula, we simulate functionals of the trajectories and solutions to fourth order PDEs that are naturally associated to a general class of iterated processes.

Equations

$$Z_t:=X(|Y_t|),t\in I.$$

$$X_{2}(t):=\left\{\begin{array}[]{lcl}X(t)&\text{if}&t\geq 0\\ X^{\prime}(-t)&\text{if}&t<0\\ \end{array}\right.$$

$$\forall 0<\alpha <1/4,\underset{n\to +\infty }{\lim }{n}^{\alpha }\underset{t\in [0,T]}{\sup }|{\tilde{Z}}^n(t)-Z(t)|=0,a.s.$$

$$\underset{n}{\mathrm{limsup}}{n}^{\alpha }|Y_T-Y_T^n|=\infty ,a.s.$$

$$d{X}_t=b_X(t,X_t)dt+{\sigma }_X(t,X_t)d{W}_t,t\ge 0,$$

$$d{Y}_t=b_Y(t,Y_t)dt+{\sigma }_Y(t,Y_t)d{W}_t^{\prime },t\in [0,T],$$

$$\forall p\ge 1,\mathbb{E}|X_0{|}^p+\mathbb{E}|Y_0{|}^p<\infty .$$

$$\forall t\ge 0,$$

$$\forall s,t\ge 0,$$

$$\forall t\in ]t_k,t_{k+1}],Y_t^n:=Y_{t_k}^n+{\sigma }_Y(t_k,Y_{t_k}^n)(W_t^{\prime }-W_{t_k}^{\prime })+b_Y(t_k,Y_{t_k}^n)(t-t_k).$$

$$\mathbb{E}(\underset{t\in [0,T]}{\sup }|Z_t-Z_t^n{|}^{2p})\le C({\Delta }_T{)}^{\rho },$$

$$\underset{n\to +\infty }{\lim }{n}^{\alpha }su{p}_{t\in [0,T]}(|Z_t^n-Z_t|+|{\tilde{Z}}_t^n-Z_t|)=0,a.s.$$

$$\mathbb{E}(\underset{t\in [0,M]}{\sup }|X_t-X_t^n{|}^{2p})\le C(1+M^p)M({\Delta }_M{)}^{\gamma }{e}^{CM}$$

$$\begin{gathered} \mathbb{E}\left(\underset{t\in [0,T]}{\sup }|Z_t^n-Z_t{|}^{2p}\right)\le 2^{2p-1}\mathbb{E}\left(\underset{t\in [0,T]}{\sup }|X(|Y_t^n|)-X(|Y_t|){|}^{2p}\right) \\ +2^{2p-1}\mathbb{E}\left(\underset{t\in [0,T]}{\sup }|X^n(|Y_t^n|)-X(|Y_t^n|){|}^{2p}\right). \end{gathered}$$

$$\Gamma ={\int }_{[0,M{]}^2}\frac{|X_a-X_b|}{|a-b{|}^{m+2}}dadb,$$

$$\mathbb{E}(\Gamma )\le H{c}_1\frac{1}{(p-1-m)}(\frac{M}{2}{)}^{p-m},$$

$$\mathbb{E}\left(\underset{t\in [0,T]}{\sup }|X(|Y_t^n|)-X(|Y_t|){|}^{2p}\right)$$

$$\mathbb{E}[{\delta }^{2m}]\le ({\Delta }_T{)}^{m\inf (2\beta ,1)}.$$

$$\mathbb{E}\left(\underset{t\in [0,T]}{\sup }|X(|Y_t^n|)-X(|Y_t|){|}^{2p}\right)\le \tilde{C}{\left(\mathbb{E}\left[P(M)e^{CM}\right]\right)}^{1/2}({\Delta }_t{)}^{(p-1)\inf \{2\beta ,1\}/2}.$$

$$\mathbb{E}\left(\underset{t\in [0,T]}{\sup }|X^n(|Y_t^n|)-X(|Y_t^n|){|}^{2p}\right)$$

$$\le n^{-\gamma }\mathbb{E}\left[C(1+M^p)M^{1+\gamma }{e}^{CM}\right].$$

$$\mathbb{E}{e}^{\lambda {\sup }_{t\in [0,T]}\max \{|Y_t|,|Y_t^n|\}}\le \mathbb{E}{e}^{\lambda {\sup }_{t\in [0,T]}|Y_t|}+\mathbb{E}{e}^{\lambda {\sup }_{t\in [0,T]}|Y_t^n|},$$

$$\mathbb{E}{e}^{\lambda {\sup }_{t\in [0,T]}|Y_t|}$$

$$\le e^{\lambda }+\mathbb{E}\left(\sum _{k\ge 0}\frac{{\lambda }^k}{k!}(\underset{t\in [0,T]}{\sup }|Y_t|{)}^{2k}\right).$$

$$e^{\lambda }+\mathbb{E}\left(\sum _{k\ge 0}\frac{{\lambda }^k}{k!}(C(1+T^k)\left\{\mathbb{E}(|Y_0{|}^{2k})+\left(1+\mathbb{E}(|Y_0{|}^{2k})\right)T^k{e}^{CT}\right\}\right),$$

$$\forall \alpha <\ell ,\underset{n\to +\infty }{\lim }{n}^{\alpha }\underset{t\in [0,T]}{\sup }|f_n(t)-f(t)|=0$$

$$\forall \alpha <\ell ,\underset{n\to +\infty }{\lim }{n}^{\alpha }\underset{t\in [0,T]}{\sup }|{\bar{f}}_n(t)-f(t)|=0.$$

$$\forall 0\le \alpha <\frac{1}{4},\underset{n\to +\infty }{\lim }{n}^{\alpha }su{p}_{t\in [0,T]}|{\bar{X}}^n(|{\bar{Y}}_t^n|)-X(|Y_t|)|=0,a.s.,$$

$$|X_{|Y_t|}-{\bar{X}}_{|{\bar{Y}}_t^n|}^n|\le C|Y_t-{\bar{Y}}_t^n{|}^{\beta }+|X_{|{\bar{Y}}_t^n|}-{\bar{X}}_{|{\bar{Y}}_t^n|}^n|.$$

$$n^{\alpha }\underset{[0,T]}{\sup }|X_{|Y_t|}-{\bar{X}}_{|{\bar{Y}}_t^n|}^n|\le C(n^{\alpha /\beta }\underset{[0,T]}{\sup }|Y_t-{\bar{Y}}_t^n|{)}^{\beta }+n^{\alpha }\underset{[0,M^{\prime }+1]}{\sup }|X_t-{\bar{X}}_t^n|).$$

$$\forall 0\le \alpha <\frac{1}{2},\underset{n\to +\infty }{\lim }{n}^{\alpha }su{p}_{t\in [0,T]}|{\tilde{X}}_t^n-X_t|=0,a.s.$$

Full text

Iterated stochastic processes : simulation and relationship with high order partial differential equations.

Michèle Thieullen Laboratoire de Probabilités et Modèles Aléatoires, Université Pierre et Marie Curie, 4, Place Jussieu, 75252 Paris cedex 05, Francemichele.thieullen@upmc.fr   

Alexis Vigot $\text{}^{\star}$ alexis.vigot@courriel.upmc.fr

In this paper, we consider the composition of two independent processes: one process corresponds to position and the other one to time. Such processes will be called iterated processes. We first propose an algorithm based on the Euler scheme to simulate the trajectories of the corresponding iterated processes on a fixed time interval. This algorithm is natural and can be implemented easily. We show that it converges almost surely, uniformly in time, with a rate of convergence of order $1/4$ and propose an estimation of the error. We then extend the well known Feynman-Kac formula which gives a probabilistic representation of partial differential equations (PDEs), to its higher order version using iterated processes. In particular we consider general position processes which are not necessarily Markovian or are indexed by the real line but real valued. We also weaken some assumptions from previous works. We show that intertwining diffusions are related to transformations of high order PDEs. Combining our numerical scheme with the Feynman-Kac formula, we simulate functionals of the trajectories and solutions to fourth order PDEs that are naturally associated to a general class of iterated processes.

Keywords : Iterated process, Euler scheme, high order partial differential equation, Feynman-Kac formula, diffusion processes, iterated Brownian motion.

1 Introduction

In the present paper we address iterated processes for which we propose a numerical scheme to simulate their trajectories and investigate their relationship with high order partial differential equations (PDEs). Then we use our scheme to approximate numerically the solution of such PDEs. Consider two independent processes $(X_{t})_{t\geq 0}$ and $(Y_{t})_{t\in I}$, where $I$ is an interval, defined on a probability space $(\Omega,\mathcal{F},\mathbb{P})$ and taking values respectively in $\mathbb{R}^{d}$ and $\mathbb{R}^{n}$. We define an interated process $Z$ by replacing the time index of $X$ by $|Y|$ where $|\cdot|$ denotes a norm on $\mathbb{R}^{n}$ as follows

[TABLE]

In the sequel we call $(X_{t})$ (resp. $(Y_{t})$) the position (resp. time) process. Different types of iterated processes have been considered recently for instance by Allouba [1, 2], Burdzy [8, 9, 10], Khoshnevisan and Lewis [17]. Some of these authors extend $X$ to a process indexed by the real line as follows

[TABLE]

where $X^{\prime}$ is an independent copy of $X$. The corresponding iterated process is then defined by $X_{2}(Y_{t})$.

In order to study properties of iterated processes, in particular to be able to simulate solutions of related PDEs, we would like to simulate their trajectories and to control the error. We haven’t found any reference on the simulation of such processes in the literature. In the first part of the paper we propose a scheme based on the classical Euler scheme to simulate the trajectories of the iterated process $(Z_{t})$ when $X$ and $Y$ are two diffusions. For every $T>0$ and positive integer $n$, we construct a process $\tilde{Z}^{n}$ piecewise linear on a regular partition of $[0,T]$ with mesh $T/n$. We will see that this scheme converges a.s. uniformly as $n$ tends to infinity on $[0,T]$ to $(Z_{t})$ with a rate of convergence of order $1/4$, which means that

[TABLE]

It seems that the order of convergence in (3) for iterated processes satisfying our assumptions cannot be better than $1/4$. Indeed for a Brownian motion $Y$, we have for all $\alpha>1/2$ (see [12])

[TABLE]

The convergence of the scheme being uniform, we can use it to simulate quantities depending on the whole trajectory of the process, not only on its value at a fixed time, like for instance the variations, the mean of a functional of the process or the measure of a subset of the function space $[0,T]^{\mathbb{R}}$ with respect to the law of the process. Our scheme is easy to implement since it only requires the simulation of $2n$ Gaussian random variables and so it is adapted to methods using simulations of many trajectories such as Monte Carlo methods. It can be extended to processes given by construction (2). We have chosen the Euler scheme because of its simplicity. The Milshtein approximation for instance, doesn’t provide a better rate of convergence and is more complex.

In the subsequent part of the paper we address the connection between iterated processes and high order PDEs. There are in the literature several papers which prove that the density of specific iterated processes satisfies a Fokker-Planck type PDE ([5], [11]), establish a Feynman-Kac type formula ([2], [14], [18]) for such processes or associate them to fractional equations ([4], [19], [20]). We extend these works to general position processes in particular not necessarily Markovian (cf. section 3.3).

When the time process is a Brownian motion or an $\alpha$-stable process, the equation we obtain reduces to the equation obtained respectively by [2] and [18]. However our result holds under more general assumptions on the initial value. Moreover in this particular case we prove that the solution of the PDE is unique. It is difficult to extend the results of this paper in a general setting when the density of the time process satisfies a PDE whose coefficients depend on the spatial variable. We provide examples of such time processes.

We also consider transformations going from one high order PDE to another one and we point out the relationship with intertwining diffusions. However, these results (as well as those in [2] and [18]) concern PDEs which contain the initial value and some of its derivatives. The PDE obtained in [14] does not have this drawback but the underlying iterated process takes values in the complex plane. Extending the construction (2), we obtain a Feynman-Kac formula where the initial value does not appear any longer in the PDE. Moreover the underlying process is real valued. An application of our result provides a stochastic representation of a solution of the Euler-Bernouilli beam equation when one considers iteration of a Brownian motion by an independent Cauchy process.

In the last part of the paper, we implement the algorithm for the iterated Brownian motion (IBM) given by (1) when $X$ and $Y$ are two independent Brownian motions. We simulate its third and fourth order variations and the solution of the corresponding fourth order PDE.

The present paper is organized as follows. In section 2 we describe our scheme, we prove its uniform convergence and study the error. Section 3 is devoted to the relationship between iterated processes and PDEs in the spirit of [2] and [18]. Transformations of PDEs are addressed in section 4. In section 5, inspired by [14], we work with construction (2) in order to obtain PDEs that do not contain the initial value. Section 6 implements our scheme in the IBM case. The Appendix (section 7) contains some auxiliary proofs.

2 Numerical scheme.

In this section we describe our numerical scheme and study its convergence. We consider two independent time inhomogeneous diffusion processes $X$ and $Y$. Our scheme is based on the idea that a pathwise approximation of $X(|Y|)$ may be obtained by the composition of the respective Euler approximations of $X$ and $Y$. More precisely for $T>0$ and a positive integer $n$, the Euler scheme for $Y$ on $[0,T]$ yields a continuous approximation $Y^{n}$. If $M_{n}$ denotes the maximum of $|Y^{n}|$ on $[0,T]$, the Euler scheme for $X$ on $[0,M_{n}]$ provides an approximation $X^{n}$. We prove below that a piecewise linear approximation $\tilde{Z}^{n}$ of the composition $Z^{n}:=X^{n}(|Y^{n}|)$ converges a.s. uniformly on $[0,T]$ as $n$ tends to infinity to $Z=X(|Y|)$ with rate of order $\frac{1}{4}$.

We start with the stochastic differential equations (SDEs) satisfied by the position and time processes. Let us assume that the position process $(X_{t})_{t\geq 0}$ starting from $X_{0}$ satisfies

[TABLE]

where $b_{X}:\mathbb{R}\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d}$, $\sigma_{X}:\mathbb{R}\times\mathbb{R}^{d}\rightarrow\mathbb{R}^{d\times p}$ and $(W_{t})$ is a standard $\mathbb{R}^{p}$ valued Brownian motion. The time process $(Y_{t})_{t\in[0,T]}$ starting from $Y_{0}$ solves

[TABLE]

with $b_{Y}:\mathbb{R}\times\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$, $\sigma_{Y}:\mathbb{R}\times\mathbb{R}^{n}\rightarrow\mathbb{R}^{n\times q}$ and $(W^{\prime}_{t})$ a standard $\mathbb{R}^{q}$ valued Brownian motion independent from $W$. Several assumptions will be used throughout the paper. Regarding the starting points, we assume that

[TABLE]

The coefficients of (4) and (5) are supposed to enjoy Lipschitz continuity in space as well as at most linear growth and Hölder continuity in time. For simplicity, we write down these assumptions for $b_{X}$ and $\sigma_{X}$ only. For $b_{Y}$ and $\sigma_{Y}$, the only difference is that the time variable in (5) belongs to the compact interval $[0,T]$. We assume that there exist two positive real numbers $K$ and $\beta$ such that

[TABLE]

[TABLE]

Under these assumptions equations (4) and (5) have unique strong solutions (cf. [16]).

Let us recall that the Euler scheme for $Y$ on $[0,T]$ given the regular subdivision $t_{0}=0<t_{1}<\ldots<t_{n}=T$ with mesh $\Delta_{T}=T/n$ is defined recursively by $Y_{0}^{n}:=Y_{0}$ and

[TABLE]

Let us define $M_{n}:=\sup_{t\in[0,T]}|Y_{t}^{n}|$, which is a.s. finite (cf. Proposition 8) and denote by $X^{n}$ the approximation resulting from the Euler scheme for $X$ on $[0,M_{n}]$. We propose to approximate $X(|Y_{t}|)$ by $X^{n}(|Y^{n}_{t}|)$. The following statement provides an estimate of the mean error associated to this approximation, uniformly on $t\in[0,T]$. In the following, $n$ is assumed to be such that $\Delta_{T}=T/n\leq 1$.

Theorem 1

Let us fix $T>0$ and assume that $X$ and $Y$ satisfy (4)-(9). Let $n\geq 1$ be an integer. For $t\in[0,T]$ we set $Z_{t}:=X(|Y_{t}|)$ and $Z_{t}^{n}:=X^{n}(|Y_{t}^{n}|)$. If $\sup_{t\in[0,T]}\max\{|Y_{t}|,|Y_{t}^{n}|\}$ has finite Laplace transform,

[TABLE]

where $C$ is a constant depending only on $K$, $p$, $T$ and $\mathbb{E}|X_{0}|^{2p}$ and $\rho$ is defined by $\rho=(p-1)\inf\{1,2\beta\}/2$. In particular (11) holds true if $Y_{0}^{2}$ has a finite Laplace transform.

In Theorem 1, we have used the exact composition of the respective Euler schemes of $X$ and $Y$ to approximate $X(|Y|)$ by $X^{n}(|Y^{n}|)$. Actually another approximation easier to implement may be chosen. Indeed, let us consider $\bar{Y}_{t}^{n}$, a step function defined on $[0,T]$ by $\bar{Y}_{t}^{n}=Y^{n}(kT/n)$ for $t\in[kT/n,(k+1)T/n]$ and $\bar{X}^{n}$, a step function defined on $[0,M_{n}]$ by $\bar{X}_{t}^{n}=X^{n}(kM_{n}/n)$ for $t\in[kM_{n}/n,(k+1)M_{n}/n]$. Let $\tilde{Z}^{n}$ be the linear interpolation between the points $\bar{X}^{n}(\bar{Y}^{n}(kT/n))$ for $k=0,1,\ldots,n$. The computation of $\tilde{Z}^{n}$ is facilitated by the use of piecewise constant processes whose composition is simpler than the composition of piecewise linear ones. We now show that $Z^{n}$ and $\tilde{Z}^{n}$ both converge uniformly to $Z$ with order $\frac{1}{4}$.

Theorem 2

Let $X$ and $Y$ satisfy (4)-(9). Then $\forall\,\,0\leq\alpha<\frac{1}{4},$

[TABLE]

The proof of Theorem 1 relies on the following proposition.

Proposition 1

*(cf. [12])

Under assumptions (4)-(9),*

[TABLE]

where $\Delta_{M}:=M/n$, $\gamma:=p\sup\{1,2\beta\}$ if $\Delta_{M}\geq 1$ and $\gamma:=p\inf\{1,2\beta\}$ otherwise. $C$ is a constant depending only on $K$, $p$, $n$ and $\mathbb{E}|X_{0}|^{2p}$.

In the sequel we need the explicit dependance of (13) on the parameter $M$, since the interval on which we approximate $X$ is random. For the sake of completeness, we provide the proof of Proposition 1 in the Appendix (section 5).

Proof of Theorem 1.

In the proof $C$ and $\tilde{C}$ denote constants depending only on $K$, $p$, $n$, $T$ and $\mathbb{E}|X_{0}|^{2p}$ which may vary from line to line. Remember that $Z_{t}=X(|Y_{t}|)$ and $Z^{n}_{t}=X^{n}(|Y_{t}^{n}|)$. We start with

[TABLE]

Let us bound the first term on the RHS. Define $M:=\sup_{t\in[0,T]}\max\{|Y_{t}|,|Y_{t}^{n}|\}$ and $\delta:=\sup_{t\in[0,T]}|Y_{t}-Y_{t}^{n}|$. Let us recall that if $M$ were fixed, Garsia-Rodemich-Rumsey Lemma (cf. [22]) would provide a random variable $\Gamma$ such that $|X_{s}-X_{t}|^{2p}\leq c\,\Gamma\,|t-s|^{m},$ for all $s,t$ in $[0,M]$ whatever $m\in]0,p-1[$. In this inequality $c$ is a constant depending on $p$ and $m$. Moreover this lemma provides the following expression for $\Gamma$,

[TABLE]

as well as the estimate

[TABLE]

where $c_{1}$ is a universal constant and $H$ denotes the right-hand side of (58). Here the pair $(M,\delta)$ is random. However by definition it is independent of $X$. Therefore we can condition on $(M,\delta)$ and write

[TABLE]

In the latter expression (obtained by Hölder inequality) we have set $P(x):=(1+x^{p})^{2}(\frac{x}{2})^{2(p-m)}$. Remember now that $\delta$ is equal to $\sup_{t\in[0,T]}|Y_{t}-Y_{t}^{n}|$. Proposition 14 in [12] implies

[TABLE]

Therefore we obtain the following upper bound for the first term:

[TABLE]

The second term can be easily dominated as follows. Indeed notice first that $\sup_{t\in[0,T]}|X^{n}(|Y_{t}^{n}|)-X(|Y_{t}^{n}|)|\leq\sup_{t\in[0,M]}|X_{t}^{n}-X_{t}|$. Conditioning by $M$ (which is independent of $X$) and applying Proposition 1 we obtain,

[TABLE]

where $\gamma=p\inf\{1,2\beta\}$. Combining (14) with inequalities (16) and (17) gives the result.

Let us now prove that the assumption $\sup_{t\in[0,T]}\max\{|Y_{t}|,|Y_{t}^{n}|\}$ has finite Laplace transform is satisfied if $Y_{0}^{2}$ has a finite Laplace transform.

Take $\lambda>0$. Then

[TABLE]

[TABLE]

Using Proposition 8 (see the Appendix), this latter term is dominated by

[TABLE]

which is finite since $Y_{0}^{2}$ has a finite Laplace transform. The same reasoning can be applied to $Y_{t}^{n}$ which satisfies (10).

In the proof of Theorem 2 we use the following lemma which enables us to approximate a given function by a sequence of step functions, instead of an arbitrary sequence, without changing the rate of convergence. This lemma is proved in the Appendix.

Lemma 1

Let $f$ and $(f_{n})_{n}$ denote functions defined on $[0,T]$. For all integer $n\geq 1$, define $\bar{f}_{n}(t):=f_{n}(kT/n),\,\,\forall t\in[kT/n,(k+1)T/n]$. Let $\ell>0$, and suppose that $f$ is Hölder continuous with exponent $\beta$ for every $\beta<\ell$. Then,

[TABLE]

is equivalent to

[TABLE]

Proof of Theorem 2. We prove the statement in detail for ${\tilde{Z}}_{n}$ (the proof for $Z_{n}$ is similar). This is equivalent to prove that

[TABLE]

thanks to Lemma 1 applied to $Z_{t}=X(|Y_{t}|)$ which is locally Hölder continuous with exponent $\beta$ for all $0<\beta<1/4$.

Let $\alpha\in]0,1/4[$. Let us take $\epsilon>0$ and $M^{\prime}>0$ such that $\mathbb{P}(\sup_{t\in[0,T]}|Y_{t}|\geq M^{\prime})\leq\epsilon$ and set $A:=\{\sup_{t\in[0,T]}|Y_{t}|<M^{\prime}\}$. Since $\bar{Y}^{n}$ converges a.s. uniformly to $Y$ on $[0,T]$, there exists $N$ such that a.s., for all $n\geq N$, $\sup_{t\in[0,T]}|\bar{Y}_{t}^{n}-Y_{t}|<1$. In particular the inequality $\sup_{t\in[0,T]}|\bar{Y}_{t}^{n}|<1+M^{\prime}$ holds a.s. on $A$. Since $X$ is Hölder continuous on $[0,M^{\prime}+1]$ of exponent $\beta$ for all $\beta\in]2\alpha,\frac{1}{2}[$, there exists a random value $C$ such that for all $n\geq N$, a.s. on $A$,

[TABLE]

Taking the supremum and multiplying by $n^{\alpha}$, we obtain that for all $n\geq N$, a.s. on $A$,

[TABLE]

It is proved in [12] that

[TABLE]

Lemma 1 applied to $X$ implies that the same convergence holds for $(\bar{X}^{n})$ instead of $\tilde{X}^{n}$ namely

[TABLE]

The same argument applies to $Y$. This leads to

[TABLE]

Therefore

[TABLE]

We conclude the proof by letting $\epsilon$ tend to [math].

3 PDEs for iterated processes with general time and position processes

In this section, we connect iterated processes and high order PDEs in the spirit of [2], [18]. However we consider a general framework where the position or time process is not necessarily Markovian. It seems difficult to prove a general statement when the density $p_{Y}$ of the time process satisfies a PDE with space dependent coefficients. In this case we treat an example. In particular we address the iteration by a Brownian motion with drift and by an Ornstein-Uhlenbeck process as time process. When the time process is a Brownian motion or an $\alpha$-stable process, the equation we obtain reduces to the equation obtained respectively by [2] and [18]. However our result holds under more general assumptions on the initial value. Moreover in this particular case we prove that the solution of the PDE is unique.

3.1 Iteration by a general time process

In this section the time process $(Y_{t})_{t\geq 0}$ is real valued, starts from [math] at time [math] and satisfies the following assumptions.

(A1) $Y_{t}$ admits a density denoted below by $p_{Y}(t,0,u)$ (or simply $p_{Y}$) for all $t>0$,

(A2) for all $t>0$, $p_{Y}(t,0,\cdot)\in\mathcal{C}^{r}(\mathbb{R})$ and $\partial^{i}_{t}\partial^{j}_{u}\,p_{Y}(t,0,\cdot)\in L^{1}(\mathbb{R})$ for all $0\leq i\leq q,\ 0\leq j\leq r$, for some integers $q,r$,

(A3) the even part $\mathcal{E}(p_{Y})$ of $p_{Y}$ satisfies

[TABLE]

for some integer $p$ and real valued functions of the time variable $\alpha_{k}$, $\beta_{i,j}$ (integers $q,r$ are those appearing in (A2)).

Examples satisfying (A1)-(A3) are provided by

• —

Brownian motion corresponding to $p=1,\,q=0,\,r=2,\,$ $\alpha_{1}=1,\,\beta_{0,1}=0,\,\beta_{0,2}=1/2$,

• —

the Cauchy process p=2,\,q=0,\,r=2,\,$$\alpha_{1}=0,\,\alpha_{2}=1,\, $\beta_{0,1}=0,\,\beta_{0,2}=-1$,

• —

the telegraph process whose density satisfies

[TABLE]

• —

the sum of a telegraph process and an independent Brownian motion. In this case $p_{Y}$ satisfies (see [6])

[TABLE]

Actually in the first three cases the density $p_{Y}(t,0,\cdot)$ itself is even for all $t$.

If $p_{Y}$ satisfies (18) and if moreover the orders of the partial derivatives w.r.t. the space variable $u$ on the right hand side of (18) are even, then $\mathcal{E}(p_{Y})$ satisfies (18) as well.

Another example is obtained when $Y$ is a Brownian motion with variance $\sigma^{2}$ and constant drift $\mu$. From the Fokker-Planck equation satisfied by $p_{Y},$

[TABLE]

we deduce that its even and odd parts satisfy the system

[TABLE]

which implies the following particular case of (18)

[TABLE]

In this subsection the position process is a Markov process $(X_{t}^{x})_{t\geq 0}$ starting from a given $x\in\mathbb{R}$ with semigroup $P^{t}$ and infinitesimal generator $\mathcal{L}$ with domain $D(\mathcal{L})$. We denote by $B(\mathbb{R}^{d},\mathbb{R}^{d})$ the set of bounded functions from $\mathbb{R}^{d}$ to $\mathbb{R}^{d}$ and define $B_{0}:=\{f\in B(\mathbb{R}^{d},\mathbb{R}^{d}),\|P^{t}f-f\|_{\infty}\rightarrow 0\text{\ as\ }t\downarrow 0\}$. $B_{0}$ is a closed subset of $B(\mathbb{R}^{d},\mathbb{R}^{d})$ containing $D(\mathcal{L})$. In the following, we will use iterations of $\mathcal{L}$:

[TABLE]

We state below the results of this section (Theorem 3 and its corollaries). Their proofs are provided in section 3.2.

Theorem 3

Let $Y$ be a real valued process independent of $X^{x}$ whose density $p_{Y}$ satisfies assumptions (A1)-(A3). Let $f\in D(\mathcal{L}^{r-1})$ where $r$ is the highest order of the partial derivatives in $u$ appearing in (18). Define $v(t,x):=\mathbb{E}\left[f(X^{x}(|Y_{t}|))\right]$. Then $v\in D(\mathcal{L}^{r})$ and satisfies

[TABLE]

on $]0,+\infty[\times\mathbb{R}^{d}$ where $\mathcal{L}$ acts on $x$ and ${\cal B}_{f}(t,x)$ is the boundary term

[TABLE]

Restricting to $Y$ a Brownian motion, we recover as a corollary the following result of [2]. Moreover we show that in this case (19) admits a unique solution under weaker assumptions on the initial condition.

Corollary 2

Let $f\in D(\mathcal{L})$ and $B$ a Brownian motion independent of $X^{x}$. Then, $v(t,x):=\mathbb{E}\left[f(X^{x}(|B_{t}|))\right]$ is an element of $D(\mathcal{L}^{2})$ and is the unique solution in $D(\mathcal{L}^{2})$ of

[TABLE]

Corollary 3

Let $f\in D(\mathcal{L}^{3})$ and $\,Y$ be a Brownian motion with drift $\mu$ and diffusion coefficient $\sigma$, independent of $X^{x}$. Then, there exist two functions of time, $\alpha$ and $\beta$, such that $v(t,x):=\mathbb{E}\left[f(X^{x}(|Y_{t}|))\right]$ satisfies

[TABLE]

with $v(t,\cdot)\in D(\mathcal{L}^{4})$, $\forall\,t>0$. By definition $v(0,x)=f(x)$.

The time processes considered in the above statements are associated to PDEs whose coefficients may depend on the time variable but not on the spatial variable. It seems difficult to extend Theorem 3 to a large class of time processes. For instance, if $Y$ is an Ornstein-Uhlenbeck process issued from [math] satisfying $dY_{t}=-\frac{1}{2}Y_{t}dt+dW_{t}$, its density satisfies the Fokker-Planck equation

[TABLE]

Since some coefficients of this equation are functions of the spatial variable $u$, we cannot apply Theorem 3 directly. However in this particular case a rewriting of the equation leads to the following result. The difficulty is that in general such a rewriting is not possible.

Proposition 4

Let $X^{x}$ be a $\mathbb{R}^{d}$-valued Markov process with infinitesimal generator $\mathcal{L}$. Let $Y$ an Ornstein-Uhlenbeck process independent of $X$, satisfying $dY_{t}=-\frac{1}{2}Y_{t}dt+dW_{t}$, $Y_{0}=0$. Then, for $f\in D(\mathcal{L})$, the function $v(t,x):=\mathbb{E}\left[f(X^{x}(|Y_{t}|))\right]$ satisfies $v(t,\cdot)\in D(\mathcal{L}^{2})$ for all $t>0$ and solves

[TABLE]

with initial condition $v(0,x)=f(x)$.

Corollary 5

Let $Y$ be a telegraph process with parameters $\lambda>0$ and $v>0$, independent of $X$ and $f\in D(\mathcal{L})$. Then, $u(t,x)=\mathbb{E}\left[f(X^{x}(|Y_{t}|))\right]$ is solution of

[TABLE]

3.2 Proofs

We now come to the proofs of Theorem 3, Corollary 2 and Proposition 4. We start with the proof of Theorem 3 which relies on the following lemma proved in the Appendix (section 7).

Lemma 2

Let $g\in L^{1}(\mathbb{R})$ such that $g$ is differentiable with $g^{\prime}\in L^{1}(\mathbb{R})$. Then for all $f\in B_{0}$, the function $F(g):x\mapsto\int_{0}^{\infty}g(s)P^{s}f(x)ds$ belongs to $D(\mathcal{L})$ and

[TABLE]

More generally, if $f\in D(\mathcal{L}^{n-1})$ and for all $k\in\{0,\ldots,n\}$, $g^{(k)}\in L^{1}(\mathbb{R})$, then $F(g)\in D(\mathcal{L}^{n})$ and

[TABLE]

Proof of Theorem 3. The independence of $X$ and $Y$ implies that

[TABLE]

and the fact that $\mathcal{E}(p_{Y})(t,0,u)=\frac{1}{2}(p_{Y}(t,0,u)+p_{Y}(t,0,-u))$ yield

[TABLE]

For fixed $t>0$ define $g(\cdot):=\mathcal{E}(p_{Y})(t,0,\cdot)$. Then $v(t,\cdot)$ coincides with $2F(g)$. From (A3) we obtain

[TABLE]

where we make explicit the fact that $\alpha_{k}$ and $\beta_{i,j}$ depend only on time. Applying (24) to each $0\leq j\leq r$, we obtain

[TABLE]

We can now conclude using the fact the even function $g\equiv\mathcal{E}(p_{Y})(t,0,\cdot)$ satisfies

[TABLE]

for all $0\leq\ell\leq j-1$.

Proof of Corollary 2 We deduce from Theorem 3 that $v\in D(\mathcal{L}^{2})$ and it satisfies (21). Let $\tilde{v}$ be another solution of (21) belonging to $D(\mathcal{L}^{2})$. Then $u=v-\tilde{v}\in D(\mathcal{L}^{2})$ and $u$ is solution of

[TABLE]

with initial condition $u(0,x)=0,\forall x\in\mathbb{R}^{d}$. Since $\mathcal{L}^{2}$ generates a strongly continuous semigroup on $B_{0}$ (see [3]), the solution of equation (26) is unique and vanishes identically, which implies that $v\equiv\tilde{v}$.

Proof of Proposition 4 Since $p_{Y}(t,0,u)=\frac{1}{\sqrt{2\pi(1-e^{-t})}}e^{-\frac{u^{2}}{2(1-e^{-t})}}$, it satisfies the two following identities,

[TABLE]

[TABLE]

The argument in the proof of Theorem 3 applies to (27) instead of the original equation (22). Thus (23) is satisfied. (23) is indeed a particular case of (19) where $p=1,q=0,r=2$, $\beta_{0,1}=0,\,\beta_{0,2}=\frac{1}{2}e^{-t},\alpha_{1}=1$. The boundary term $\mathcal{B}_{f}(t,x)$ is equal to $\frac{1}{2}e^{-t}p_{Y}(t,0,0)=\frac{e^{-t}}{2\sqrt{2\pi(1-e^{-t})}}.$

Remark. Let $r:=\inf\{k\in\mathbb{N},\,\partial^{k}_{u}\,p_{Y}(t,0,u)|_{u=0}\not\equiv 0\}$. Then the assumption $f\in D(\mathcal{L}^{n_{j}-1})$ can be weakened and replaced by $f\in D(\mathcal{L}^{n_{j}-1-r})$ by a slight modification of equality (3.2). For instance, if $0\not\in supp\,p_{Y}(t,0,\cdot)$, for all $t>0$, , then we only need $f\in B_{0}$ and equation (19) becomes

[TABLE]

Obviously if $Y$ is a non-negative process and $p_{Y}$ satisfies an equation of the type (18), then removing the absolute value, $v(t,x):=\mathbb{E}\left[f(X^{x}(Y_{t}))\right]$ is solution of (19).

3.3 General position processes

In this subsection, the position process $X^{x}$ is not necessarily markovian. Neither is the time process for which we assume (A1)-(A3). The position process $X^{x}$ is an Itô process, it can be written as

[TABLE]

where $(W(t),\mathcal{F}_{t})$ is a Wiener process, $\delta$ and $\beta$ are bounded $\mathcal{F}_{t}$-nonanticipative processes such that $\delta\delta^{T}$ is uniformly positive definite.

Theorem 4

Let us keep the assumptions of Theorem 3 on $Y$. Let $\mathcal{L}$ be a differential operator defined by

[TABLE]

with

[TABLE]

Let $f\in D(\mathcal{L}^{r-1})$ where $r$ is the highest order of the partial derivatives in $u$ appearing in (18). Define $v(t,x):=\mathbb{E}\left[f(X^{x}(|Y_{t}|))\right]$. Then $v\in D(\mathcal{L}^{r})$ and satisfies

[TABLE]

on $]0,+\infty[\times\mathbb{R}^{d}$ where $\mathcal{L}$ acts on $x$ and ${\cal B}_{f}(t,x)$ is the boundary term

[TABLE]

Proof. Under assumptions on $\delta$ and $\beta$, there exists from [15] a process $({\cal X}_{t})$ weak solution of

[TABLE]

such that $\forall\,t>0$, ${\cal X}_{t}^{x}$ and $X_{t}^{x}$ are identically distributed and coefficients $b$ and $\sigma$ are given by equations (29) and (30).

Then for all $t$, $\mathbb{E}(f(X_{t}^{x}))=\mathbb{E}(f({\cal X}^{x}_{t}))$. For $(Y_{t})$ independent of $X$ and of ${\cal X}$, we have

[TABLE]

therefore

[TABLE]

where $P^{u}$ in (35) denotes the semigroup of ${\cal X}$ at time $u$. Hence under assumptions (A1)-(A3), the result of Theorem 3 applies to ${\cal X}$, $Y$ and the function $v$ which satisfies PDE (31) with ${\cal L}$ the infinitesimal generator of ${\cal X}$.

4 Transformations of high order PDEs

In this section we consider transformations between two high order PDEs. We start with the use of intertwining diffusions to build such a transformation.

4.1 Intertwining diffusions

Consider two diffusions $(X_{t})_{t\geq 0}$ and $(U_{t})_{t\geq 0}$ which take values in $\mathbb{R}^{d}$ and $\mathbb{R}^{n}$ respectively, with respective semigroups $P_{t}^{X}$ and $P_{t}^{U}$. The diffusions $X$ and $U$ are intertwining if there exists a density function $\Lambda$ such that the operator $L$ defined by

[TABLE]

satisfies $P_{t}^{U}\,L=L\,P_{t}^{X}$ for all $t>0$.

Assume moreover that there exists a $\mathbb{R}^{d}\times\mathbb{R}^{n}$-valued diffusion process $(Z_{1},Z_{2})$ satisfying:

(i) $X=Z_{1}$ and $U=Z_{2}$ in law,

(ii) $\mathbb{E}(f(Z_{1}(0))|Z_{2}(0)=z_{2})=Lf(z_{2})$,

(iii) for any $t>0$ the random variables $Z_{2}(0)$ and $Z_{1}(t)$ are conditionally independent given $Z_{1}(0)$,

(iv) for any $t>0$ the random variables $Z_{2}(0)$ and $Z_{1}(t)$ are conditionally independent given $Z_{2}(t)$.

In this case it is proved in [21] that $\Lambda$ is a distributional solution of the hyperbolic PDE

[TABLE]

where $\mathcal{L}_{X}$ and $\mathcal{L}_{U}$ denote the infinitesimal generators of $X$ and $U$.

In the following theorem we show that intertwining preserves the structure of (3), provided that we iterate $X$ and $U$ by the same time process. Let us recall that if $X^{x}$, $Y$ and $f$ satisfy the assumptions of Theorem 3, then $v(t,x):=\mathbb{E}\left[f(X^{x}(|Y_{t}|))\right]$ is a solution of (19) with initial condition $v(0,\cdot)\equiv f(\cdot)$. We rewrite (19) as

[TABLE]

where the polynomial $Q(\frac{\partial}{\partial t},\mathcal{L}_{X})$ coincides with $\sum_{i=0}^{q}\sum_{j=1}^{r}\beta_{i,j}(-1)^{j}\frac{\partial^{i}}{\partial t^{i}}\mathcal{L}_{X}^{j}$.

Theorem 5

Let $X^{x}$ be a diffusion. Assume that $X^{x},Y$ and $f$ satisfy the assumptions of Theorem 3. Let $U^{x}$ be another diffusion independent of $Y$ such that $X$ and $U$ are intertwining. Define $g:=Lf$ and $h(t,x):=\mathbb{E}(g(U^{x}_{|Y_{t}|}))$. Then $h$ belongs to $\mathcal{D}((\mathcal{L}_{U})^{r-1})$ ($r$ is the highest order of the partial derivatives in $u$ appearing in (18)) and satisfies the PDE

[TABLE]

on $]0,+\infty[\times\mathbb{R}^{n}$ with initial condition $h(0,\cdot)\equiv g(\cdot)$ by definition.

Remark. An interesting point is that this theorem maps a solution defined on $\mathbb{R}_{+}\times\mathbb{R}^{d}$ into a solution defined on $\mathbb{R}_{+}\times\mathbb{R}^{n}$. A version without the boundary term $\mathcal{B}_{f}$ can be derived from Theorem 8 by a slight modification of the proof.

Proof of Theorem 5. From the relation $P_{t}^{U}L=LP_{t}^{X}$, we have for all $w\in\mathcal{D}(\mathcal{L}_{X})$,

[TABLE]

and using the boundedness of $L$, the right-hand side converges. Therefore so does the left-hand side which implies that $Lw\in\mathcal{D}(\mathcal{L}_{U})$. Hence, $f\in\mathcal{D}((\mathcal{L}_{X})^{r-1})$ implies $Lf\in\mathcal{D}((\mathcal{L}_{U})^{r-1})$. Moreover by definition

[TABLE]

Using first that $P^{U}L=LP^{X}$, and then the kernel $\Lambda$, we obtain

[TABLE]

We can apply Fubini Theorem which yields

[TABLE]

and by definition of $v$, we conclude that

[TABLE]

We apply $\sum_{k=1}^{p}\alpha_{k}\frac{\partial^{k}}{\partial t^{k}}$ to (40), exchange the order of $\frac{\partial^{k}}{\partial t^{k}}$ and the integral, and use (39) to obtain

[TABLE]

In order to simplify the proof, we consider $Q(M,N)=M^{i}N^{j}$ and $\mathcal{B}_{f,\mathcal{L}_{X}}(t,\rho)=(\mathcal{L}_{X})^{k}f(\rho)\alpha(t)$. We now decompose the right hand side of (41) into two parts.

[TABLE]

Taking the adjoint brings the term $({\cal L}_{X})^{*}\Lambda(\cdot,\rho)$ under the integral. By assumption (37), applied $j$-times, we conclude that

[TABLE]

Inverting the order between $(\mathcal{L}_{U})^{j}$ and $\int$, we obtain

[TABLE]

So by definition of $h$, we have shown

[TABLE]

[TABLE]

This combined with (41) and (42) gives the announced PDE for $h$ thanks to (40). In all the identities above ${\cal L}_{X}$ acts on $\rho$ whereas ${\cal L}_{U}$ acts on the other variable represented by a $\cdot$ (as in $(\mathcal{L}_{U})^{k}\Lambda(\cdot,\rho))$).

As an example let $X=X_{\alpha}$ and $U=X_{\alpha}+X_{\beta}$ with $X_{\alpha}$ and $X_{\beta}$ two independent squared-Bessel processes of respective dimension $2\alpha>0$ and $2\beta>0$ (cf. [21]). Their semigroups are intertwining with

[TABLE]

Here $B$ denotes the Beta function. The infinitesimal generators are given by

[TABLE]

Let $f$ be a real function of class $C^{2}$ with compact support so that $f\in\mathcal{D}(\mathcal{L}_{X})$ and let $W$ be a Brownian motion independent of $X$ and $U$. Then from Corollary 2, $v(t,x):=\mathbb{E}\left[f(X^{x}(|W_{t}|))\right]$ is the unique solution in $\mathcal{D}((\mathcal{L}_{X})^{2})$ of

[TABLE]

and from Theorem 5, $h(t,x):=\mathbb{E}\left[g(U^{x}(|W_{t}|))\right]$ is the unique solution in $\mathcal{D}((\mathcal{L}_{U})^{2})$ of

[TABLE]

where $g(x):=Lf(x)=\int_{0}^{1}f(x\rho)\rho^{\alpha-1}(1-\rho)^{\beta-1}d\rho.$

4.2 Mapping a high order PDE into another one

In this section, we consider two initial-boundary value problems of high order (cf.[7]). We will construct a mapping that transforms one of these problems into the other one. We show that this mapping can be expressed using a Feynman-Kac formula. If $x=(x_{1},\ldots,x_{n})$ is a point in ${\mathbb{R}}^{n}$ and $\phi$ is a smooth function, we write $D_{i}\phi(x):=\partial\phi(x)/\partial x_{i}$. We need compositions of the form $D^{\alpha}:=D_{1}^{\alpha_{1}}D_{2}^{\alpha_{2}}\cdots D_{k}^{\alpha_{k}}$ where $k$ and the $\alpha_{i}$ are integers. Finally we denote by $P(x,D)$ any polynomial in the partial derivatives as follows

[TABLE]

where $a_{\alpha}(x)$ are given functions of $x$. The first initial-boundary value problem is

[TABLE]

where $B(x,D)$ is a non tangential boundary operator and ${\cal S}:=\{x;S(x)=0\}$ for some function $S$, denotes a cylindrical surface.

Theorem 6

Suppose that (45) admits a solution $v(x,t)$. Let $B$ be a Brownian motion and suppose that the function $u(t,x):=\frac{\partial^{1/2}}{\partial t^{1/2}}\mathbb{E}(v(x,|B_{2t}|))$ is well defined. Then $u(t,x)$ solves the following problem

[TABLE]

where $f(x,t):=\frac{1}{2\sqrt{\pi}t^{3/2}}\int_{0}^{\infty}\xi e^{-\xi^{2}/4t}g(x,\xi)d\xi$.

In Theorem 6 the derivative in $u(t,x)$ is a particular case of the Caputo fractional derivative. Take a function $f$ and a positive real $\gamma$. If $m-1<\gamma<m$ for some integer $m$, the Caputo derivative of $f$ at order $\gamma$ is

[TABLE]

If $\gamma$ is an integer $\frac{\partial^{\gamma}}{\partial t^{\gamma}}f(t)$ is the usual derivative $f^{(\gamma)}$.

Proof of Theorem 6 Remember that $\Gamma(1/2)=\pi$. Then for $\gamma=1/2$,

[TABLE]

Let $w(t,x):=\mathbb{E}(v(x,|B_{2t}|))=\frac{2}{\sqrt{4\pi t}}\int_{0}^{\infty}e^{-\xi^{2}/4t}v(x,\xi)d\xi$. Then,

[TABLE]

From [7], $u(x,t)=\frac{1}{2\sqrt{\pi}t^{3/2}}\int_{0}^{\infty}\xi e^{-\xi^{2}/4t}v(x,\xi)d\xi$ is solution of (46).

4.3 Time change.

Theorem 7

Let $\alpha:\mathbb{R}\rightarrow{\mathbb{R}}^{+}$ be a positive increasing and differentiable function. Let $X$ be an $\mathbb{R}^{d}$ valued process and $Y$ a real valued process satisfying the assumptions of Section 2 such that $\sigma_{Y}$ and $b_{Y}$ in $(\ref{edsY})$ depend only on time. Let us define

[TABLE]

and $\tilde{a}_{X}:=\tilde{\sigma}_{X}\tilde{\sigma}_{X}^{T}$. Consider the iterated process $Z:=X(\alpha(Y_{t}))$. Its density $p_{Z}(t,x,z)$ satisfies

[TABLE]

where the differential operator $\Gamma$ acts on smooth functions $\varphi$ by

[TABLE]

For instance, we can take $\alpha(x):=e^{x}$.

Proof of Theorem 7 By independence of $X$ and $Y$, the density of $Z$ is given by

[TABLE]

We then have

[TABLE]

Now, since $\alpha^{\prime}(y)a^{i,j}_{X}(\alpha(y),z)=\tilde{a}^{i,j}_{X}(z)$ and $\alpha^{\prime}(y)b^{i}_{X}(\alpha(y),z)=\tilde{b}^{i}_{X}(z)$, using the forward Kolmogorov (or Fokker-Planck) equation for $p_{X}$, we are able to conclude since

[TABLE]

5 Position process indexed by the real line.

The PDEs that we have associated to iterated processes so far exhibit terms depending on the initial value (cf. Theorem 3 where $f(\cdot)\equiv v(0,\cdot)$ and some of its derivatives appear on the right-hand side of (19)). The PDEs obtained in [2] and [18] for iterated processes also have this drawback due to the use of the absolute value of the time process. Another type of PDE without this drawback is obtained in [14] but the underlying iterated process takes values in the complex plane. Extending the construction (2), we obtain in this section a Feynman-Kac formula where the initial value does not appear any longer in the PDE with a real valued underlying process.

Let us consider $(X_{+}(t))_{t\geq 0}$ and $(X_{-}(t))_{t\geq 0}$ two Markov processes with infinitesimal generators ${\cal L}_{+}$ (resp. ${\cal L}_{-}$) such that $X_{+}$ (resp. $X_{-}$) takes values in $\mathbb{R}^{d}$ (resp. $\mathbb{R}^{n}$). Inspired by Funaki’s construction [14], we define the $\mathbb{R}^{d+n}$-valued process $(X^{(x_{1},x_{2})}_{t})_{t\in\mathbb{R}}$ starting from $(x_{1},x_{2})\in\mathbb{R}^{2d}$ and given for any real time index by

[TABLE]

Note that in [14] the resulting process takes values in the complex plane whereas each component of our $X^{(x_{1},x_{2})}$ is real valued.

We will be interested in bounded functions $f:\mathbb{R}^{d}\rightarrow\mathbb{R}$ which admit an extension $\tilde{f}:\mathbb{R}^{d+n}\rightarrow\mathbb{R}$ satisfying

[TABLE]

In (iv), the operator $\mathcal{L_{+}}$ (resp. $\mathcal{L_{-}}$) acts on $x_{1}$ (resp. $x_{2}$). Let us stress the resemblance between (iv) and the intertwining identity (37) in section 4.1.

5.1 Main statement

Theorem 8

Let $X$ be defined by (47). Consider $Y$ a real valued continuous process, independent of $X_{+}$ and $X_{-}$, such that $Y_{0}=0$. Let us assume that $Y$ admits a density $p_{Y}(t,0,y)$ which satisfies

[TABLE]

for some polynomial $P$ with constant coefficients. Suppose moreover that $\frac{\partial^{k}}{\partial t^{k}}p_{Y}(t,0,\cdot)$ is integrable for all $k$ up to the order of $P$ w.r.t. its second variable.

Then

[TABLE]

is solution of the PDE

[TABLE]

with initial condition $v(0,x)=f(x)$, where $\mathcal{L}_{+}$ acts on $v$ as a function of $x$.

It can be noticed that if $X_{+}$ and $X_{-}$ are independent and $f$ satisfies assumptions (48), then for all $t>0$ the function $x\mapsto\mathbb{E}\tilde{f}(X_{+}^{x}(t),X_{-}^{0}(t))$ satisfies (48) too, with an extension given by $(x,y)\mapsto\mathbb{E}\tilde{f}(X_{+}^{x}(t),X_{-}^{y}(t))$.

Theorem 8 applies if we choose for instance $f(x):=e^{-x^{2}}$ with its extension $\tilde{f}(x,y):=e^{-x^{2}}\cdot\arctan(y)$ (cf. (48)), and independent diffusions $X_{+}$ and $X_{-}$ with infinitesimal generators $\mathcal{L}_{+}g(x)=\frac{1}{2}\partial_{xx}g(x)+x\partial_{x}g(x)$ and $\mathcal{L}_{-}g(y)=(1+y^{2})\partial_{yy}+3\partial_{y}$.

Let us mention that the assumptions can be weakened and $\tilde{f}$ can be unbounded when $X_{-}$ or $X_{+}$ is a diffusion process. For instance, let $X_{+}$, $X_{-}$, $Y$ be independent, $X_{+}$ be an Ornstein-Uhlenbeck process and $X_{-}$ and $Y$ be two Brownian motions. In this case $\mathcal{L}_{+}=\frac{1}{2}\frac{\partial^{2}}{\partial x^{2}}-x\frac{\partial}{\partial x}$ and $\mathcal{L}_{-}=\frac{1}{2}\frac{\partial^{2}}{\partial y^{2}}$. Let $f(x):=x$ extended in $\tilde{f}(x,y):=x\cosh(y)$. Then $(\mathcal{L}_{+}+\mathcal{L}_{-})\tilde{f}(x,y)=0$ and $\tilde{f}(x,0)=f(x)=x$, $\forall\,x,y\in\mathbb{R}$. Then for $t>0$,

[TABLE]

If $t<0$,

[TABLE]

Therefore

[TABLE]

Then iterating by $Y_{t}$, and setting $v(t,x):=\mathbb{E}\left[\tilde{f}(X^{(x,0)}(Y_{t}))\right]]$, we obtain $v(t,x)=\mathbb{E}\left[xe^{-Y_{t}/2}\right]=xe^{t/8}$ which indeed satisfies

[TABLE]

forall $x\in\mathbb{R}^{d}$ and $t>0$, with initial condition $v(0,x)=x$, as stated in (49).

Proof of Theorem 8. For $t\in\mathbb{R}$ and $x\in\mathbb{R}^{d}$, define $\psi(x,t):=\mathbb{E}\left[\tilde{f}(X^{(x,0)}_{t})\right]$. Remember that the notation $X^{(x,0)}$ implies that $X_{-}$ starts at [math] this is why we write $X_{-}^{0}$ below. We start by proving that

[TABLE]

If $t>0$, then $\psi(x,\cdot)=\mathbb{E}(f(X_{+}^{x}(\cdot)))$ on an open neighborhood of $t$ and therefore $\frac{\partial}{\partial t}\psi(x,t)=\mathcal{L}_{+}\psi(x,t).$

If $t<0$, for all $s$ in an open interval containing $t$, we have $\psi(x,s)=\mathbb{E}(\tilde{f}(x,X_{-}^{0}(-s)))=\mathbb{E}(g(X_{-}^{0}(-s)))=P_{-}^{-s}g\,(0)$ with $g:x_{2}\mapsto\tilde{f}(x,x_{2})$ where $x$ plays the role of a parameter. Therefore $\frac{\partial}{\partial t}\psi(x,t)=-\,P_{-}^{-t}\,{\cal L}_{-}\,g\,(0).$ The operator ${\cal L}_{-}$ acts on the second variable $x_{2}$. With the notations of (48), ${\cal L}_{-}\,g(x,x_{2})$ is also ${\cal L}_{-}\,\tilde{f}$ so using $(iv)$ of (48), we obtain $-{\cal L}_{-}\,g(x,x_{2})={\cal L}_{+}\tilde{f}(x,x_{2})$ where ${\cal L}_{+}$ acts only on the first variable $x$. We conclude that $\frac{\partial}{\partial t}\psi(x,t)=P_{-}^{-t}\,{\cal L}_{+}\,g\,(0)={\cal L}_{+}\,P_{-}^{-t}\,g\,(0)$, the latter identity being true since ${\cal L}_{+}$ acts only on $x$.

It remains to study the case $t=0$. Previously we obtained that $\frac{\partial}{\partial t}\psi(x,t)=\mathcal{L}_{+}\,P_{+}^{t}\,f\ (x)$ for $t>0\,$ and $\frac{\partial}{\partial t}\psi(x,t)=\mathcal{L}_{+}\,P_{-}^{-t}\,g\,(0)$ for $t<0\,$. The left-hand sides admit the same limit when $t\rightarrow 0$ equal to $\mathcal{L}_{+}\,f\,(x)$ since $\mathcal{L}_{+}$ acts only on $x$. Hence $t\rightarrow\psi(x,t)$ is differentiable at $t=0$ with derivative $\mathcal{L}_{+}\,f\,(x)$ which coincides with $\mathcal{L}_{+}\,\psi\,(x,0)$. Hence (50) is proved.

We now consider $v(t,x):=\mathbb{E}\left[\tilde{f}(X^{(x,0)}_{Y_{t}}\right]$. Then $v(t,x)=\int_{\mathbb{R}}p_{Y}(t,0,s)\psi(x,s)ds$ by independence of $X$ and $Y$ and the following identities hold,

[TABLE]

If we perform successive integrations by parts on each term involving $\frac{\partial}{\partial s}$ and its powers that appear in $P$, we see that ${(-1)}^{k}\,\int_{\mathbb{R}}\frac{\partial^{k}}{\partial s^{k}}p_{Y}(t,0,s)\psi(x,s)ds=\int_{\mathbb{R}}p_{Y}(t,0,s)\frac{\partial^{k}}{\partial s^{k}}\psi(x,s)ds$ for all integer $k$. Then we conclude using (50) that

[TABLE]

where ${\cal L_{+}}$ acts on $x$. In this way we have been able to separate the variables $t$ and $x$. Therefore

[TABLE]

This is (49) that we wanted to prove.

5.2 Application of Theorem 8

Corollary 6

Let us keep the assumptions and notations of Theorem 8. Furthermore, we assume $X_{-}$ is such that the set $\{s\in[0,t],X_{-}(s)=0\}$ has Lebesgue measure zero $\mathbb{P}$-almost surely for all $t>0$. Let $c_{+}:\mathbb{R}^{d}\rightarrow\mathbb{R}_{+}$ and $c_{-}:\mathbb{R}^{n}\rightarrow\mathbb{R}_{+}$ be two continuous bounded functions. Let $f$ satisfying assumptions (i) to (iii) of (48) and the relation

[TABLE]

*Set $c(x_{1},x_{2})=c_{+}(x_{1})\mathbf{1}_{x_{2}=0}-c_{-}(x_{2})\mathbf{1}_{x_{2}\neq 0}$, then

$\nu(t,x):=\mathbb{E}\left[\exp\left\{-\int_{0}^{Y_{t}}c(X^{(x,0)}_{s})ds\right\}\tilde{f}(X^{(x,0)}_{Y_{t}})\right]$ is solution of the PDE*

[TABLE]

with initial condition $\nu(0,x)=f(x)$.

In the following statement we show that the Euler-Bernoulli beam equation

[TABLE]

where $g(x)>0$ is the flexural rigidity and $m(x)>0$ the lineal mass density, can be obtained as a consequence of Theorem 8, by considering the iteration of a Brownian motion by an independent Cauchy process.

Corollary 7

Let $g$ and $m$ be two positive functions. Let $X_{+}$ be a diffusion process with infinitiesmal generator $\mathcal{L}_{+}:=g(x)\partial_{xx}^{2}$. Consider $X_{-}$ (resp. C) a Markov (resp. Cauchy) process such that $X_{+}$, $X_{-}$ and $C$ are independent. Define $\gamma_{t}:=C_{t/\sqrt{g(x)m(x)}}$. Then for any $f$ extendable in $\tilde{f}$ in the sense of (48), the function

[TABLE]

satisfies equation $(\ref{eq:47})$ with initial condition $u(0,x)=f(x),\forall\,x\in\mathbb{R}$.

As an example, let be $Y$ a Brownian motion with drift $\mu$. Corollary 6 provides a probabilistic representation of the solution to equation

[TABLE]

by $u(t,x)=\mathbb{E}\left[\exp\left\{\int_{0}^{Y_{t}}c(X^{(x,0)}(s))ds\right\}\tilde{f}(X^{(x,0)}_{Y_{t}})\right]$. A solution to such an equation can be processed using the algorithm developed in this paper and the remark following Theorem 2.

Proof of Corollary 6: Let us define $(\tilde{X}_{+}(t))_{t\geq 0}$ (resp. $(\tilde{X}_{-}(t))_{t\geq 0}$) two Markov processes with infinitesimal generators ${\cal L}_{+}-c_{+}$ (resp. ${\cal L}_{-}-c_{-}$) and the process

[TABLE]

From Theorem 8, $\nu(t,x):=\mathbb{E}\left[\tilde{f}(\tilde{X}^{(x,0)}_{Y_{t}})\right]$ is solution of

[TABLE]

For $t\geq 0$, using the definition (47) of $X^{(x_{1},x_{2})}$ and the Feynman-Kac formula,

[TABLE]

When $t<0$, the following identities hold

[TABLE]

Since $c(X^{(x,0)}(s)=-c_{-}(X_{-}^{0}(-s))$ if $X_{-}^{0}(-s)\neq 0$, using assumption on the zero set of $X_{-}$, we have

[TABLE]

So

[TABLE]

Proof of Corollary 7. Since the density of $C$, $p_{C}(t,0,y):=\mathbb{P}(C_{t}\in dy)=\frac{t}{\pi(t^{2}+y^{2})}$ satisfies

[TABLE]

$w(t,y):=p_{C}(t/\sqrt{g(x)m(x)},0,y)$ is solution of

[TABLE]

Applying Theorem 8 gives $w$ is solution of

[TABLE]

Factorizing by $g(x)$ ends the proof.

6 Application of the numerical scheme to the iterated Brownian motion

In this section, we illustrate the algorithm proposed in Section 2. First we simulate a trajectory of an iterated Brownian motion (IBM) $Z_{t}=X(|Y_{t}|)$ on $[0,T]$ where $X$ and $Y$ are two independent Brownian motions. Then we approximate numerically the function $v(t,x):=\mathbb{E}\left[f(Z^{x}_{t}))\right]$ and the variations of order three and four of $(Z_{t})$.

For a fixed positive integer $n$, the Brownian motion $Y$ is evaluated at times $kT/n$ and we define our piecewise constant processes $(\bar{Y}_{t}^{n})_{0\leq t\leq T}$ recursively by

[TABLE]

where $\xi_{k}^{n}$ are independent centered Gaussian random variables with variance $T/n$. This determines $M_{n}=\sup_{t\in[0,T]}|\bar{Y}_{t}^{n}|$ which is a.s. finite. The same construction for $\bar{X}_{t}^{n}$ on $[0,M_{n}]$ can be performed,

[TABLE]

where $(\zeta_{k}^{n})$ are independent centered Gaussian random variables with variance $M_{n}/n$.

The composition $\bar{X}^{n}(|\bar{Y}^{n}(kT/n)|)$ for $k=0,1,\ldots,n$, is given by

[TABLE]

with $\lfloor\cdot\rfloor$ the floor function. The continuous approximation $\tilde{Z}^{n}$ of $Z$ is the linear interpolation between the points defined in (54). A trajectory of $\tilde{Z}^{n}$ can be seen in Figure 2 presenting huge variations.

As mentioned in section 2 the scheme (54) is interesting since it converges uniformly as $n$ tends to infinity, it simply requires to simulate $2n$ independent Gaussian random variables and the composition of $\bar{X}^{n}$ by $|\bar{Y}^{n}|$ is facilitated by the use of step functions. This makes it possible to use methods such as Monte Carlo methods which require the simulation of thousands of trajectories (see Figure 2). Moreover the almost sure uniform convergence proved in this paper for $Z^{n}$ and its version $\tilde{Z}^{n},$ which is more convenient for implementation, renders possible all types of numerical studies requiring an approximation of a whole trajectory and not only of the value at a fixed time, like for instance the variations of various order. In Figure 3 we apply this remark to the third and fourth order variations of $(Z_{t})$, illustrating the following results of [9]

[TABLE]

[TABLE]

where $\Lambda$ denotes an arbitrary subdivision $t_{0}=0\leq t_{1}\leq\ldots\leq t_{n}=t$ of $[0,t]$ with mesh $|\Lambda|:=\max_{1\leq k\leq n}|t_{k}-t_{k-1}|$:

The algorithm can also be used to simulate the solution to a fourth-order PDE of the type we studied in the previous sections. Figure 4 shows an approximation of $v(t,x):=\mathbb{E}\left[f(Z^{x}_{t}))\right]$ corresponding to $f(x)=e^{-x^{2}}$. From Corollary 2 we know that $v(t,x)$ is the unique solution of

[TABLE]

with initial condition $v(0,x)=f(x)$.

7 Appendix

7.1 Classical results for section 2

Proposition 8

*(cf. [13]).

Let $T>0$ and $(U_{t})_{0\leq t\leq T}$ be the solution of (4) under assumptions (8)-(9). Then for all $t,s\in[0,T]$ and $p\geq 1$ such that $\mathbb{E}|U_{0}|^{2p}<\infty$,*

[TABLE]

[TABLE]

[TABLE]

for some constant $C>0$ depending only on $K$ and $p$.

Proof of Proposition 1. $C$ can denote different constants from lines to lines depending only on $K$, $p$, $n$ and $\mathbb{E}|X_{0}|^{2p}$. Keeping the notations in [12], we have for $t\in[t_{k},t_{k+1}]$:

[TABLE]

where $\epsilon_{t_{k}}^{n}$ denotes the error process defined by $\epsilon_{t}^{n}:=X_{t}-X_{t}^{n}$.

[TABLE]

From Proposition 8 and the assumptions of Section 2, we have

[TABLE]

Combining these inequalities with (60) yields

[TABLE]

for $\gamma=\max(p,2p\beta)$ if $\Delta_{t}\geq 1$ and $\gamma=\min(p,2p\beta)$ otherwise. The same inequality holds for $\sigma$. By writing $\epsilon_{s}^{n}=\epsilon_{t_{k}}^{n}+(X_{s}-X_{t_{k}})+(\tilde{X}_{t_{k}}^{n}-\tilde{X}_{s}^{n})$, and using inequality (56), we obtain

[TABLE]

Thus, inequality (59) for $t=t_{k+1}$ becomes

[TABLE]

Consequently,

[TABLE]

Using

[TABLE]

and BurkholderDavisGundy inequality, we have

[TABLE]

Using this latter inequality as well as (61) and (62) we complete the proof.

Proof of Lemma 1. We suppose $\forall\alpha<l,n^{\alpha}\sup_{t\in[0,T]}|f_{n}(t)-f(t)|\rightarrow_{n}0$. The second implication can be shown similarly. Let $\alpha>0$ and $t\in[0,T]$, then

[TABLE]

By construction,

[TABLE]

Let $t,t^{\prime}\in[0,T]$,

[TABLE]

$f$ being $\beta$-Hölder continue for all $\beta<l$, we choose $\beta$ such that $l>\beta>\alpha$ and inequality (64) becomes for some $C>0$

[TABLE]

Thus,

[TABLE]

Consequently, the right side of inequality (63) multiplied by $n^{\alpha}$ tends to [math] as $n$ tends to infinity and the result follows.

7.2 Classical results for Theorem 3

Proof of Lemma 2. We first assume that $g$ is infinitely differentiable with compact support. Let $h>0$,

[TABLE]

We have to show the convergence of the right hand side as $h\downarrow 0$ in

$(B(\mathbb{R},\mathbb{R}),{||\cdot||_{\infty}})$. Since $P^{t}f$ converges to $f$ uniformly as $t\downarrow 0$ and $g$ is continuous, $\exists\delta>0$ such that for all $s\in[0,\delta]$, $\|P^{s}f-f\|_{\infty}<\epsilon/\|g\|_{\infty}$ and $|g(s)-g(0)|<\epsilon/\|f\|_{\infty}$ so

[TABLE]

Thus, $\forall\epsilon>0$, $\exists\delta$ such that $\forall\ 0<h<\delta$,

[TABLE]

By dominated convergence theorem, using the bound $\|P^{s}f\|_{\infty}\leq\|f\|_{\infty}$ and the regularity of $g$, we have

[TABLE]

This ensures that $F(g)\in D(\mathcal{L})$ and $\mathcal{L}F(g)=-F(g^{\prime})-g(0)f$. Let $(\alpha_{n})_{n\in\mathbb{N}}$ be an infinitely differentiable mollifier. We define $g_{m}(x):=g(x)\mathbf{1}_{|x|\leq m}$ and $\tilde{g}_{m}(x):=g^{\prime}(x)\mathbf{1}_{|x|\leq m}$. Then for all $n,m\in\mathbb{N},(g_{m}\star\alpha_{n})$ is $\mathcal{C}^{\infty}$ with compact support. Since $g_{m}\in L^{1}(\mathbb{R})$ and $g^{\prime}_{m}\in L^{1}$, $(g_{m}\star\alpha_{n})\rightarrow g_{m}$ in $L^{1}(\mathbb{R})$ and $(g_{m}\star\alpha_{n})^{\prime}=(\tilde{g}_{m}\star\alpha_{n})\rightarrow\tilde{g}_{m}$ in $L^{1}(\mathbb{R})$ $g$ being continuous, $g_{m}\in L^{\infty}(\mathbb{R})$ for all $m\in\mathbb{N}$, so $|(g_{m}\star\alpha_{n})(0)-g(0)|=|(g_{m}\star\alpha_{n})(0)-g_{m}(0)|\rightarrow_{n}0$. This and the convergence of $g_{m}$ to $g$ and $\tilde{g}_{m}$ to $g^{\prime}$ in $L^{1}(\mathbb{R})$ imply that we can construct a sequence $(\phi_{n})_{n}$ of $\mathcal{C}^{\infty}$ functions with compact support such that

[TABLE]

For all $n\in\mathbb{N}$, $F(\phi_{n})\in D(\mathcal{L})$ and $\mathcal{L}F(\phi_{n})=-F(\phi_{n}^{\prime})-\phi_{n}(0)f$. From the inequality

[TABLE]

we have $\|F(\phi_{n})-F(g)\|_{\infty}\rightarrow_{n}0$ and similarly $\|F(\phi^{\prime}_{n})-F(g^{\prime})\|_{\infty}\rightarrow_{n}0$.

[TABLE]

We have shown that $(F(\phi_{n}))$ is then a sequence of $D(\mathcal{L})$ converging to $F(g)$ and $\mathcal{L}F(\phi_{n})=-F(\phi^{\prime}_{n})-\phi_{n}(0)f$ is converging to $-F(g^{\prime})-g(0)f$. Since $\mathcal{L}$ is a closed operator, we conclude that $F(g)\in D(\mathcal{L})$ and $\mathcal{L}F(g)=-F(g^{\prime})-g(0)f$. The second part of the Lemma is proved by induction.

Lemma 3

Let $g$ be a continuous density function. For all $f$ in the domain of $\mathcal{L}$,

[TABLE]

Proof of Lemma 3. Let $h>0$ and $p(h,x,dy):=\mathbb{P}(X_{h}^{x}\in dy)$

[TABLE]

So

[TABLE]

where $I$ stands for the identity. Since, $(P^{h+s}f(x)-P^{s}f(x))/h\rightarrow\mathcal{L}P^{s}f(x)$ uniformly in $x$ as $h\downarrow 0$, letting $h$ decreasing to [math] in equality (66) ends the proof.