Stability analysis of stochastic systems with fractional dynamics

Abstract

This paper investigates the stability properties of stochastic functional differential equations driven by fractional Brownian motion (FBM), a natural extension of classical Brownian motion that incorporates memory effects through the Hurst parameter. By constructing appropriate Lyapunov functionals, we derive sufficient conditions for global exponential mean square stability and asymptotic stochastic stability. The analysis is applied to a delayed stochastic fractional Black–Scholes model, where explicit criteria for both stability types are established. Numerical simulations are presented to validate the theoretical results and illustrate the influence of the Hurst index on the stability behavior of the system.

Introduction

Numerous real-world applications are described by stochastic differential equations (SDEs) and stochastic functional differential equations (SFDEs), which are affected by the classical standard Brownian motion. The traditional principles of stochastic calculus can be applied to any kind of uncertainty or variation. Numerous authors have thoroughly examined the theory of stochastic processes and stochastic calculus; for example, see Refs. 1–4. In 1944, ² published a famous study that provided the well-known formula known as the Ito formula. In classical calculus, this formula is equivalent to the Libniz-Newton chain rule. The dynamics of numerous well-known mathematical models in fields like ecology, epidemiology, and economics can be shown by incorporating Brownian motion into the mathematical model; for instance, Refs. 1–4. Now, we can reveal the dynamics of many known mathematical models in economics, epidemiology, ecology, etc., see Refs. 5–9.

Itô stochastic calculus was developed for semimartingale processes to make sense of the stochastic integrals, but the fractional Brownian motions B ^H (t), H ≠ 1/2 are not semimartingale processes. Consequently, describing the phenomena by the theory of semimartingale is insufficient to describe models involving long memory. Fractional Brownian motion was first introduced by Hurst,^10,11 and it has played a vital role in many applications in economics and communications, see Refs. 12–16.

The author in Ref. 17 has developed stochastic differential equations for fBm-driven systems which extend stochastic calculus to model processes with memory effects, relevant for financial markets and physics. Existence and uniqueness of a class of stochastic differential equations driven by fractional Brownian motions are studied in Ref. 18. Some studies advance stochastic analysis for systems with memory effects driven by fractional Brownian motion. Agram et al. ¹⁹ developed optimal control methods for delayed systems using novel fractional calculus techniques. Xu et al. ²⁰ established averaging principles that simplify complex delayed systems through convergence approximations. Together, they provide powerful tools for analyzing non-Markovian processes in economics, finance, and engineering. Our work makes several key advances in the stability analysis of fractional Brownian motion (FBM)-driven systems because it bridges a critical gap by extending Lyapunov-based stability methods to stochastic functional differential equations with FBM noise, addressing the non-Markovian and non-semimartingale challenges inherent to such systems. The derived sufficient conditions for global exponential mean square stability and asymptotic stochastic stability provide the first explicit criteria tailored to delayed fractional Black–Scholes models, a previously underexplored area. The numerical analysis reveals how the Hurst parameter quantitatively modulates stability thresholds, offering new insights into memory-dependent dynamics. The methodology allows rigorous stability control in financial and engineering applications with memory effects. This work stands out by unifying theoretical advances (Lyapunov functionals for FBM) with actionable criteria for real-world systems, filling a gap between abstract theory and practical stability guarantees.

Define the Hurst index 0 < H < 1, the fractional Brownian motion B ^H (t), t ≥ 0 is a continuous Gaussian process with zero means, continuous trajectories, and a correlation function

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

The fractional Brownian motion has the properties

1. $\mathbb{E}\left[B^H(0)\right]=0$ , and $\mathbb{E}\left[{\left(B^H(t)\right)}^2\right]=t^{2H}$ .

2. B ^H (t) has homogenous increments.

For H = 1/2, the B ^H (t) is a semimartingale process. Many problems described by differential equations driven by fractional Brownian motion have been studied by Refs. 21–23 in various fields.

Delay differential equations consider the past states of the system, and it is an important feature in the population and financial dynamics. In our work, we consider the delay and the uncertainty in the stability analysis. Neglecting small delays in the mathematical models does not lead to the rightful conclusions. Also, neglecting the uncertainty does not capture the right behavior of the solution. Many stochastic delayed mathematical models have been studied by Refs. 24–27.

Define J≔C([−τ, 0], L₂), a Banach space of mean square continuous functionals φ defined on [−τ, 0] with the norm

$\Vert \phi {\Vert }_J=\underset{-\tau \le s\le 0}{\sup }\Vert \phi {\Vert }_2=\underset{-\tau \le s\le 0}{\sup }{\left(\mathbb{E}\left[{\phi }^2(s)\right]\right)}^{1/2}.$

Consider the following nonlinear delay differential equation driven by fractional Brownian motion

$dx(t)=f(t,x_t)dt+g(t,x_t)dB(t),x\in L_2(\mathrm{\Omega }),f,g:{\mathbb{R}}_{+}\times J\to \mathbb{R},$

(1)

with the initial history function

$x(t_0)=\varphi ,\varphi \in J,$

(2)

where x _t = x(t + s), s ≤ 0. With the Hurst parameter $H>\frac{1}{2}$ , the standard fractional Brownian motion $B^H(t):[t_0,\infty )\to \mathbb{R}$ is defined on the complete probability space $(\mathrm{\Omega },\mathcal{F},{\mathbb{P}}^H)$ . The continuous functionals f, g are supposed to satisfy the local mean square Lipschitz condition, that is, for C₁, C₂ > 0

$\Vert f(t,x_t)-f(t,x_t^{\prime }){\Vert }_2\le C_1\Vert x_t-x_t^{\prime }{\Vert }_2,$

$\Vert g(t,x_t)-g(t,x_t^{\prime }){\Vert }_2\le C_2\Vert x_t-x_t^{\prime }{\Vert }_2.$

