Lyapunov stability analysis for nonlinear delay systems under random effects and stochastic perturbations with applications in finance and ecology

This manuscript is involved in the study of stability of the solutions of functional differential equations (FDEs) with random coefficients and/or stochastic terms. We focus on the study of different types of stability of random/stochastic functional systems, specifically, stochastic delay differential equations (SDDEs). Introducing appropriate Lyapunov functionals enables us to investigate the necessary conditions for stochastic stability, asymptotic stochastic stability, asymptotic mean square stability, mean square exponential stability, global exponential mean square stability, and practical uniform exponential stability. Some examples with numerical simulations are presented to strengthen the theoretical results. Using our theoretical study, important aspects of epidemiological and ecological mathematical models can be revealed. In ecology, the dynamics of Nicholson’s blowflies equation is studied. Conditions of stochastic stability and stochastic global exponential stability of the equilibrium point at which the blowflies become extinct are investigated. In finance, the dynamics of the Black–Scholes market model driven by a Brownian motion with random variable coefficients and time delay is also studied.


Introduction
Many applications in many disciplines assume that the future state of the system is independent of the past states and is only determined by the current state. This is known as the causality principle which is only an approximation to the true situation of the system.
Consequently, it is practical to consider the past states of the system. A real system should be modeled by differential equations involving delays in time [1,2]. The past dependence has been assumed in the differential equation through the state variable and not its derivative. Hence, this paper focuses on the retarded (delayed) differential equations (DDEs) in © The Author(s) 2021. This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/. the form For h > 0, let C := C([-h, 0], R n ) denote a Banach space of continuous functions: [-h, 0] → R n with respect to the supremum norm φ := sup s∈ [-h,0] |φ(s)|. In (1), x t := x(t + s), -h ≤ s ≤ 0, and the initial history function is Equation (1) with (2) is involved in many interesting mathematical problems. Different delay differential equations are primarily taken from ecological and financial science literature. There are many disciplines such as, but not limited to, population dynamics, ecology, economy, and neural networks [2,3]. Many recent works dedicated to the delay differential equations and the fractional differential equations can be found in [4][5][6][7][8][9][10][11][12][13]. Many numerical studies dedicated to fractional models in epidemiology can be found in [14][15][16][17]. The deterministic models do not preserve the natural uncertainty of the dynamics of the system, whereas the random/stochastic systems preserve all types of the uncertainty. Stochastic delay differential equations can be considered as DEs with random elements and/or stochastic terms with time delays. In this paper, we propose the uncertainty through perturbing the system by the Brownian motion or through assuming that the coefficients, initial conditions, and forcing terms are random variables or stochastic terms. We shall insert the stochastic term involving the Brownian motion, and this technique is a kind of non-parametric stochastic perturbation, then (1)-(2) become ⎧ ⎨ ⎩ dx(t) = μ(t, x t ) dt + σ (t, x t ) dB(t), x t := x(t + s), -h < s ≤ t 0 , This system describes the real world problems in an effective way. SDDEs have been applied in biology, finance, neural networks, ecology, etc. Many works can be found in [18,19]. The probabilistic behavior of x(t) in (3) is restricted to a specific pattern which is Gaussian distribution as the source of randomness is the Wiener process. In (3), μ, σ are continuous functionals defined on [t 0 , ∞) × C( [-h, 0]) and are assumed to satisfy the local Lipschitz condition, i.e., for L > 0, |μ 1 (t, x t )μ 1 (t, x * t )| ≤ L|x tx * t |, obviously, | ∂μ(t,x t ) ∂x | < L, and also for σ . The stochastic process B(t) : [t 0 , ∞) → R is a one-dimensional Brownian motion defined on the complete probability space ( , F, F B t , P), where F B t is the filtration generated by it up to time t [20].
We will also discuss the random delay differential equations (RDDEs) with discrete delay τ > 0 in the form where ω belongs to the underlying probability space ( , F, P). Via the outcome ω, the uncertainty is introduced. This is another way to introduce the uncertainty by assuming that the inputs (coefficients, initial conditions, forcing terms, etc.) are random variables and/or stochastic processes. In this case, we have a wider type of probability distributions such as Gaussian, gamma, beta, binomial, etc. And this provides a great flexibility in dealing with the real life applications. Many applications regarding the RDDEs can be found in [21][22][23][24]. Many authors have studied the existence and uniqueness of the stochastic delay differential equations, for more details, see [25]. The SDDE has a unique global solution x(t, t 0 , φ) for any given initial history function x 0 = φ ∈ C([-h, 0]). Our study assumes two important assumptions as follows.
Assumption 1 Assume that the two functionals μ and σ satisfy the following conditions: Then system (3) admits the trivial solution, and hence stochastic stability, mean square stability, and global mean square exponential stability have received very much attention.
If the origin is not the equilibrium point, then practical stability guarantees the stability and boundedness of all paths in a certain ball of the center at the origin.

