Backward Stochastic Differential Equations Driven by Fractional Brownian Motion: Theory and Applications

Bou Diop

Department of Mathematics and Computer Science, Faculty of Technology and Computer Science, University Iba Der Thiam of Thiès, Thiès, Senegal.

DOI: 10.4236/jamp.2025.1311212                

Abstract

This article develops the theory of backward stochastic differential equations (BSDEs) governed by a fractional Brownian motion with Hurst parameter $H\in \left(1/2,1\right)$ . We establish existence, uniqueness, and regularity results for solutions in appropriate Sobolev spaces. In particular, we examine the behavior of solutions depending on the value of the parameter $H$ . We also provide an application to stochastic optimal control through a precise formulation of the maximum principle. Numerical simulations illustrate the dependence of solutions on model parameters.

Keywords

BSDE, Fractional Brownian Motion, Malliavin Calculus, Optimal Control, Long-Range Dependence

1. Introduction

Backward stochastic differential equations (BSDEs) were introduced by Pardoux and Peng [1] to solve stochastic control problems and financial asset pricing. Fractional Brownian motion, whose increments are correlated, allows modeling phenomena with long memory, such as certain financial time series or physical processes [2] [3].

We have three main contributions:

1) Establishment of sufficient conditions for existence and uniqueness for fractional BSDEs.

2) Analysis of solution regularity via Malliavin calculus.

3) Application to optimal control of dynamical systems perturbed by fractional noise.

The option to work with backward stochastic differential equations (BSDEs) was favored because they allow for the natural integration of terminal conditions and the formulation of the maximum principle in a fractional framework, which is not achievable with strictly forward SDEs.

2. Mathematical Framework and Notations

We consider the complete probability space generated by the Brownian motion $B$ , namely: $\left(\Omega ,ℱ,{\left\{{ℱ}_t\right\}}_{t\ge 0},\mathbb{P}\right)$ where $B^H$ is a fractional Brownian motion with parameter $H\in \left(1/2,1\right)$ adapted to the filtration $\left\{{ℱ}_t\right\}$ . When $H>1/2$ , the covariance takes the form:

$\mathbb{E}\left[B_t^H{B}_s^H\right]=\frac{1}{2}\left(t^{2H}+s^{2H}-{\left|t-s\right|}^{2H}\right)$

We define the stochastic integral in the Skorokhod sense. For an adapted process $u$ , the isometry is written as:

$\mathbb{E}\left[{\left({\displaystyle {\int }_0^T{u}_s\text{d}{B}_s^H}\right)}^2\right]=H\left(2H-1\right){\displaystyle {\int }_0^T{\displaystyle {\int }_0^T\mathbb{E}}}\left[u_s{u}_t\right]{\left|s-t\right|}^{2H-2}\text{d}s\text{d}t$

The studied BSDE is:

$Y_t=\xi +{\displaystyle {\int }_t^{T}f\left(s,Y_s,Z_s\right)\text{d}s}-{\displaystyle {\int }_t^T{Z}_s\text{d}{B}_s^H}$ (1)

with:

• $\xi :\Omega \to ℝ$ is ${ℱ}_T$ -measurable (terminal condition).

• $f:\Omega \times \left[0,T\right]\times ℝ\times ℝ\to ℝ$ is progressively measurable.

• $\left(Y,Z\right)$ is the adapted solution pair with values in $ℝ\times ℝ$ .

The Hurst coefficient $H$ regulates memory: $H>0.5$ $\to$ persistence (long memory), $H=0.5$ $\to$ classical Brownian motion, $H<0.5$ $\to$ anti-persistence. This is why it is used to model phenomena with extended temporal dependence [4] [5].

3. Existence and Uniqueness of the Solution

Theorem 3.1 (Existence and uniqueness) Under the hypotheses:

(H1) $\xi \in L^2\left(\Omega ,{ℱ}_T,\mathbb{P}\right)$

(H2) $\exists K>0$ such that $\forall t\in \left[0,T\right],\forall y,y^{\prime },z,z^{\prime }$ :

$\left|f\left(t,y,z\right)-f\left(t,y^{\prime },z^{\prime }\right)\right|\le K\left(\left|y-y^{\prime }\right|+\left|z-z^{\prime }\right|\right)$

(H3) $\mathbb{E}\left[{\displaystyle {\int }_0^T{\left|f\left(s,0,0\right)\right|}^2\text{d}s}\right]<\infty$

Then the BSDE (1) admits a unique solution $\left(Y,Z\right)\in S^2\left(0,T\right)\times H_H^2\left(0,T\right)$ where:

$S^2\left(0,T\right)=\left\{Y\text{adapted}\text{continuous}:\mathbb{E}\left[\underset{t\in \left[0,T\right]}{\text{sup}}{\left|Y_t\right|}^2\right]<\infty \right\}$