We assume that the two functionals f, g satisfy

$f(t,0)=0,g(t,0)=0.$

Consequently, system (1) with (2) admits the zero solution. Moreover, $f,g:[t_0,\infty )\times \mathcal{C}[-\tau ,0]\to \mathbb{R}$ are Borel measurable and belong to the Hilbert space J _ψ , where $\psi (t,s)=H(2H-1){\left|s-t\right|}^{2H-2}$ , then

$\begin{array}{l}\hfill {\left|f(t,x_t)\right|}_{\psi }^2={\int }_{t_0}^{\infty }{\int }_{t_0}^{\infty }f(s,x_s)f(t,x_t)\psi (t,s)dsdt<\infty ,\\ \hfill {\left|g(t,x_t)\right|}_{\psi }^2={\int }_{t_0}^{\infty }{\int }_{t_0}^{\infty }g(s,x_s)g(t,x_t)\psi (t,s)dsdt<\infty .\end{array}$

Now, if g ∈ J _ψ , then the stochastic integral ${\int }_{t_0}^{\infty }g(s,x_s)d{B}^H(s)$ is well-defined with zero mean, and this can be accomplished using the Riemann sums using the Wick product as ${\sum }_{i=1}^{n-1}{\mathcal{X}}_{t_i}^{\pi }\diamond \left(B^H(t_{i+1}^n)-B^H(t_i^n)\right)$ , where the Wick product of two random variables is denoted by ⋄, and $\mathcal{X}$ is a second-order stochastic process.

For the Hurst index $H>\frac{1}{2}$ , we have the fractional Itô formula that can be applied to the Lyapunov functional $V(t,x_t)\in C^{\mathrm{2,1}}([0,\infty )\times J,{\mathbb{R}}_{+})$ , and the stochastic differential of V will be in the form

$\begin{array}{cc}\hfill & dV(t,x_t)\hfill \\ \hfill & =\underset{L^{H}V(t,x_t)}{\underbrace{\left[\frac{\partial V(t,x_t)}{\partial t}+\frac{\partial V(t,x_t)}{\partial x}f(t,x_t)+H(2H-1)\frac{{\partial }^{2}V(t,x_t)}{\partial x^2}g(t,x_t){\int }_{t_0}^{t}g(s,x_s){\left|t-s\right|}^{2H-2}ds\right]}}dt\hfill \\ \hfill & +\frac{\partial V(t,x_t)}{\partial x}g(t,x_t)d{B}^H(t).\hfill \end{array}$

The main contribution of this manuscript is the stability analysis of stochastic delay systems influenced by stochastic fractional Brownian motion with Hurst parameter H > 1/2. The paper addressed the global mean square exponential stability which makes any trajectory tends to the attractor of the system at an exponential rate regardless of the initial history function. Consequently, it is very important to address this type of stability whereas few researchers do. Moreover, as an example, in finance, the Black–Scholes market model is very effective in describing stock prices, and we shall study the behavior of the solution of the stochastic fractional version of this model. It should be noted that the effect of the Hurst index on the stability of this model. We introduce the stochastic fractional Black–Scholes model with delay and without delay with some computer simulations.

This paper is organized as follows. Sect. 2 is devoted to the main results, and our two main theorems of stability analysis of the zero solution of the general equations (1) and (2) are introduced. As an application, in Sect. 3, the stochastic fractional Black–Scholes market model is studied with delay τ and without delay.

Stability analysis

Stochastic stability, mean square stability, and global exponential mean square stability have received very much attraction. By introducing suitable Lyapunov functionals, we can investigate the necessary criteria of asymptotic stability in probability and global mean square exponential stability of the trivial equilibrium of the fractional stochastic functional differential system in the form (1, 2). Introducing different choices of Lyapunov functionals implies more stability conditions, for more details with regard to Lyapunov functionals and their role in the study of the stability of deterministic and stochastic systems, see Refs. 28–32. We mainly focus on stochastic stability and global mean square exponential stability of the general nonlinear delay differential equation driven by the fractional Brownian motion.

Theorem 2.1

The zero solution of (1) is stochastically asymptotically stable.

Proof

Let the stochastic process {x(t), t ≥ 0} is ${\mathcal{F}}_t$ -measurable, then for s < t, x(t) is super-martingale if

$\mathbb{E}\left[x(t)|{\mathcal{F}}_s\right]\le x(s).$

For ‖ϕ‖ < δ, δ > 0, assume that

$V(t,x_t)\ge c_1|x(t){|}^2,c_1>0$

The positive definite Lyapunov functional V(t, x _t ) is super-martingale as L ^H V(t, x _t ) ≤ 0.

From the boundedness of the Lyapunov functional at t = t₀, we have

$V(t_0,\varphi )\le c_2\Vert \varphi {\Vert }^2,c_2>0.$

Based on the facts above, we have

$\begin{array}{c}\hfill \mathbb{P}\left[\underset{t}{\sup }|x(t;t_0,\varphi )|>{\epsilon }_1|{\mathcal{F}}_{t_0}\right]\le \mathbb{P}\left[\underset{t}{\sup }\frac{\sqrt{\mathcal{V}(t,x(t))}}{c_1}>{\epsilon }_1|{\mathcal{F}}_{t_0}\right]\\ =\mathbb{P}\left[\underset{t}{\sup }\mathcal{V}(t,x(t))>c_1{\epsilon }_1^2|{\mathcal{F}}_{t_0}\right]\hfill \\ \le \frac{\mathcal{V}(t_0\varphi )}{c_1{\epsilon }_1^2}\le \frac{c_2\Vert \varphi {\Vert }^2}{c_1{\epsilon }_1^2}\le \frac{c_2{\delta }^2}{c_1{\epsilon }_1^2}<{\epsilon }_2.\hfill \end{array}$