Assumption 2
For practical global uniform stability, if μ and σ satisfy μ(t, 0) = 0, σ (t, 0) = 0, then our study of stability of SDDE will be in a small neighborhood (ball B r ) of the origin.
Constructing appropriate Lyapunov functions/functionals is an effective way in the study of stability of deterministic and stochastic systems, more details can be found in [18,25,26]. Introducing suitable Lyapunov functionals implies various types of stability with different stability conditions. This manuscript uses the Lyapunov functionals to investigate the asymptotic stochastic stability, global mean square exponential stability of the zero solution of (3). And practical stability of (3) is also studied when the origin is not necessarily an equilibrium point. It is worthy to study the global exponential mean square stability while many researchers did not. Regardless of the initial history function of the system, global exponential stability makes any trajectory tend to the attractor of the system. Moreover, the paper discusses the asymptotic mean square stability of (3) with random coefficients. Many works regarding the stability of random/stochastic dynamical systems are considered in [27][28][29].
With regard to the applications, the main application in our manuscript is Nicholson's blowflies equation perturbed by white noise. This model is one of the major problems in ecology which describes the dynamics of the population of the Australian sheep blowfly, Lucilia cuprina, which is known as Nicholson's blowfly. The author in [30] introduced the differential equation that models the population dynamics of this blowfly with delay in the form where x(t) is the population of the mature adults at time t. All parameters p, a, δ ∈ [0, ∞), where p represents the maximum daily production rate of eggs per capita, a -1 is the size at which the population reproduces at the maximum rate, δ is the adult death rate of per capita daily, and τ > 0 is the delay in the production process. Model (5) admits only the trivial solution (x = 0) if the basic reproduction ratio of the system R 0 < 1, where R 0 = p δ . The species become extinct due to the global stability of the zero solution of (5), see [31]. Many authors have studied the dynamics of this equation [31][32][33][34][35], and stochastically [36][37][38].
We shall introduce the necessary and sufficient conditions for stochastic stability and stochastic global exponential stability of the trivial equilibrium of (5) under the influence of Brownian motion.
Moreover, we study the Black-Scholes delay market model driven by a Brownian motion with random variable coefficients. This model has become the most known way to model pricing options in financial markets. It was shown firstly in 1973 by [39]. Many authors have studied the stochastic volatility of this model, see [40][41][42]. Generally, the stability of a stochastic fractional Black-Scholes market model has been studied by [43], for example.
For the numerical simulation, this paper uses the algorithm of an Euler-Maruyama scheme which was shown by [44]. The solution of (3) can be written in the integral form The last integral is called an Itô stochastic integral.
. The step size t must be small enough for better numerical approximations. From solution (6), Then, the Euler-Maruyama scheme has the form x(t n ) is the data required in the scheme and x(t n+1 ) is the resulting process at t n+1 . The EM scheme is strongly convergent with order 0.5, i.e., if we want to decrease the error 10 times, the step size should be smaller 100 times. t cannot be too small because of time and the computational errors. [45] studied the convergence rate of this scheme. There are alternative numerical schemes like the Milstein scheme (which is more accurate because of the second order "correction" term added to the scheme), firstly shown by [46], and the Runge-Kutta scheme. With order 1, these schemes are strongly convergent, i.e., by reducing the error 10 times, t must be smaller 10 times. The Runge-Kutta scheme requires computing the derivatives, so we avoid it at present. This manuscript is organized as follows. Section 2 is devoted to some important preliminaries. A rigorous mathematical study of stability with proofs of some important theorems is presented in Sect. 3. Stability of the stochastic Nicholson's blowflies model and the Black-Scholes model with random coefficients is presented in Sect. 4. Some numerical examples with stability regions and computer simulations are shown in Sect. 5. Results and discussions are presented in Sect. 6.