$H_H^2\left(0,T\right)=\left\{Z\text{adapted}:\mathbb{E}\left[{\displaystyle {\int }_0^T{\left|Z_s\right|}^2\text{d}s}+{\displaystyle {\int }_0^T{\displaystyle {\int }_0^T\left|Z_s{Z}_t\right|{\left|s-t\right|}^{2H-2}\text{d}s\text{d}t}}\right]<\infty \right\}$

Proof. The proof is established by the fixed point method in the appropriate Banach space. Define the operator $\Phi :S^2\left(0,T\right)\times H_H^2\left(0,T\right)\to \left(U,V\right)\mapsto \left(Y,Z\right)$ in the following sense:

$Y_t=\mathbb{E}\left[\xi +{\displaystyle {\int }_t^{T}f}\left(s,U_s,V_s\right)\text{d}s|{ℱ}_t\right]$

with $Z$ determined by the martingale representation theorem.

Step 1: Observation.

In the fractional context, even though $B^H$ is not a martingale, a Clark-Ocone-type representation remains accessible within the framework of the Skorokhod integral (refer to Biagini et al. [4], Nualart [5]). This characteristic is essential, as it ensures the existence of a process $Z$ . By verifying the stochastic decomposition of the solution, it facilitates the efficient use of the fixed-point argument.

Preliminary estimate: Let $\left(Y,Z\right)=\Phi \left(U,V\right)$ . We have the classical estimate:

$\begin{array}{c}\mathbb{E}\left[\underset{t\in \left[0,T\right]}{\text{sup}}{\left|Y_t\right|}^2\right]\le C_1\mathbb{E}\left[{\left|\xi \right|}^2+{\displaystyle {\int }_0^T{\left|f\left(s,U_s,V_s\right)\right|}^2\text{d}s}\right]\\ \le C_2\left(\mathbb{E}\left[{\left|\xi \right|}^2\right]+\mathbb{E}\left[{\displaystyle {\int }_0^T{\left|f\left(s,0,0\right)\right|}^2\text{d}s}\right]+K^2{\displaystyle {\int }_0^T\mathbb{E}}\left[{\left|U_s\right|}^2+{\left|V_s\right|}^2\right]\text{d}s\right)\end{array}$

Regarding the term $Z$ , we apply the fractional isometry:

$\mathbb{E}\left[{\displaystyle {\int }_0^T{\left|Z_s\right|}^2\text{d}s}\right]\le C\left(H\right)\mathbb{E}\left[{\left|\xi \right|}^2+{\displaystyle {\int }_0^T{\left|f\left(s,U_s,V_s\right)\right|}^2\text{d}s}\right]$

Step 2: Contraction.

Consider two pairs $\left(U^1,V^1\right)$ and $\left(U^2,V^2\right)$ in $S^2\left(0,T\right)\times H_H^2\left(0,T\right)$ and denote $\left(Y^1,Z^1\right)=\Phi \left(U^1,V^1\right)$ , $\left(Y^2,Z^2\right)=\Phi \left(U^2,V^2\right)$ . We then have:

$\begin{array}{c}\left|Y_t^1-Y_t^2\right|\le \mathbb{E}\left[{\displaystyle {\int }_t^T\left|f\left(s,U_s^1,V_s^1\right)-f\left(s,U_s^2,V_s^2\right)\right|\text{d}s}|{ℱ}_t\right]\\ \le K\mathbb{E}\left[{\displaystyle {\int }_t^T\left(\left|U_s^1-U_s^2\right|+\left|V_s^1-V_s^2\right|\right)\text{d}s}|{ℱ}_t\right]\end{array}$

Taking the expectation of the supremum and using Doob’s inequality:

$\mathbb{E}\left[\underset{t\in \left[0,T\right]}{\text{sup}}{\left|Y_t^1-Y_t^2\right|}^2\right]\le C_3{K}^{2}T\mathbb{E}\left[{\displaystyle {\int }_0^T\left({\left|U_s^1-U_s^2\right|}^2+{\left|V_s^1-V_s^2\right|}^2\right)\text{d}s}\right]$

For the $Z$ term, we again use the isometry:

$\mathbb{E}\left[{\displaystyle {\int }_0^T{\left|Z_s^1-Z_s^2\right|}^2\text{d}s}\right]\le C_4{K}^{2}T\mathbb{E}\left[{\displaystyle {\int }_0^T\left({\left|U_s^1-U_s^2\right|}^2+{\left|V_s^1-V_s^2\right|}^2\right)\text{d}s}\right]$

Thus, for $T$ sufficiently small such that $C_5{K}^{2}T<1$ , Φ is indeed a contraction.

Step 3: Extension to the entire interval.

When $T$ is not sufficiently small, we partition the interval $\left[0,T\right]$ into $n$ subintervals $\left[t_i,t_{i+1}\right]$ , $t_0=0$ , $t_n=T$ with $t_{i+1}-t_i$ small enough so that Φ is a contraction on each subinterval. We then solve the BSDE on each interval taking as terminal condition the one from the previous interval.

Step 4: Uniqueness. Uniqueness is an immediate consequence of the contraction principle. If $\left(Y^1,Z^1\right)$ and $\left(Y^2,Z^2\right)$ are solutions, then:

$\Vert \left(Y^1,Z^1\right)-\left(Y^2,Z^2\right)\Vert \le C_5{K}^{2}T\Vert \left(Y^1,Z^1\right)-\left(Y^2,Z^2\right)\Vert$

which implies $\left(Y^1,Z^1\right)=\left(Y^2,Z^2\right)$ if $C_5{K}^{2}T<1$ . For an arbitrary interval $T$ , we partition again. □

4. Regularity of the Solution and Malliavin Calculus

Theorem 4.1 (Hölder regularity). If in addition $\xi \in {\mathbb{D}}^{1,2}$ and $f$ is continuously differentiable with bounded derivatives, then:

1) $Y_t\in {\mathbb{D}}^{1,2}$ for all $t\in \left[0,T\right]$ .

2) $Z_t=D_t{Y}_t$ a.e. where $D$ is the Malliavin derivative.

3) $Y$ admits a $\gamma$ -Hölder continuous version for any $\gamma <H$ .

The Malliavin derivative satisfies the BSDE:

$D_t{Y}_s=D_t\xi +{\displaystyle {\int }_s^T\left[{\partial }_{y}f\left(r\right)D_t{Y}_r+{\partial }_{z}f\left(r\right)D_t{Z}_r\right]\text{d}r}-{\displaystyle {\int }_s^T{D}_t{Z}_r\text{d}{B}_r^H}$

Proof. Proof structure

We want to show, under the hypotheses:

• $\xi \in {\mathbb{D}}^{1,2}$ (stochastic Sobolev space).

• $f$ continuously differentiable with bounded derivatives.

that the solution $\left(Y,Z\right)$ of the BSDE possesses additional regularity.

Malliavin differentiability of Y

We begin by approximating $\xi$ and $f$ by regular sequences ${\xi }_n$ and $f_n$ . For each approximation, the corresponding solution $\left(Y_n,Z_n\right)$ is Malliavin differentiable. By the BSDE equation:

$Y_n\left(t\right)={\xi }_n+{\displaystyle {\int }_t^T{f}_n\left(s,Y_n\left(s\right),Z_n\left(s\right)\right)\text{d}s}-{\displaystyle {\int }_t^T{Z}_n\left(s\right)\text{d}{B}_s^H}$

Applying the Malliavin derivative operator $D$ which commutes with the Skorokhod integral, we obtain:

$D_r{Y}_n\left(t\right)=D_r{\xi }_n+{\displaystyle {\int }_t^T\left[{\partial }_y{f}_n{D}_r{Y}_n\left(s\right)+{\partial }_z{f}_n{D}_r{Z}_n\left(s\right)\right]\text{d}s}-{\displaystyle {\int }_t^T{D}_r{Z}_n\left(s\right)\text{d}{B}_s^H}$

We observe that $D_r{Y}_n$ satisfies a linear BSDE.

Limit passage

Under our regularity hypotheses, we show that:

$Y_n\to Y\text{in}{L}^2\left(\Omega \times \left[0,T\right]\right)$

$D_r{Y}_n\to D_{r}Y\text{in}{L}^2\left(\Omega \times {\left[0,T\right]}^2\right)$

$Z_n\to Z\text{in}{L}^2\left(\Omega \times \left[0,T\right]\right)$

$D_r{Z}_n\to D_{r}Z\text{in}{L}^2\left(\Omega \times {\left[0,T\right]}^2\right)$

The limit then verifies the BSDE:

$D_{r}Y\left(t\right)=D_r\xi +{\displaystyle {\int }_t^T\left[{\partial }_{y}f{D}_{r}Y\left(s\right)+{\partial }_{z}f{D}_{r}Z\left(s\right)\right]\text{d}s}-{\displaystyle {\int }_t^T{D}_{r}Z\left(s\right)\text{d}{B}_s^H}$

Identification $Z_t=D_t{Y}_t$

By the Clark-Ocone representation (adapted to the fractional case), we have:

$Y\left(t\right)=\mathbb{E}\left[Y\left(t\right)\right]+{\displaystyle {\int }_0^T\mathbb{E}}\left[D_{s}Y\left(t\right)|{ℱ}_s\right]\text{d}{B}_s^H$

But from the original BSDE, we also have:

$Y\left(t\right)=\mathbb{E}\left[\xi +{\displaystyle {\int }_t^{T}f}\left(s,Y\left(s\right),Z\left(s\right)\right)\text{d}s\right]-{\displaystyle {\int }_0^{T}Z}\left(s\right)\text{d}{B}_s^H$

By uniqueness of the representation, we identify:

$Z\left(s\right)=-\mathbb{E}\left[D_{s}Y\left(t\right)|{ℱ}_s\right]\text{for}s\le t$

A continuity argument allows us to conclude that $Z_t=D_t{Y}_t$ a.e.

Hölder-type regularity

To establish this regularity, we adopt the estimate:

$\mathbb{E}\left[{\left|Y\left(t\right)-Y\left(s\right)\right|}^2\right]\le C{\left|t-s\right|}^{2H}$