This completes the proof of stochastic stability of the zero solution of (1). Assume that the positive definite Lyapunov functional $V(t,x_t)\in {\mathcal{C}}^{\mathrm{2,1}}\left([t_0,\infty )\times J,{\mathbb{R}}_{+}\right)$ is continuous and bounded at t = t₀, then

$\underset{\Vert x_t\Vert \in J_{\delta }}{\sup }V(t_0,\varphi )\le {\epsilon }_1\nu (r),\epsilon \in (\mathrm{0,1}),r>0,$

i.e., δ < r. Define the stopping time

$\kappa =\inf \{t\ge t_0;\Vert x(t)\Vert \notin J_r\}.$

Then by applying the fractional Itô formula to V(t, x _t ), we get

$\begin{array}{l}\hfill d\mathcal{V}(\kappa \wedge t,x_{\kappa \wedge t})={\mathrm{L}}^H\mathcal{V}(t,x_t)dt+g(t,x_t)\frac{\partial \mathcal{V}(t,x_t)}{\partial x}d{B}^H(t),\\ \hfill \mathcal{V}(\kappa \wedge t,x_{\kappa \wedge t})-\mathcal{V}(t_0,\varphi )={\int }_{t_0}^{\kappa \wedge t}{\mathrm{L}}^H\mathcal{V}(s,x_s)ds+{\int }_{t_0}^{\kappa \wedge t}\sigma (s,x_s)\frac{\partial \mathcal{V}(s,x_s)}{\partial x}d{B}^H(s).\end{array}$

$\mathbb{E}\mathcal{V}(\kappa \wedge t,x_{\kappa \wedge t})=\mathbb{E}\mathcal{V}(t_0,\varphi )+\mathbb{E}{\int }_{t_0}^{\kappa \wedge t}{\mathrm{L}}^H\mathcal{V}(s,x_s)ds.$

As $\mathcal{V}(t,x_t)$ is a suitable Lyapunov functional with ${\mathrm{L}}^H\mathcal{V}(t,x_t)\le 0$ , then

$\mathbb{E}\mathcal{V}(\kappa \wedge t,x_{\kappa \wedge t})\le \mathcal{V}(t_0,\varphi ),$

(3)

and if κ < t, then

$\begin{array}{c}\hfill \epsilon \upsilon (r)\le \mathbb{E}\mathcal{V}(\kappa ,x_{\kappa })\le \mathcal{V}(t_0,\varphi ),\\ \hfill \mathbb{P}(\kappa \le t)\upsilon (r)\le \mathbb{E}\mathcal{V}(\kappa ,x_{\kappa })\le \mathcal{V}(t_0,\varphi ),\end{array}$

we obtain $\mathbb{P}\left(\kappa \le t\right)\le \epsilon$ , then $\mathbb{P}(\kappa <\infty )\le \epsilon$ as t → ∞, then for t ≥ t₀

$\mathbb{P}(|x(t;t_0,\varphi )|<r)=1-\mathbb{P}(\kappa <\infty )\ge 1-\epsilon .$

(4)

By choosing δ ∈ [0, ℓ/2], for ℓ > 0 and 0 < ɛ < 1, then (4) implies

$\mathbb{P}\left(|x(t;t_0,\varphi )|<\frac{\ell }{2}\right)\ge 1-\frac{\epsilon }{2}.$

(5)

Define κ _α = inf(t ≥ 0; ‖x(t)‖ ≤ α) and ${\kappa }_{\ell }=\inf (t\ge 0;\Vert x(t)\Vert \ge \frac{\ell }{2})$ . By fractional Itô formula, we have

$\begin{array}{cc}\hfill d\mathcal{V}({\kappa }_{\alpha }\wedge t\wedge {\kappa }_{\ell },x_{{\kappa }_{\alpha }\wedge t\wedge {\kappa }_{\ell }})& ={\mathrm{L}}^H\mathcal{V}(t,x_t)dt+g(t,x_t)\frac{\partial \mathcal{V}(t,x_t)}{\partial x}d{B}^H(t)\\ \hfill \mathcal{V}({\kappa }_{\alpha }\wedge t\wedge {\kappa }_{\ell },x_{{\kappa }_{\alpha }\wedge t\wedge {\kappa }_{\ell }})-\mathcal{V}(t_0,\varphi )& ={\int }_{t_0}^{{\kappa }_{\alpha }\wedge t\wedge {\kappa }_{\ell }}{\mathrm{L}}^H\mathcal{V}(s,x_s)ds+{\int }_{t_0}^{{\kappa }_{\alpha }\wedge t\wedge {\kappa }_{\ell }}\sigma (s,x_s)\frac{\partial \mathcal{V}(s,x_s)}{\partial x}d{B}^H(s).\end{array}$

$\mathbb{E}\mathcal{V}({\kappa }_{\alpha }\wedge t\wedge {\kappa }_{\ell },x_{{\kappa }_{\alpha }\wedge t\wedge {\kappa }_{\ell }})=\mathcal{V}(t_0,\varphi )+\mathbb{E}{\int }_{t_0}^{{\kappa }_{\alpha }\wedge t\wedge {\kappa }_{\ell }}{\mathrm{L}}^H\mathcal{V}(s,x_s)ds.$

Assume that $\mathcal{V}(t,x_t)$ has infinitesimal upper bound with ${\mathrm{L}}^H\mathcal{V}(t,x_t)\le -\upsilon (\Vert x(t)\Vert )$ , then