Preliminaries
Notations ( [20, 47-49]) 1. Define S h = {x ∈ R n , x < h} and K as the family of all continuous nondecreasing functions υ : The zero-mean property of the stochastic integral in (6) is useful, so we state without proof this property in the next theorem, and the proof is stated in [50].
Define the differential operator L associated with (3) by This operator can be applied to the Lyapunov function which is called an Itô formula, T is the transposition, and Tr is the trace of the matrix.
which is called Hölder's inequality. In particular, for n = m = 2, we have the Cauchy- and let α, β, T be arbitrary positive numbers, then this is known as the exponential martingale inequality.

Stochastically asymptotically stable if it is stochastically stable and
Asymptotically mean square stable if it is mean square stable and lim t→∞ E|x(t, φ)| 2 = 0. 7. Globally exponentially mean square stable if it is mean square stable and ∀δ > 0, As we mentioned earlier, we will investigate the stability of the solutions with respect to small neighborhood of the origin. Sometimes the state of a system is probably unstable and then the system may oscillate near this state. It is more suitable to address another notion of stability. Practical stability concepts have been introduced by [53][54][55]. Often, the practical stability is referred to as ultimate roundness with a fixed bound.

Definition 2.5
The ball B r is stochastically practically uniformly stable if for each 0 < ε < 1, r ≥ 0, and k > r there exist δ = δ(ε, k) and φ < δ such that (3) is practically uniformly exponentially mean square stable if there exists r > 0 such that the ball B r is practically uniformly exponentially stable a.s., i.e., (3) is practically uniformly exponentially unstable if there exists r > 0 such that the ball B r is practically uniformly exponentially unstable a.s., i.e.,

Stability theory of random/ stochastic DDE
In this section, we show the proofs of stability theorems. Based on Assumption 1, Sects. 3.1, 3.2, and 3.3 investigate theorems of asymptotic stochastic stability, global mean square exponential stability of (3) under the influence of white noise, initial stochastic process, and/or random variable coefficients. Based on Assumption 2, Sect. 3.4 investigates the practical stability of the solution of (3).

Stochastic stability and asymptotic stochastic stability
Assume that 0 < ε < 1 and > 0, one can find by the continuity and the boundedness of and for κ < t we have

Global exponential mean square stability
Theorem 3. 4 The zero solution of (3) is globally exponentially mean square stable.
Then the zero solution is stable in mean square as For mean square exponential stability, for λ > 0, we have For global stability, for t ≥ t 0 , choose a proper constant α ∈ [0, 1) such that Now, define the Lyapunov functional which implies From Itô formula (8) we get The last term vanishes by the zero-mean property of stochastic integral (Theorem 2.1). Therefore, .
Then (16) holds. Consequently, the zero solution of (3) is globally exponentially mean square stable.
Theorem 3.5 Exponential mean square stability implies almost sure exponential stability of the zero solution of (3).

DDE with random/stochastic inputs
Consider the Itô type DDE with discrete delay and an initial function as a stochastic process: The initial function is assumed to be a second-order stochastic process defined on the probability space ( , F, P) and x(t, ω) := x(t) is the solution stochastic process. The continuous functionals are assumed to satisfy the following conditions, for φ 1 (t), φ 2 (t) from C[-h, t 0 ]: with ∞ 0 dr i (s) < ∞, r i (t) for i = 1, 2 are nondecreasing bounded functions.

Theorem 3.6
The zero solution of (18) is asymptotically mean square stable.
Consider the Itô type SDDE with random coefficients a(ω), b(ω) are assumed to be fourth-order random variables defined on the probability space ( , F, P) and φ(t) is an initial deterministic function.

Theorem 3.7
The zero solution of (21) is asymptotically mean square stable if 1. a(ω), b(ω) are fourth-order random variables.

Practical uniform stability and instability
This subsection is devoted to the study of stability of the SDDEs when the origin is not an equilibrium point, we can study the stability of all solution paths in the neighborhood of a ball that is centered at the origin. υ 1 , υ 2 , υ 3 , and V(t, ·) are required to be defined in the neighborhood of the origin. Under Assumption 2, and the ball B r := {x ∈ R : x(t) ≤ r}, we state the following theorems.

Theorem 3.10 The solution of system (3) is almost sure practically uniformly exponentially unstable if we reverse the inequalities in the previous theorem.
Proof The same as in the previous theorem, one gets The quadratic variation of by the law of large numbers for martingales. Therefore, and this proves the a.s. globally uniformly exponential instability in mean square for 2c 2 ≥ c 1 .