coming from:

• The form of the BSDE.

• The regularity of fractional Brownian motion.

• The estimates concerning $f$ and $\xi$ .

More precisely, we demonstrate that:

$\mathbb{E}\left[{\left|{\displaystyle {\int }_t^{T}Z}\left(r\right)\text{d}{B}^H\left(r\right)\right|}^2\right]\le C{\left|t-s\right|}^{2H}\left(\text{by fractional isometry}\right)$

$\mathbb{E}\left[{\left|{\displaystyle {\int }_t^{T}f}\left(r,Y\left(r\right),Z\left(r\right)\right)\text{d}r\right|}^2\right]\le C{\left|t-s\right|}^2\left(f\text{is}\text{Lipschitz}\right)$

Which in combination with $\gamma <H$ gives us the result.

Summary of results

All these combined arguments demonstrate that:

• $Y_t\in {\mathbb{D}}^{1,2}$ for all $t$ .

• $Z_t=D_t{Y}_t$ a.e.

• $Y$ is $\gamma$ -Hölder for $\gamma <H$ .

• The Malliavin derivative satisfies the aforementioned linear BSDE.

□

5. Application to Optimal Control

Consider the dynamical system:

$\{\begin{array}{l}\text{d}{X}_t=b\left(t,X_t,u_t\right)\text{d}t+\sigma \left(t,X_t,u_t\right)\text{d}{B}_s^H\hfill \\ X_0=x_0\hfill \end{array}$ (2)

with cost:

$J\left(u\right)=\mathbb{E}\left[{\displaystyle {\int }_0^{T}g}\left(t,X_t,u_t\right)\text{d}t+h\left(X_T\right)\right]$

Assumptions (H*).

To properly ensure the adjoint BSDE (3), we assume that:

• The functions $b\left(t,x,u\right)$ , $\sigma \left(t,x,u\right)$ , and $g\left(t,x,u\right)$ are measurable and Lipschitz with respect to $x$ , uniformly in $t$ and $u$ .

• (A2) They are bounded, as are their partial derivatives with respect to $x$ and $u$ .

• (A3) The terminal cost function $h\left(x\right)$ belongs to the class $C^1$ and has at most polynomial growth.

In this context, the associated BSDE admits a unique solution in the space $S^2\left(0,T\right)\times H^2\left(0,T\right)$ .

Theorem 5.1 (Stochastic maximum principle). Let $u^{*}$ be an optimal control and $X^{*}$ the associated optimal trajectory. Then there exists an adapted process $\left(p,q\right)$ solution of the adjoint BSDE:

$\{\begin{array}{l}\text{d}{p}_t=-{\partial }_{x}H\left(t,X_t^{\text{*}},u_t^{\text{*}},p_t,q_t\right)\text{d}t+q_t\text{d}{B}_t^H\hfill \\ p_T={\partial }_{x}h\left(X_T^{\text{*}}\right)\hfill \end{array}$ (3)

where the Hamiltonian is defined by:

$H\left(t,x,u,p,q\right)=g\left(t,x,u\right)+b\left(t,x,u\right)p+\sigma \left(t,x,u\right)q$

Moreover, we have the minimization condition:

$H\left(t,X_t^{\text{*}},u_t^{\text{*}},p_t,q_t\right)=\underset{u\in U}{\text{min}}H\left(t,X_t^{\text{*}},u,p_t,q_t\right)$ a.e.

Proof. The proof of the stochastic maximum principle for fractional Brownian motion follows the classical steps but requires technical adaptations due to the non-martingale nature of the noise [6] [7].

Step 1: Control variation. Let $u^{\text{*}}$ be the optimal control and $u^{ϵ}=u^{\text{*}}+ϵv$ a variation ($v$ adapted, $ϵ>0$ small). We set the trajectory variation:

$X_t^{ϵ}=X_t^{\text{*}}+ϵY_t+o\left(ϵ\right)$

where $Y_t$ satisfies the linearized BSDE:

$\text{d}{Y}_t=\left[b_x\left(t\right)Y_t+b_u\left(t\right)v_t\right]\text{d}t+\left[{\sigma }_x\left(t\right)Y_t+{\sigma }_u\left(t\right)v_t\right]\text{d}{B}_t^H$

$Y_0=0$

Step 2: Cost variation. The cost variation is written as:

$J\left(u^{ϵ}\right)-J\left(u^{\text{*}}\right)=ϵ\delta J+o\left(ϵ\right)$

with:

$\delta J=\mathbb{E}\left[{\displaystyle {\int }_0^T\left(g_x\left(t\right)Y_t+g_u\left(t\right)v_t\right)\text{d}t}+h_x\left(X_T^{\text{*}}\right)Y_T\right]$

Step 3: Introduction of adjoint variables. We introduce the adjoint process $\left(p_t,q_t\right)$ solution of the adjoint BSDE:

$\text{d}{p}_t=-\left[b_x\left(t,X_t^{\text{*}},u_t^{\text{*}}\right)p_t+{\sigma }_x\left(t,X_t^{\text{*}},u_t^{\text{*}}\right)q_t+g_x\left(t,X_t^{\text{*}},u_t^{\text{*}}\right)\right]\text{d}t+q_t\text{d}{B}_t^H$

$p_T=h_x\left(X_T^{\text{*}}\right)$

Step 4: Itô’s formula for $p_t{Y}_t$ . Using Itô’s formula for fractional Brownian motion:

$\text{d}\left(p_t{Y}_t\right)=p_t\text{d}{Y}_t+Y_t\text{d}{p}_t+\text{d}{\left[p,Y\right]}_t$

The bracket $\text{d}{\left[p,Y\right]}_t$ depends on the regularity of the processes and the index $H$ .

Step 5: Variation calculation. Integrating and taking expectation:

$\begin{array}{l}\mathbb{E}\left[p_T{Y}_T\right]=\mathbb{E}\left[{\displaystyle {\int }_0^T{p}_t}\left(b_u\left(t\right)v_t\right)\text{d}t+{\displaystyle {\int }_0^T{Y}_t}\left(-b_x\left(t\right)p_t-{\sigma }_x\left(t\right)q_t-g_x\left(t\right)\right)\text{d}t\right]\\ +\text{bracket terms}\end{array}$

The bracket terms are expressed in terms of $H$ and must be treated separately depending on whether $H>1/2$ or $H<1/2$ .

Step 6: Simplification. Using $p_T=h_x\left(X_T^{\text{*}}\right)$ , we obtain:

$\mathbb{E}\left[h_x\left(X_T^{\text{*}}\right)Y_T\right]=\mathbb{E}\left[{\displaystyle {\int }_0^T\left(p_t{b}_u\left(t\right)v_t-Y_t\left(b_x\left(t\right)p_t+{\sigma }_x\left(t\right)q_t+g_x\left(t\right)\right)\right)\text{d}t}\right]+R\left(H\right)$

where $R\left(H\right)$ represents additional terms due to the fractional nature.

Step 7: Optimality condition. Since $u^{\text{*}}$ is optimal, $\delta J\ge 0$ for any perturbation $v$ , so:

$\mathbb{E}\left[{\displaystyle {\int }_0^T\left(g_u\left(t\right)+p_t{b}_u\left(t\right)+q_t{\sigma }_u\left(t\right)\right)v_t\text{d}t}\right]\ge 0$

This means that for almost every $t$ :

$g_u\left(t,X_t^{\text{*}},u_t^{\text{*}}\right)+p_t{b}_u\left(t,X_t^{\text{*}},u_t^{\text{*}}\right)+q_t{\sigma }_u\left(t,X_t^{\text{*}},u_t^{\text{*}}\right)=0\text{a}.\text{e}.$

Step 8: Minimization condition. The last equality can also be expressed as:

$H\left(t,X_t^{\text{*}},u_t^{\text{*}},p_t,q_t\right)=\underset{u\in U}{\text{min}}H\left(t,X_t^{\text{*}},u,p_t,q_t\right)\text{a}.\text{e}.$

where the Hamiltonian is given by:

$H\left(t,x,u,p,q\right)=g\left(t,x,u\right)+b\left(t,x,u\right)p+\sigma \left(t,x,u\right)q$

Important remarks

• Fractional case: The complete proof requires an in-depth analysis of fractional stochastic calculus, particularly the manipulation of the bracket which differs from the standard Brownian case.

• Regularity: For $H>1/2$ , we often use Malliavin calculus, while for $H<1/2$ , approximation techniques are employed.

• Existence: The existence and uniqueness of the adjoint backward stochastic differential equations must be demonstrated separately.

□

6. Application to Optimal Control of Systems Driven by Fractional Brownian Motion

Backward stochastic differential equations (BSDEs) find a natural application in stochastic optimal control theory via the dynamic programming principle. Consider a dynamical system whose state $X_t\in {ℝ}^n$ is governed by an SDE driven by a fractional Brownian motion (fBm) $B_t^H$ with Hurst exponent $H\in \left(0,1\right)$ :

$\{\begin{array}{l}\text{d}{X}_t=b\left(t,X_t,u_t\right)\text{d}t+\sigma \left(t,X_t,u_t\right)\text{d}{B}_t^H\hfill \\ X_0=x_0\hfill \end{array}$ (4)

where $u_t\in U$ is an admissible control, belonging to a class of processes adapted to the fBm filtration. Our goal is to minimize the following cost criterion:

$J\left(u\right)=\mathbb{E}\left[{\displaystyle {\int }_0^{T}g\left(t,X_t,u_t\right)\text{d}t}+h\left(X_T\right)\right]$

Under standard regularity conditions on the coefficients $b$ , $\sigma$ , $g$ and $h$ ,

the optimal cost $V\left(t,x\right)={\inf }_u{J\left(u\right)|}_{X_t=x}$ is the solution of a stochastic Hamilton-Jacobi-Bellman (HJB) equation. The solution pair $\left(Y_t,Z_t\right)$ of the associated BSDE:

$\{\begin{array}{l}-\text{d}{Y}_t=f\left(t,Y_t,Z_t\right)\text{d}t-Z_t\text{d}{B}_t^H\hfill \\ Y_T=h\left(X_T\right)\hfill \end{array}$ (5)

plays a crucial role. In particular, we have the relation $Y_t=V\left(t,X_t\right)$ , and the process $Z_t$ is related to the derivative of the optimal value. The Hamiltonian function $f$ is defined by:

$f\left(t,x,y,z\right)=\underset{u\in U}{\text{inf}}\left\{g\left(t,x,u\right)+z\cdot {\sigma }^{-1}\left(t,x,u\right)b\left(t,x,u\right)\right\}$

Thus, the numerical resolution of the BSDE (5) allows us not only to evaluate the optimal cost $Y_0=V\left(0,x_0\right)$ but also to characterize the optimal control policy $u_t^{\text{*}}$ via the diffusion term $Z_t$ .

Challenges and Difficulties Encountered listing the obstacles (non-martingale, Skorokhod, Malliavin, simulation, sensitivity to $H$ ).

7. Numerical Simulations

In this section, we define a numerical analysis to verify the proposed discretization scheme and demonstrate the effect of the Hurst exponent $H$ on the solution of the BSDE. We consider an academic test case for which the analytical solution is known, allowing us to precisely quantify the discretization error [8] [9].

7.1. Model and Parameters

We consider the following linear BSDE:

$\{\begin{array}{l}-\text{d}{Y}_t=\left[-\beta Y_t+\gamma Z_t\right]\text{d}t-Z_t\text{d}{B}_t^H\hfill \\ Y_T=\text{sin}\left(B_T^H+T\right)\hfill \end{array}$

with parameters $T=1$ , $\beta =0.2$ , and $\gamma =0.3$ . For this problem, the exact solution is given by:

$Y_t={\text{e}}^{-\beta \left(T-t\right)}\mathbb{E}\left[\text{sin}\left(B_T^H+T\right)|{ℱ}_t\right]$

Explanation of the fBm Simulation Algorithm

Simulation of fBm Increments

To generate correlated increments of fractional Brownian motion, we use the Davies-Harte method (see [8]), based on the circular embedding of the covariance matrix. This method provides accurate and fast simulation of trajectories, with a complexity of $O\left(N\log N\right)$ . We have verified that the results match those obtained using the Cholesky method, although the latter involves a cost of order $O\left(N^2\right)$ . We therefore prefer the Davies-Harte method for our numerical experiments.

7.2. Discretization Scheme

For a regular partition of $\left[0,T\right]$ with $N$ time steps ($\Delta t=T/N$ ), the pair $\left(Y_{t_i},Z_{t_i}\right)$ is approximated by $\left(Y_i,Z_i\right)$ according to the following scheme:

$Y_i=\mathbb{E}\left[Y_{i+1}|{ℱ}_{t_i}\right]+\Delta tf\left(t_i,Y_i,Z_i\right)$

$Z_i=\frac{1}{\Delta t}\mathbb{E}\left[Y_{i+1}\Delta B_{t_i}^H|{ℱ}_{t_i}\right]$

The conditional expectations are predicted using a least squares regression methodology based on a set of Legendre polynomials.

7.3. Results and Analysis

Figure 1(b) shows a typical trajectory of the solution $Y_t$ for $H=0.7$ . We observe that the process follows a rather smooth dynamics, with equally continuous trends, a typical characteristic of long-memory noise. The estimated solution fits exactly with the reference solution, demonstrating the accuracy of the scheme.

On the other hand, Figure 1(a) for $H=0.3$ corresponds to a much more erratic trajectory, with frequent oscillations and anti-persistence. Our scheme can also reproduce this more realistic behavior, even though the instantaneous error is slightly larger due to the greater volatility.

(a) Simulation for $H=0.3$ (b) Simulation for $H=0.7$

Figure 1. Simulations for different H values.

Table 1 presents the accuracy of the scheme and its computational effort as a function of $H$ . The $L^2$ error is calculated with respect to a reference solution obtained using a very small time step ($\Delta t={10}^{-4}$ ). We observe that the error decreases as $H$ increases. This type of behavior is expected because increasing $H$ improves the regularity of the driver trajectories, making the problem intrinsically more stable and easier to approximate.

The computation time also increases with $H$ , reflecting the greater complexity of exactly simulating correlated increments for strong long memory.

Table 1. Discretization error and computation time as a function of the Hurst exponent $H$ ($\Delta t=0.01$ , $T=1$ ).

$H$	0.6	0.7	0.8	0.9
$L^2$ Error	0.15	0.12	0.08	0.05
Computation time (s)	45	52	68	85

7.4. Discussion

The numbers obtained in the presented results confirm the viability and effectiveness of the solution of BSDEs driven by fractional Brownian motion. The detailed analysis of the simulations highlights several key elements:

Influence of the Hurst parameter on solution regularity

The strong dependence of solutions on the parameter $H$ highlights the crucial importance of an exact representation of noise in practical applications. For $H>0.5$ , the persistence of trajectories results in greater regularity of solutions, thus facilitating their numerical approximation. On the other hand, for $H<0.5$ , anti-persistence creates additional complexity resulting from the more erratic nature of trajectories and the need for more sophisticated numerical schemes.

Performance of the numerical scheme

The decrease in $L^2$ error observed for large values of $H$ (Table 1) is a particularly encouraging result. This improvement in accuracy is justified by:

• The better regularity of fBm increments for high $H$ .

• The reduction of covariance of regression estimators.

• The greater numerical stability of approximation methods.

Comparison with related works

The attempt of this paper that discusses fractional BSDEs is still evolving in publications. The focus here is on the theoretical aspects as well as the practical computation aspects.

• Zhang (2017) [9] provided existence and uniqueness of certain BSDEs with fractional Brownian motion and established some of its qualitative properties. His work focused on the analytical aspects of the problem. Theory suggests that one can augment such works with the Malliavin Regularity analysis ($Z_t=D_t{Y}_t$ ) as well as an explicit numerical scheme and convergence tests. These results bridge the gap between theory and practice.

• Nostratipour and Hamdi (2022) [10] proposed numerical approaches based on polynomial chaos expansions. These approaches are highly accurate but very expensive and unfriendly in high dimensions. In the high dimensional cases, my polynomial based (Legendre, Hermite) regression approach is easier and much less expensive, particularly with nonliner generators. More importantly, it is applicable to general problems and more practical than other methods.

• Bender and Zhang (2008) [8] investigate time discretization methods specific to fBm.

Computational complexity

The fact that computation time increases with $H$ reflects the increase in complexity of exact simulation of correlated increments. This complexity is essentially due to:

• The need to compute and store large covariance matrices

• The increasing number of terms in polynomial expansions

• The management of long memory requiring extended time windows

Comparison with Other Numerical Techniques

To situate our model relative to current works, it is relevant to mention recently developed numerical techniques for fractional BSDEs. For instance, Nostratipour and Hamdi [10] proposed a model based on polynomial chaos expansions, while Bender and Zhang [8] explored time discretization methods specific to fBm. In comparison, our approach using polynomial regression on simulated increments offers an interesting trade-off:

• A convergence rate comparable to that mentioned in [10].

• But at a reduced computational cost thanks to the use of least squares regression.

• A more adaptable implementation for nonlinear generators.

This comparison highlights the competitiveness of our method, while leaving the door open to potential hybridizations with multi-scale or chaos-based techniques.

Implications for optimal control

These results offer promising horizons for the application of these methods to optimal control problems where exogenous perturbations have long memory, such as:

• Portfolio management with long-range dependence in returns.

• Control of physical systems with hysteresis phenomena.

• Network optimization with self-similar traffic.

• Control of chemical processes with slow relaxation.

Current limitations and challenges

Despite promising performances, some challenges persist:

• Hypersensitivity to regularity assumptions of coefficients.

• Computational complexity for high-dimensional problems.

• Accurate estimation of the parameter $H$ in practice.

• Adaptation to situations where $H$ may vary over time.

Improvement perspectives

Some research directions for improvement seem particularly promising:

• Multi-scale schemes based on fractal structure.

• Integration of machine learning techniques for approximation of conditional expectations.

• Adaptation of variance reduction techniques specific to the fractional context.

• Extension to cases of non-linear coefficients and irregular terminal conditions.

These results confirm that fractional BSDEs constitute a rich theoretical framework and offer effective numerical tools for modeling and optimizing systems subjected to long-memory perturbations. Mastering these techniques opens the way to numerous applications in domains where classical Brownian hypotheses prove insufficient.

Convergence and Robustness Analysis

Robustness and Stability of the Scheme

We studied the sensitivity of the numerical approach based on the choice of basis functions (Legendre polynomials vs. Hermite polynomials), as well as regarding discretization parameters ($\Delta t$ , $N$ ).

The results indicate that $L^2$ -norm convergence is robust with respect to the choice of polynomials, although higher degrees lead to increased computational cost without significant improvement in accuracy.

The scheme demonstrates robustness over a wide range of $H$ values, including in the anti-persistent regime ($H<0.5$ ).

Time-step refinement tests demonstrate linear convergence in $\Delta t$ , which is consistent with theoretical predictions.

This research supports the validity of our method by proving that it remains stable and accurate even under parameter variations.

8. Conclusions and Perspectives

8.1. Summary of Contributions

This study has established a complete theory of backward stochastic differential equations driven by fractional Brownian motion. These conclusions follow and extend the research of Zhang [9], Bender-Zhang [8], and Nostratipour-Hamdi [10], by introducing a Malliavin regularity analysis and a numerical validation based on the Hurst parameter. Our main contributions can be summarized as follows:

1) Theoretical foundations: We demonstrated existence and uniqueness results for solutions in adapted Sobolev spaces, establishing general sufficient conditions on coefficients and the terminal condition.

2) Regularity analysis: Through Malliavin calculus, we characterized the Hölder regularity of solutions and established the fundamental relation $Z_t=D_t{Y}_t$ , thus linking the diffusion process to stochastic derivatives.

3) Application framework: We formulated and proved a stochastic maximum principle for optimal control problems, providing necessary optimality conditions for systems perturbed by long-memory noise.

4) Numerical validation: The development and implementation of an efficient numerical scheme allowed us to precisely quantify the impact of the Hurst parameter on solution behavior and to experimentally validate our theoretical results.

8.2. Theoretical and Practical Implications

The obtained results present several significant implications.

8.2.1. On the Theoretical Level

• Extension of the classical BSDE framework to the fractional case, thus enriching stochastic calculus theory.

• Better understanding of the links between noise regularity and solution properties.

• Development of analysis tools adapted to long-memory processes.

8.2.2. On the Practical Level

• Provision of robust mathematical tools for modeling real phenomena exhibiting long memory.

• Opening of new perspectives for optimal control of complex systems.

• Proposal of efficient numerical methods for concrete implementation.

8.3. Research Perspectives

Several research directions deserve to be explored as extensions of this work.

8.3.1. Theoretical Extensions

• Rough case ($H<1/2$ ): Development of a theory for $H$ values less than 1/2, requiring the use of rough stochastic calculus and advanced approximation techniques.

• Quadratic BSDEs: Study of BSDEs with generator having quadratic growth in $z$ , presenting substantial mathematical challenges but offering extended applications.

• Infinite-dimensional systems: Extension to the framework of BSDEs in Hilbert spaces, to model distributed systems and stochastic partial differential equations.

8.3.2. Numerical Developments

• Higher-order schemes: Design of high convergence order algorithms, exploiting in particular multi-step methods and extrapolation techniques.

• Complexity reduction: Development of model reduction methods and sparse approximation to treat high-dimensional problems.

• Artificial intelligence integration: Combination of traditional methods with deep learning approaches for approximation of conditional expectations.

8.3.3. Advanced Applications

• Mathematical finance: Pricing and hedging of exotic options depending on market memory, particularly in markets presenting persistence anomalies.

• Systems engineering: Optimal control of physical systems exhibiting hysteresis phenomena or long relaxation times.

• Epidemiology and biology: Modeling and control of epidemic or biological processes exhibiting long temporal dependence.

• Energy and environment: Optimization of renewable energy production and distribution, taking into account the persistence of meteorological conditions.

8.4. Scientific and Societal Impact

The advances presented in this work contribute to expanding the field of financial mathematics and stochastic control. On the societal level, the fine modeling of long-memory phenomena opens perspectives for better financial risk management, more efficient optimization of energy systems, and improvement of control strategies in various application domains.

The robustness of the theoretical results coupled with the efficiency of the proposed numerical methods positions this work as a solid foundation for future research, both fundamental and applied, in the domain of fractional BSDEs.

Conflicts of Interest

The author declares no conflicts of interest regarding the publication of this paper.

References

[1]	Pardoux, E. and Peng, S.G. (1990) Adapted Solution of a Backward Stochastic Differential Equation. Systems & Control Letters, 14, 55-61. [Google Scholar] [CrossRef]
[2]	Peng, S. (1991) Probabilistic Interpretation for Systems of Quasilinear Parabolic Partial Differential Equations. Stochastics and Stochastics Reports, 37, 61-74.
[3]	Hu, Y. and Øksendal, B. (1996) Fractional White Noise Calculus and Applications to Finance. Infinite Dimensional Analysis, Quantum Probability and Related Topics, 6, 1-32.
[4]	Biagini, F., Hu, Y., Øksendal, B. and Zhang, T. (2008) Stochastic Calculus for Fractional Brownian Motion and Applications. Springer. [Google Scholar] [CrossRef]
[5]	Nualart, D. (2006) The Malliavin Calculus and Related Topics. 2nd Edition, Springer.
[6]	Mishura, Y. (2008) Stochastic Calculus for Fractional Brownian Motion and Related Processes. Springer. [Google Scholar] [CrossRef]
[7]	Coutin, L. (2007) An Introduction to (stochastic) Calculus with Respect to Fractional Brownian Motion. In: Donati-Martin, C., Émery, M., Rouault, A. and Stricker, C., Eds., Séminaire de Probabilités XL, Springer, 3-65. [Google Scholar] [CrossRef]
[8]	Bender, C. and Zhang, T. (2008) Time-Discretization of Stochastic Differential Equations with Fractional Brownian Motion. Mathematics of Computation, 77, 1587-1610.
[9]	Zhang, X. (2017) Backward Stochastic Differential Equations with Fractional Brownian Motion. Stochastic Processes and their Applications, 127, 4035-4055.
[10]	Nostratipour, M. and Hamdi, M. (2022) Numerical Methods for Fractional Backward Stochastic Differential Equations. Journal of Computational and Applied Mathematics, 405, Article 113972.