$\begin{array}{cc}\hfill 0\le \mathbb{E}\mathcal{V}({\kappa }_{\alpha }\wedge t\wedge {\kappa }_{\ell },x_{{\kappa }_{\alpha }\wedge t\wedge {\kappa }_{\ell }})& \le \mathcal{V}(t_0,\varphi )-\mathbb{E}\left[\upsilon (\alpha ){\int }_{t_0}^{{\kappa }_{\alpha }\wedge t\wedge {\kappa }_{\ell }}ds\right]\hfill \\ \hfill & =\mathcal{V}(t_0,\varphi )-\upsilon (\alpha )\mathbb{E}\left({\kappa }_{\alpha }\wedge t\wedge {\kappa }_{\ell }-t_0\right),\hfill \end{array}$

(6)

implies $\upsilon (\alpha )\mathbb{E}\left({\kappa }_{\alpha }\wedge t\wedge {\kappa }_{\ell }-t_0\right)\le \mathcal{V}(t_0,\varphi )$ . From (5), $\mathbb{P}({\kappa }_{\ell }<\infty )\le \frac{\epsilon }{2}$ and as t → ∞

$\mathbb{P}({\kappa }_{\alpha }\wedge {\kappa }_{\ell }<\infty )=1\le \mathbb{P}({\kappa }_{\alpha }<\infty )+\mathbb{P}({\kappa }_{\ell }<\infty )\le \mathbb{P}({\kappa }_{\alpha }<\infty )+\frac{\epsilon }{2}.$

Now, let $0<\alpha <\beta <\frac{\ell }{2}$ , (β is arbitrary), then $\mathbb{P}(|x(t;t_0,\varphi )|<\beta )\ge 1-\frac{\epsilon }{2}$ and

$\begin{array}{cc}\hfill \mathbb{P}(\underset{t\to \infty }{\lim }|x(t,\varphi )|\le \beta )& \ge \mathbb{P}\left(\{{\kappa }_{\alpha }<\infty \}\cap \{|x(t,\varphi )|\le \beta \}\right)\hfill \\ \hfill & =\mathbb{P}({\kappa }_{\alpha }<\infty )\mathbb{P}(|x(t,\varphi )|\le \beta \mid {\kappa }_{\alpha }<\infty )\hfill \\ \hfill & \ge \left(1-\frac{\epsilon }{2}\right)\left(1-\frac{\epsilon }{2}\right)\ge 1-\epsilon .\hfill \end{array}$

β is arbitrary, so $\mathbb{P}({\lim }_{t\to \infty }|x(t;t_0,\varphi )|=0)\ge 1-\epsilon$ . Then the zero solution is stochastically asymptotically stable.

Theorem 2.2

The zero solution of (1) is globally exponentially mean square stable.

Proof

Assume that there exists a Lyapunov functional V(t, x _t ) such that

1. c₁‖x(t)‖² ≤ V(t, x _t ) ≤ c₂‖x(t)‖², c₁, c₂ > 0

2. ${\mathrm{L}}^H\mathcal{V}(t,x_t)\le -c_{3}V(t,x_t),c_3>0$

Applying the fractional Itô formula to the functional

$V(t,x_t)=e^{c_3(t-t_0)}V(t,x_t),$

which is sufficiently smooth, then for H > 1/2

$\begin{array}{cc}\hfill e^{c_3(t-t_0)}V(t,x_t)& =V(t_0,\varphi )+{\int }_{t_0}^t{e}^{c_3(s-t_0)}V(s,x_s)ds+{\int }_{t_0}^t{e}^{c_3(s-t_0)}\frac{\partial c_{3}V(s,x_s)}{\partial s}ds\hfill \\ \hfill & +{\int }_{t_0}^t{e}^{c_3(s-t_0)}\frac{\partial V(s,x_s)}{\partial x}f(s,x_s)ds+{\int }_{t_0}^t{e}^{c_3(s-t_0)}\frac{\partial V(s,x_s)}{\partial x}g(s,x_s)d{B}^H(s)\hfill \\ \hfill & +H(2H-1){\int }_{t_0}^t{e}^{c_3(s-t_0)}\frac{{\partial }^{2}V(s,x_s)}{\partial x^2}g(s,x_s)\left[{\int }_{t_0}^{s}g(r,x_r){\left|s-r\right|}^{2H-2}dr\right]ds.\hfill \end{array}$

Define the stopping time

${\kappa }_h=\inf \{t\ge t_0;\Vert x(t)\Vert \ge h\},$

where κ _h → ∞, as h → ∞. Then

$\begin{array}{cc}\hfill e^{c_3({\kappa }_h\wedge t-t_0)}& V({\kappa }_h\wedge t,x_{{\kappa }_h\wedge t})\hfill \\ \hfill & =V(t_0,\varphi )+{\int }_{t_0}^{{\kappa }_h\wedge t}{e}^{c_3(s-t_0)}{c}_{3}V(s,x_s)ds+{\int }_{t_0}^{{\kappa }_h\wedge t}{e}^{c_3(s-t_0)}\frac{\partial V(s,x_s)}{\partial s}ds\hfill \\ \hfill & +{\int }_{t_0}^{{\kappa }_h\wedge t}{e}^{c_3(s-t_0)}\frac{\partial V(s,x_s)}{\partial x}f(s,x_s)ds+{\int }_{t_0}^{{\kappa }_h\wedge t}{e}^{c_3(s-t_0)}\frac{\partial V(s,x_s)}{\partial x}g(s,x_s)d{B}^H(s)\hfill \\ \hfill & +H(2H-1){\int }_{t_0}^{{\kappa }_h\wedge t}{e}^{c_3(s-t_0)}\frac{{\partial }^{2}V(s,x_s)}{\partial x^2}g(s,x_s)\left[{\int }_{t_0}^{{\kappa }_h\wedge s}g(r,x_r){\left|s-r\right|}^{2H-2}dr\right]ds.\hfill \end{array}$

Taking the expectation implies