Stochastic Nicholson's blowflies model
System (5) is exposed to stochastic perturbation in the form of white noise which is assumed to be proportional to the deviation of the current state of the system from the zero solution. Therefore, the stochastic version of (5) will be in the form The parameter σ represents the intensity of the noise. Here, we will investigate the necessary conditions of mean square stability of the zero solution of the corresponding linear system which are the same conditions of stochastic stability of the nonlinear system (24). Moreover, stochastic global exponential stability of the zero equilibrium is discussed.
, then the zero solution of (24) is stochastically stable.
Then, by the induction and (30), we have to prove the global stability by proving This can be done by following the same argument of Theorem 3.4 with |α| = | p δ |. Therefore, the zero solution of (25) is globally exponentially mean square stable under condition (29) which is the same condition of stochastic stability of the nonlinear system (24).

The random delay Black-Scholes market model
Consider the stochastic Black-Scholes model driven by random variable coefficients in the forṁ this is a model of volatility in finance of the process of asset prices. The process {x(t), t ≥ 0} is the price of a stock with a nonnegative random variable drift A(ω), the volatility random variable σ (ω), and a Wiener process B(t).

Proposition 4.3 Equation (32) is:
1. Stable in probability if only A(ω) is a nonnegative random variable.

Asymptotically stable in probability if
Proof Consider the random Lyapunov function V(t, x(t, ω)) = x 2 (t, ω)e -λht with λ > 0. Then For stochastic stability, the last term vanishes for sufficiently large λ in (33). Therefore, according to inequality (11) of Theorem 3.1, LV(t, x t , ω) ≤ 0 and equation (32) is stochastically stable if only A(ω) is a nonnegative random variable. For asymptotic stochastic stability, consider the Lyapunov functional V(t, It is enough to prove ELV ≤ -υ( x(t) ), υ( x(t) ) ∈ K according to Theorem 3.2. It is easy to verify the decrescent property of the considered Lyapunov functional as for M 1 , M 2 > 0, Then, according to Theorem 3.7, A(ω) 2 + σ σ (ω) 2 4 < λh -1 is a condition of the asymptotic stochastic stability which is the same condition of mean square exponential stability according to Theorem 3.4. Now, from Theorem 3.3, we can see that V is radially unbounded as Therefore, (32) is stable in the large. For mean square stability, consider the Lyapunov The condition of mean square stability is σ (ω) 2 4 < 2 A(ω) 2 according to Theorems 3.4, 3.7.
For practical stability, consider the Lyapunov functional V(t, x t , ω) = x 2 (t, ω). Then Hence, And According to the conditions of Theorem 3.9 and the conditions on the random variables in Theorem 3.7, the constants c 1 = 0 and c 2 = -2 A(ω) 2 . Then Consequently, for i.e., if A(ω) 2 > 0.5, the solution of (32) is globally practically uniformly exponentially stable.
By introducing a different Lyapunov functional a different stability condition is obtained in the form σ 2 2 + |b|( |b|h 2 + 1) + c < 0 which is the delay-dependent condition.  Fig. 2(a), and also when c = -7 in Fig. 2(b). Better stability regions are obtained by decreasing c with the delay h. Numerical simulation of the solution of (34) is shown in Fig. 3 for h = 2, all trajectories of the solution are convergent to the zero solution.
Example 2 Consider the stochastic scalar equation with two discrete time delays h, τ > 0 The Lyapunov functional leads to a + |b| + σ 2 2 < 0 as a condition of the mean square asymptotic stability of the zero solution. Stability regions are shown in Fig. 7(a) for different values of a, (1) a = -2, (2) a = -3, (3) a = -5, (4) a = -10, and (5) a = -15, and the numerical simulation of the zero solution is shown in Fig. 7(b).
Example 3 Consider the two-dimensional system of stochastic delay differential equation B(t) = (B 1 (t), B 2 (t)) T is a standard Wiener process, x(t) = (x 1 (t), x 2 (t)) T is a two-dimensional vector function where T is the transposition. The system in the matrix form The Lyapunov functional V(t, Hence, the zero solution of (38) is asymptotically mean square stable if a < -1, b < -1 for arbitrary constants a, b. Fig. 10 gives the mean square stability region in a (a, b)-space of parameters and the numerical simulation of the solution is shown in Fig. 11. Consider the random system The Lyapunov functional x 2 2 (s, ω) ds Figure 11 Ten stable and unstable (unbounded) trajectories for each process x 1 (t) and x 2 (t) of the solution of (38) Figure 12 Ten stable solutions for each process x 1 (t) and x 2 (t) of (39) Figure 13 Ten unstable solutions for each process x 1 (t) and x 2 (t) of (39) leads to the conditions of asymptotic mean square stability in the form 2a + A(ω) 2 + B(ω) 2 4 < 0 and 2b + A(ω) 2 + C(ω) 2 4 < 0, A(ω) is a nonnegative random variable. Figs. 12, 13 give the numerical simulation of the solutions for different random probability distributions.
Example 4 Consider a stochastic one-dimensional Nicholson's blowflies equation (24) with σ = 0.8 and φ(s) = 0.1 cos(s). In Fig. 14 and according to condition (27), the stochastic stability regions are obtained for different values of the delay τ in a (p, δ)-space of parameters. Numerical simulations of the solution are obtained in Fig. 15 for p = 0.5, δ = 1.2, and τ = 0.1.

Figure 14
Regions of stochastic stability of the zero solution of (24) Figure 15 Numerical simulations of the solution of (24) Figure 16 Practical uniform exponential solution of (32) Example 5 Consider equation (32) with A(ω) ∼ Bin(1, 0.6), σ (ω) ∼ Bin(1, 0.001), and τ = 0.001. Fig. 16 shows that the trajectories are bounded and remain in a sufficiently small neighborhood of the origin for small enough τ .

Results and discussions
1. Stability of the zero solution of a delay differential equation driven by a Brownian motion and/or random variable coefficients has been studied through Theorems 3.1-3.7. Asymptotic stochastic stability and global mean square exponential stability are investigated. Examples 1-4 support the theoretical results. Introducing more Lyapunov functionals implies various types of stability with different stability conditions. By Lyapunov functionals, we have obtained delay-dependent and delay-independent conditions. Our study shows that the delay-dependent conditions are very important and have a great effect on the regions of stability. SDDEs involving random coefficients have been also discussed. Examples in one and two dimensions are shown. Some computer simulations are carried out for various random probability distributions. These random variables provide various patterns to the solution stochastic process and make the model more attractive. The main application in this paper is Nicholson's blowflies model. Propositions 4.1, 4.2 have studied the stochastic stability and stochastic global exponential stability for the zero solution, respectively. Regions of stability of the zero solution at which these blowflies become extinct are shown. These regions become better for small values of delay τ . Regarding the financial market application, Proposition 4.3 has studied the stability of a stochastic Black-Scholes model with random variable coefficients. 2. Theorems 3.8, 3.9 have studied the practical stability of the solution when the origin is not necessarily an equilibrium point. By these theorems, we have concluded that all trajectories of the solution of the Black-Scholes model are practically uniformly exponentially stable. According to Example 5, for small enough τ and A(ω) 2 ≥ 1 2 , almost all trajectories are bounded and remain in a sufficiently small neighborhood of the origin. By stochastic mathematical epidemiological and ecological models, we can reveal the mechanisms that influence the transmission and control the disease. Stochastic models are proposed to capture the uncertainty and variations in the transmission of the diseases. Many recent works dedicated to stochastic epidemiological and ecological models can be found in [37,38,40,41,50,52,[58][59][60][61]. There is still considerable uncertainty with regard to the study of stability, therefore we have addressed the stochastic global exponential and practical stability. Mathematically, we have focused on the qualitative theory (stability of equilibria) of stochastic/random systems in a rigorous way. From ecological point of view, dynamics of stochastic Nicholson's blowflies models is studied. We claim that our approach can be extended to various biological, ecological, and economical models. This approach will reveal some important aspects of these models.

Conclusion and further directions
In this paper, we have studied the stability of Nicholson's blowflies model and Black-Scholes market model based on some important proven theorems in asymptotic stochastic stability, global exponential mean square stability, and practical uniform exponential stability. Our results on Nicholson's blowflies model can be extended to its general form known as the stochastic neoclassical growth model, for γ > 0, γ = 1 dx(t) = px γ (tτ )e -ax(t-τ )δx(t) dt + σ x(t) dB(t).
Moreover, we hope that our findings will be useful for models driven by fractional Brownian motion (FBM).