$\begin{array}{cc}\hfill e^{c_3({\kappa }_h\wedge t-t_0)}& \mathbb{E}\left[V({\kappa }_h\wedge t,x_{{\kappa }_h\wedge t})\right]\hfill \\ \hfill & =\mathbb{E}\left[V(t_0,\varphi )\right]+\mathbb{E}\left[{\int }_{t_0}^{{\kappa }_h\wedge t}{e}^{c_3(s-t_0)}\left(c_{3}V(s,x_s)+\frac{\partial V(s,x_s)}{\partial s}+\frac{\partial V(s,x_s)}{\partial x}f(s,x_s)\right.\right.\hfill \\ \hfill & \left.\left.+H(2H-1)\frac{{\partial }^{2}V(s,x_s)}{\partial x^2}g(s,x_s){\int }_{t_0}^{{\kappa }_h\wedge s}g(r,x_r){\left|s-r\right|}^{2H-2}dr\right)\right]dt\hfill \\ \hfill & =\mathbb{E}\left[V(t_0,\varphi )\right]+\mathbb{E}\left[{\int }_{t_0}^{{\kappa }_h\wedge t}{e}^{c_3(s-t_0)}\left(c_{3}V(s,x_s)+{\mathrm{L}}^{H}V(s,x_s)\right)ds\right]\hfill \end{array}$

From the assumption 2, we have

$e^{c_3({\kappa }_h\wedge t-t_0)}\mathbb{E}\left[V({\kappa }_h\wedge t,x_{{\kappa }_h\wedge t})\right]\le V(t_0,\varphi ).$

Then assumption 1 implies

$c_1{e}^{c_3(t-t_0)}\Vert x(t){\Vert }^2\le e^{c_3({\kappa }_h\wedge t-t_0)}V({\kappa }_h\wedge t,x_{{\kappa }_h\wedge t})\le V(t_0,\varphi )\le c_2\Vert \varphi {\Vert }^2.$

That is, as n → ∞

$c_1{e}^{c_3(t-t_0)}\Vert x(t){\Vert }^2\le c_2\Vert \varphi {\Vert }^2.$

Hence,

$\mathbb{E}\left[|x(t;t_0,\varphi ){|}^2\right]\le \frac{c_2}{c_1}{e}^{-c_3(t-t_0)}\Vert \varphi {\Vert }^2,$

which proves the exponential mean square stability.

Now, choose a proper constant 0 ≤ ℓ < 1 such that

$|\ell |e^{-\epsilon \tau }<1.$

(7)

The Lyapunov functional

$\mathcal{V}(t,x_t)=(1+|\ell |)e^{\epsilon t}{\left[x(t)\right]}^2+{\int }_{t-\tau }^t{e}^{\epsilon s}{x}^2(s)ds,$

implies

$\mathcal{V}(t_0,\varphi )\le (1+|\ell |)e^{\epsilon t_0}\Vert \varphi {\Vert }^2+\tau e^{\epsilon t_0}\Vert \varphi {\Vert }^2.$

From fractional Itô formula

$d\mathcal{V}(t,x_t)={\mathrm{L}}^H\mathcal{V}(t,x_t)dt+{\mathcal{V}}_x(t,x_t)g(t,x_t)d{B}^H(t),$

then

$\begin{array}{c}\hfill \mathbb{E}\left[\mathcal{V}(t,x_t\right]-\mathbb{E}\left[\mathcal{V}(t_0,\varphi )\right]=\mathbb{E}{\int }_{t-\tau }^t{\mathrm{L}}^H\mathcal{V}(s,x_s)ds+\mathbb{E}\left[{\int }_{t-\tau }^t{\mathcal{V}}_x(s,x_s)\sigma (s,x_s)d{B}^H(s)\right].\end{array}$

Taking the expectation implies

$\mathbb{E}\left[\mathcal{V}(t,x_t\right]-\mathbb{E}\left[\mathcal{V}(t_0,\varphi )\right]\le 0.$

Then

$(1+|\ell |)e^{\epsilon t}\mathbb{E}{\left[x(t)\right]}^2\le \mathbb{E}\left[\mathcal{V}(t,x_t)\right]\le \mathbb{E}\left[\mathcal{V}(t_0,\varphi )\right]\le (1+|\ell |)e^{\epsilon t_0}\Vert \varphi {\Vert }^2+\tau e^{\epsilon t_0}\Vert \varphi {\Vert }^2.$

Let M = (1 + |ℓ| + τ) ≥ 1, then

$\mathbb{E}\left|x^2(t)+|\ell |x^2(t)\right|\le M\Vert \varphi {\Vert }^2{e}^{-\epsilon (t-t_0)}.$

Consequently,

$\begin{array}{cc}\hfill \mathbb{E}|x(t){|}^2& \le |\alpha |\mathbb{E}{\left[x(t)\right]}^2+\mathbb{E}\left|x^2(t)+|\ell |x^2(t)\right|\hfill \\ \hfill & \le |\ell |\mathbb{E}{\left[x(t-\tau )\right]}^2+M\Vert \varphi {\Vert }^2{e}^{-\epsilon (t-t_0)}.\hfill \end{array}$

Then by the induction and (7), we have to show for global stability.

$\mathbb{E}|x(t){|}^2\le \frac{M\Vert \varphi {\Vert }^2}{1-|\ell |e^{\epsilon \tau }}{e}^{-\epsilon (t-t_0)},t\ge t_0,$

(8)

First, for $\left[t_0,t_0+\tau \right)$ ,

$\begin{array}{cc}\hfill \mathbb{E}|x(t){|}^2& \le |\ell |\mathbb{E}{\left[x(t-\tau )\right]}^2+M\Vert \varphi {\Vert }^2{e}^{-\epsilon (t-t_0)}\hfill \\ \hfill & \le \Vert \varphi {\Vert }^2\left[|\ell |+M{e}^{-\epsilon (t-t_0)}\right]\hfill \\ \hfill & \le M\Vert \varphi {\Vert }^2\left[|\ell |+e^{-\epsilon (t-t_0)}\right]\hfill \\ \hfill & \le M\Vert \varphi {\Vert }^2{e}^{-\epsilon (t-t_0)}\left[1+|\ell |e^{\epsilon \tau }\right]\hfill \\ \hfill & \le \frac{M\Vert \varphi {\Vert }^2}{1-|\ell |e^{\epsilon \tau }}{e}^{-\epsilon (t-t_0)}.\hfill \end{array}$

(9)

For $\left[t_0+\tau ,t_0+2\tau \right)$ and using (9), we deduce

$\begin{array}{cc}\hfill \mathbb{E}|x(t){|}^2& \le |\alpha |\mathbb{E}{\left[x(t-\tau )\right]}^2+M\Vert \varphi {\Vert }^2{e}^{-\epsilon (t-t_0)}\hfill \\ \hfill & \le M\Vert \varphi {\Vert }^2{e}^{-\epsilon (t-t_0)}+|\ell |\left[M\Vert \varphi {\Vert }^2{e}^{-\epsilon (t-\tau -t_0)}\left(1+|\ell |e^{\epsilon \tau }\right)\right]\hfill \\ \hfill & =M\Vert \varphi {\Vert }^2{e}^{-\epsilon (t-t_0)}+M\Vert \varphi {\Vert }^2{e}^{-\epsilon (t-t_0)}\left[|\ell |e^{\epsilon \tau }+|\ell {|}^2{e}^{2\epsilon \tau }\right]\hfill \\ \hfill & =M\Vert \varphi {\Vert }^2{e}^{-\epsilon (t-t_0)}\left[1+|\ell |e^{\epsilon \tau }+|\ell {|}^2{e}^{2\epsilon \tau }\right]\hfill \\ \hfill & \le \frac{M\Vert \varphi {\Vert }^2}{1-|\ell |e^{\epsilon \tau }}{e}^{-\epsilon (t-t_0)}.\hfill \end{array}$

Similarly, for $\left[t_0+n\tau ,t_0+(n+1)\tau \right),n=\mathrm{1,2},\cdots$ , we get

$\begin{array}{cc}\hfill \mathbb{E}|x(t){|}^2& \le |\ell |\mathbb{E}{\left[x(t-\tau )\right]}^2+M\Vert \varphi {\Vert }^2{e}^{-\epsilon (t-t_0)}\hfill \\ \hfill & \le M\Vert \varphi {\Vert }^2{e}^{-\epsilon (t-t_0)}\left[1+|\ell |e^{\epsilon \tau }+|\ell {|}^2{e}^{2\epsilon \tau }+\cdots +|\alpha {|}^{n+1}{e}^{(n+1)\epsilon \tau }\right]\hfill \\ \hfill & \le \frac{M\Vert \varphi {\Vert }^2}{1-|\ell |e^{\epsilon \tau }}{e}^{-\epsilon (t-t_0)}.\hfill \end{array}$

Hence, (8) is satisfied. Consequently, the zero solution of (1) is globally exponentially stable in the sense of the mean square.

Application: Fractional Black–Scholes market model

The Black–Scholes market model in its classical form established by Ref. 33 plays a vital role in describing the financial pricing option. In this part, we study this model driven by fractional Brownian motion with discrete delay τ > 0. Stochastic stability, asymptotic stochastic stability, and global mean square exponential stability of this model are investigated. Moreover, we consider the fractional Black–Scholes model without delay and investigate the effect of the Hurst parameter on the stability. Some numerical simulations are carried out. Consider the fractional delayed Black–Scholes model in the form

$dx(t)=\alpha x(t)dt+\beta x(t-\tau )d{B}^H(t).$

(10)

The parameters α, β are arbitrary. B ^H (t) is the fractional Brownian motion with the Hurst index H ∈ (0, 1).

Proposition 3.1

Equation (10) is

1. Stochastically stable if L ^H V(t, x _t ) ≤ 0 and α ≤ λτ/2, λ > 0.

2. Stochastically asymptotically stable if ${\mathrm{L}}^{H}V(t,x_t)\le -\nu \left(\Vert x(t)|\right)$ and α ≤ 1/2(λτ + ɛ), ɛ is arbitrary positive and $\nu \in \mathcal{K}$ .

3. Globally exponentially mean square stable if

(a) c₁‖x(t)‖² ≤ V(t, x _t ) ≤ c₂‖x(t)‖², c₁, c₂ > 0

(b) ${\mathrm{L}}^H\mathcal{V}(t,x_t)\le -c_{3}V(t,x_t),c_3>0$

Proof

By introducing the stochastic Lyapunov functional

$V(t,x_t)=e^{-\lambda \tau t}{x}^2(t),\lambda >0.$

(11)

This is a positive definite functional as

$M_1\le e^{-\lambda \tau t}\le M_2,$

for positive constants M₁, M₂ and V(0, t) = 0. Then

$d\mathcal{V}(t,x_t)={\mathrm{L}}^H\mathcal{V}(t,x_t)dt+\beta x(t-\tau ){\mathcal{V}}_x(t,x_t)d{B}^H(t),$

with

$\begin{array}{cc}\hfill {\mathrm{L}}^H\mathcal{V}(t,x_t)& =-\lambda \tau e^{-\lambda \tau t}{x}^2(t)+2\alpha e^{-\lambda \tau t}{x}^2(t)+2H(2H-1){\beta }^2{x}^2(t-\tau )e^{-\lambda \tau t}{x}^2(t){\int }_{t_0}^t|t-r{|}^{2H-2}dr\hfill \\ \hfill & =\left(2\alpha -\lambda \tau \right)e^{-\lambda \tau t}{x}^2(t)+2{\beta }^{2}H{x}^2(t-\tau )t^{2H-1}{e}^{-\lambda \tau t}.\hfill \end{array}$

(12)

The last term vanishes for sufficiently large values of λ, then ${\mathrm{L}}^H\mathcal{V}(t,x_t)\le 0$ if 2α − λτ ≤ 0.

According to Theorem (2.1), equation (10) is stochastically stable.

For the asymptotic stochastic stability, it is sufficient to prove the decreasing property of the Lyapunov functional. It is known $V(t,x_t)\in \mathcal{K}$ such that

$V(t,x_t)\le \nu \left(\Vert x(t)\Vert \right)=(\lambda \tau -2\alpha )e^{-\lambda \tau t}{x}^2(t)-2{\beta }^{2}H{x}^2(t-\tau )t^{2H-1}{e}^{-\lambda \tau t}.$

From (12),

$\begin{array}{cc}\hfill {\mathrm{L}}^H\mathcal{V}(t,x_t)& \le -\nu \left(\Vert x(t)\Vert \right)\hfill \\ \hfill & =\left(2\alpha -\lambda \tau +\epsilon \right)e^{-\lambda \tau t}{x}^2(t)+2{\beta }^{2}H{x}^2(t-\tau )t^{2H-1}{e}^{-\lambda \tau t},\hfill \end{array}$

if only 2α ≤ λτ + ɛ. So equation (10) is stochastically asymptotically stable.

Regarding the exponential stability, consider the same Lyapunov functional (11). Following the same argument of Theorem (2.2), condition (a) is satisfied and from fractional Itô formula, we get

$\begin{array}{cc}\hfill {\mathrm{L}}^H\mathcal{V}(t,x_t)& \le (2\alpha -\lambda \tau )e^{-\lambda \tau t}{x}^2(t)\hfill \\ \hfill & =(2\alpha -\lambda \tau )V(t,x_t),\hfill \end{array}$

which implies condition (b). Therefore, equation (10) is mean square exponential stable for 2α − λτ + ɛ ≤ 0. We can prove the global stability of (10) by following the same argument of Theorem (2.2).

Impact of the Hurst index H

Consider the fractional Black–Scholes market model without delay in the form

$dx(t)=\left(\alpha dt+\beta d{B}^H(t)\right)x(t)$

(13)

Proposition 3.2

With probability 1, equation (13) is stable if

1. α < 0 and H ∈ (0, 1/2)

2. α < β²/2 and H = 1/2

3. α, β arbitrary and H ∈ (1/2, 1)

Moreover, equation (13) is mean square stable if

1. α < 0 and H ∈ (0, 1/2)

2. α < − β²/2 and H = 1/2

Proof

Assume the stochastic process

$y(t)=\ln (x(t)).$

By the fractional Itô formula and (13), we have

$dy(t)=\alpha dt+\beta d{B}^H(t)-{\beta }^{2}H{t}^{2H-1}.$

Then

$y(t)=y(t_0)+\alpha t+\beta b^H(t)-\frac{1}{2}{\beta }^2{t}^{2H},$

That is,

$y(t)=\ln (x_0)+\alpha t+\beta b^H(t)-\frac{1}{2}{\beta }^2{t}^{2H}.$

Consequently,

$x(t)=e^{y(t)}=x_0{e}^{\alpha t+\beta b^H(t)-\frac{1}{2}{\beta }^2{t}^{2H}}.$

The negativeness of the large Lyapunov exponent,

$\underset{t\to \infty }{\lim \sup }\frac{1}{t}\ln |x(t;t_0,x_0)|$

Is a necessary and sufficient condition for stability with probability 1, then

$\begin{array}{cc}\hfill \underset{t\to \infty }{\lim \sup }\frac{1}{t}\ln |x(t;t_0,x_0)|& =\underset{t\to \infty }{\lim \sup }\frac{1}{t}\ln |x_0{e}^{\alpha t+\beta b^H(t)-\frac{1}{2}{\beta }^2{t}^{2H}}|\hfill \\ \hfill & =\underset{t\to \infty }{\lim }\alpha -\frac{1}{2}{\beta }^2{t}^{2H-1}.\hfill \end{array}$

Therefore, for 0 < H < 1/2, the Lyapunov exponent equals α which should be negative for stability, for H = 1/2, the Lyapunov exponent equals α − 1/2β² and it is less than zero for stability, and finally for 1/2 < H < 1, the solution will be stable for all values of α, β

For mean square stability, let

$y(t)=\ln (x^2(t)).$

Then

$\begin{array}{ll}\hfill dy(t)& =2\alpha dt+2\beta d{B}^H(t)-2{\beta }^{2}H{t}^{2H-1},\hfill \\ \hfill y(t)& y(t_0)+2\alpha t+2\beta B^H(t)-{\beta }^2{t}^{2H}=\ln (x^2(t_0)))+2\alpha t+2\beta B^H(t)-{\beta }^2{t}^{2H},\hfill \\ \hfill x^2(t)& ey(t)=x_0^2{e}^{2\alpha t+2\beta B^H(t)-{\beta }^2{t}^{2H}}.\hfill \end{array}$

Then

$\begin{array}{cc}\hfill \underset{t\to \infty }{\lim \sup }\frac{1}{t}\ln |x(t;t_0,x_0){|}^2& =\underset{t\to \infty }{\lim \sup }\frac{1}{t}\mathbb{E}\left[x_0^2{e}^{2\alpha t+2\beta B^H(t)-{\beta }^2{t}^{2H}}\right]\hfill \\ \hfill & =\underset{t\to \infty }{\lim }{x}_0^2\mathbb{E}\left[e^{2\alpha t+2\beta B^H(t)-{\beta }^2{t}^{2H}}\right].\hfill \end{array}$

Now, applying the fractional Itô formula to $e^{2\beta B^H(t)}$ , it is known that

$f(B^H(t))=f(0)+{\int }_0^t{f}^{\prime }(B^H(s))d{B}^H(s)+H{\int }_0^t{f}^{″}(B^H(s))s^{2H-1}ds.$

Then

$e^{2\beta B^H(t)}=1+2\beta {\int }_0^t{e}^{2\beta B^H(t)}d{B}^H(s)+4{\beta }^{2}H{\int }_0^t{e}^{2\beta B^H(t)}{s}^{2H-1}ds.$

Taking the expectation implies

$\mathbb{E}\left[e^{2\beta B^H(t)}\right]=1+4{\beta }^{2}H\mathbb{E}\left[{\int }_0^t{e}^{2\beta B^H(t)}{s}^{2H-1}ds\right].$

We have $e^{2\beta B^H(0)}=1$ , then the solution of the integral equation

$\mathbb{E}\left[e^{2\beta B^H(t)}\right]=2{\beta }^2{t}^{2H}.$

Then

$\begin{array}{cc}\hfill \underset{t\to \infty }{\lim \sup }\frac{1}{t}\ln |x(t;t_0,x_0){|}^2& =\underset{t\to \infty }{\lim }\frac{1}{t}\ln |x_0^2{e}^{2\alpha t+2\beta B^H(t)-{\beta }^2{t}^{2H}}|\hfill \\ \hfill & =\underset{t\to \infty }{\lim }\left(2\alpha +{\beta }^2{t}^{2H-1}\right).\hfill \end{array}$

Consequently,

$\begin{array}{cc}\hfill \underset{t\to \infty }{\lim \sup }\frac{1}{t}\ln |x(t;t_0,x_0){|}^2& =2\alpha <0,H\in \left(0,\frac{1}{2}\right)\text{stable}\hfill \\ \hfill & =2\alpha +{\beta }^2<0,H=\frac{1}{2}\text{stable}\hfill \\ \hfill & =\infty ,H\in \left(\frac{1}{2},1\right)\text{unstable}\hfill \end{array}$

The previous proposition indicates the impact of the Hurst parameter on the stability of (13), and we have some cases

Case 1

1. α = −1, β = 0.05 and H = 0.4, (α < 0, H < 1/2), stable solution, Figure 1(a).

2. α = 1.5, β = 0.1 and H = 0.4, (α > 0), unstable solution, Figure 1(b).

Figure 1.

Numerical simulation of the zero solution for case 1.

Case 2

1. α = 0.004, β = 0.2 and H = 1/2, (α < β²/2), stable solution, Figure 2(a).

2. α = 0.1, β = 0.02 and H = 1/2, (α > β²/2), unstable solution, Figure 2(b).

Figure 2.

Numerical simulation of the zero solution for case 2.

Case 3

1. α = −1, β = 2 and H = 0.8, stable solution, Figure 3(a).

2. α = 1.5, β = −0.5 and H = 0.6, stable solution, Figure 3(b).

3. α = −0.7, β = −0.9 and H = 0.7, stable solution, Figure 3(c).

4. α = 0.02, β = 0.3 and H = 0.95, stable solution, Figure 3(d)

5. α = 0.02, β = 0.3 and H = 0.35, (H < 1/2), unstable solution, Figure 3(e).

Figure 3.

Numerical simulation of the zero solution for case 3.

The exponential stability criteria we derive are fundamentally stronger than the asymptotic or moment stability conditions found in previous work, as they guarantee not only convergence but exponential decay rates for all moments. These conditions automatically imply almost-sure stability due to the exponential mean square bounds. The resulting stability bounds emerge directly from solving the fractional differential inequality, which makes them naturally tight for the chosen functional class. Our Lyapunov functional candidate was specifically designed to capture both the fractional noise and delay effects through exact weighting terms. Although we have derived sufficient conditions for global exponential mean square stability and asymptotic stochastic stability of the system, the tightness or optimality of these conditions remains an open question. In particular, the gap between sufficiency and necessity has not been closed. Future work could focus on investigating whether these conditions are also close to being necessary or if they can be relaxed without compromising the stability guarantees. Some preliminary numerical evidence suggests that the derived conditions are conservative in certain parameter regimes, especially for extreme values of the Hurst index, indicating potential room for refinement.

The parameters used in the simulations were selected to reflect realistic values in financial and physical systems modeled by fractional Brownian motion. In particular, the Hurst index H ∈ (0, 1) where long-range dependence becomes significant. The coefficients were chosen to satisfy the sufficient conditions derived in our stability theorems. These values were also chosen to highlight different stability scenarios and explore the sensitivity of the system to key parameters, as used in Ref. 5.

The simulation results provide a realistic illustration of how memory effects and delay terms influence stability in systems driven by fractional Brownian motion. For example, higher Hurst index values H > 0.5 often correspond to greater persistence and can destabilize solutions when combined with larger delay or drift parameters. Conversely, smaller H values tend to promote stability under appropriate damping coefficients. These patterns align with behaviors observed in real-world models, such as in finance, where long-range dependence can amplify volatility. The simulation cases also demonstrate that even small changes in parameters can lead to dramatic shifts from stability to instability, underscoring the importance of precise model calibration.

Conclusion and further directions

In the present study, we have investigated the dynamics of the solution of the general nonlinear delay differential equation driven by a fractional Brownian motion. Using appropriate Lyapunov functionals, we have derived the necessary and sufficient criteria for the stability in probability and global exponential mean square stability. The dynamics of the fractional delayed stochastic market model is also studied. We have investigated the impact of the Hurst parameter H ∈ (0, 1) on the behavior of the solution. We have supported our findings with some numerical simulations of the trajectory of the solution. We hope that our work can be applied to many applications in epidemiology and ecology. Moreover, a challenging point for us is to study the mean square stability of systems involving distributed delays.