Asymptotic Behavior of Three Connected Stochastic Delay Neoclassical Growth Systems Using Spectral Technique

Abstract

In this study, we consider a nonlinear system of three connected delay differential neoclassical growth models along with stochastic effect and additive white noise, which is influenced by stochastic perturbation. We derived the conditions for positive equilibria, stability and positive solutions of the stochastic system. It is observed that when a constant delay reaches a certain threshold for the steady state, the asymptotic stability is lost, and the Hopf bifurcation occurs. In the case of the finite domain, the three connected, delayed systems will not collapse to infinity but will be bounded ultimately. A Legendre spectral collocation method is used for the numerical simulations. Moreover, a comparison of a stochastic delayed system with a deterministic delayed system is also provided. Some numerical test problems are presented to illustrate the effectiveness of the theoretical results. Numerical results further illustrate the obtained stability regions and behavior of stable and unstable solutions of the proposed system.

1. Introduction

In mathematical economics, the examination of the stochastic delay differential neoclassical growth model (NGM) plays a key role. In general, this model is constructed with two very simple assumptions; one is capital and full-time labor hiring, while the other is the immediate adjustment in the market, which helps in the long-run behavior of the economy [1,2,3]. The main advantage of these models is that they are well-behaved and are usually asymptotically stable for the steady state, but in reality, these growth path models constantly exhibit fluctuations. For this reason, the neoclassical model could be a good alternative to show how such persistent behavior can emerge when nonlinearities and a production delay are present. Since the neoclassical growth model is always affected by environmental noises, the stochastic models are more suitable in the real world [4,5,6,7].

In economics, the most frequently discussed issue is to test the economic growth models. Many researchers have investigated these models for various population data and complex behaviors. Day studied a neoclassical growth model with time delay and noticed that despite its simple structure, the resulting dynamic system shows the emergence of erratic fluctuations in the capital accumulation process when the production function is unimodal and the delay in production is explicitly considered [8,9]. However, his models were totally occupying discrete time and a mound-shaped function that described the negative effect of subsequent pollution from increasing fundamentals. It was identified by numerical approaches that such models could achieve periodic and even chaotic behavior. Following the pioneering work of Day, Matsumoto and Szidarovszky, an economics-based model for understanding the complex dynamics of economics was created [10,11,12]. In earlier work in this field, most of the researchers only considered discrete time scales [13,14,15]. For the detection of chaos, the three period condition, introduced by Li and Yorke, has many applications in nonlinear differential equations of the first-order, followed by the work of Rosser, which offers many applications [16,17]. Very little work has been performed that is committed to the case of continuous time scales due to the fact that there is no preferred criterion to detect chaos and the system has three dimensions.

In this article, we will examine an extension of the NGM to the early works of Swan and Solow [18,19]. The new NGM is constructed with three connected assumptions; one is the permanent labor employment, the second is the continual adjustment in the market output and the third one is the instantaneous growth of the products. Thus, it is very convenient to describe the long-term behavior of the economy due to the well-operating production function. We provide a detailed stability analysis of the steady state in the continuous time structure with time delays. We further investigate the equilibrium points of a system of three connected NGMs, positive equilibria and conditions for stability with stochastic-type effects that are directly proportional to the obtained equilibrium from deviation of the system state. For the numerical simulations, we use spectral methods based on Legendre polynomials [20,21,22,23,24,25,26,27].

The remaining structure of the article is: In Section 2, the mathematical model is formulated with time delays, followed by the description of the method in Section 3. Some preliminary results are given in Section 4. A stability analysis is presented in Section 5. For the confirmation of theoretical results, some numerical tests are performed in Section 6, and Section 7 concludes the article.

2. Model Description

To study the three connected NGMs and discuss the stability of the zero equilibrium under stochastic effects, the stochastic three connected NGMs have the following form:

$$\begin{array}{c}\hfill \left\{\begin{array}{c}d{x}_1\left(t\right)=\left[-{\alpha }_1{x}_1\left(t\right)+{\beta }_1{x}_2\left(t\right)+{\gamma }_1{x}_3\left(t\right)+{\delta }_1{x}_1^{{\nu }_1}(t-{\tau }_1)e^{-{\rho }_1{x}_1(t-{\tau }_1)}\right]dt+{\sigma }_1{x}_1\left(t\right)d{B}_1\left(t\right)\hfill \\ d{x}_2\left(t\right)=\left[-{\alpha }_2{x}_2\left(t\right)+{\beta }_2{x}_3\left(t\right)+{\gamma }_2{x}_1\left(t\right)+{\delta }_2{x}_2^{{\nu }_2}(t-{\tau }_2)e^{-{\rho }_2{x}_2(t-{\tau }_2)}\right]dt+{\sigma }_2{x}_2\left(t\right)d{B}_2\left(t\right)\hfill \\ d{x}_3\left(t\right)=\left[-{\alpha }_3{x}_3\left(t\right)+{\beta }_3{x}_1\left(t\right)+{\gamma }_3{x}_2\left(t\right)+{\delta }_3{x}_3^{{\nu }_3}(t-{\tau }_3)e^{-{\rho }_3{x}_3(t-{\tau }_3)}\right]dt+{\sigma }_3{x}_3\left(t\right)d{B}_3\left(t\right),\hfill \end{array}\right.\end{array}$$

(1)

with initial values:

$$\begin{array}{c}\hfill x_i\left(s\right)={\varrho }_i\left(s\right);s\in [-\tau ,0];\varrho \in C\left([-\tau ,0],{\mathbb{R}}_{+}\right);i=1,2,3.\end{array}$$

(2)

x is the capital per labor, where ${\mathbb{R}}_{+}=(0,+\infty ),$ and ${\alpha }_i,(i=1,2,3)$ are each positive. Moreover, ${\beta }_i$ and ${\gamma }_i$ are the coupling coefficients, where all the remaining parameters ${\delta }_i,{\nu }_i,{\rho }_i$ and $\tau =max\{{\tau }_1,{\tau }_2,{\tau }_3\}$, are greater than zero. ${\beta }_i\left(t\right)(i=1,2,3)$ are independent white noises and ${\sigma }_i^2(i=1,2,3)$ denote noises intensities. For brief details of the above parameters backdrop, we refer the readers to [28]. The neoclassical growth differential system with a delay and with variable coefficients is investigated in [29,30,31]. Shaikhet studies the two connected NGMs with stochastic perturbation and investigates the stability of equilibrium [32]. Some research work related to the stochastic delay system, stochastic fractional delay system, stochastic complex network with delay and stochastic highly nonlinear coupled system with delays can be found in [33,34,35,36,37]. In the literature, to the best of our knowledge, no one has considered the three connected stochastic systems. The main motivation of this work is to consider the three connected NGMs and to apply a high-order numerical scheme based on Legendre polynomials along with theoretical justifications.

3. Description of the Method

This section incorporated the spectral method (SM) for solving the stochastic neoclassical growth model given by Equation (1). In the present method, we used Legendre Gauss quadrature along with the weight function. For the SM, we consider the Legendre Gauss Lobatto points ${\left\{t_j\right\}}_{j=0}^N$.

Our aim in this study is to develop an approximate solution to Equation (1). We apply the integral of Equation (1) from $[0,t]$, then:

$$\begin{array}{cc}\hfill x_1\left(t\right)& =x_1\left(0\right)+{\int }_0^t\left(-{\alpha }_1{x}_1\left(s\right)+{\beta }_1{x}_2\left(s\right)+{\gamma }_1{x}_3\left(s\right)+{\delta }_1{x}_1^{{\nu }_1}(s-{\tau }_1)e^{-{\rho }_1{x}_1(s-{\tau }_1)}\right)ds\hfill \\ \hfill & +{\int }_0^t{\sigma }_1{x}_1\left(s\right)dB\left(s\right),\hfill \\ \hfill x_2\left(t\right)& =x_2\left(0\right)+{\int }_0^t\left(-{\alpha }_2{x}_2\left(s\right)+{\beta }_2{x}_3\left(s\right)+{\gamma }_2{x}_1\left(s\right)+{\delta }_2{x}_2^{{\nu }_2}(s-{\tau }_2)e^{-{\rho }_2{x}_2(s-{\tau }_2)}\right)ds\hfill \\ \hfill & +{\int }_0^t{\sigma }_2{x}_2\left(s\right)dB\left(s\right),\hfill \\ \hfill x_3\left(t\right)& =x_3\left(0\right)+{\int }_0^t\left(-{\alpha }_3{x}_3\left(s\right)+{\beta }_3{x}_1\left(s\right)+{\gamma }_3{x}_2\left(s\right)+{\delta }_3{x}_3^{{\nu }_3}(s-{\tau }_3)e^{-{\rho }_3{x}_3(s-{\tau }_3)}\right)ds\hfill \\ \hfill & +{\int }_0^t{\sigma }_3{x}_3\left(s\right)dB\left(s\right),\hfill \end{array}$$

(3)

where $x_1\left(0\right),x_2\left(0\right)$ and $x_3\left(0\right)$ are the initial values for the functions $x_1\left(t\right),x_2\left(t\right)$ and $x_3\left(t\right)$, respectively. Taking linear transformation $s=\frac{t}{2}(1+\theta )=\eta$(say) to analyze the SM over standard interval $[-1,1]$ in Equation (3):

$$\begin{array}{cc}\hfill x_1\left(t\right)& =x_1\left(0\right)+\frac{t}{2}{\int }_{-1}^1\left(-{\alpha }_1{x}_1\left(\eta \right)+{\beta }_1{x}_2\left(\eta \right)+{\gamma }_1{x}_3\left(\eta \right)+{\delta }_1{x}_1^{{\nu }_1}(\eta -{\tau }_1)e^{-{\rho }_1{x}_1(\eta -{\tau }_1)}\right)d\theta \hfill \\ \hfill & +\frac{t}{2}{\int }_{-1}^1{\sigma }_1{x}_1\left(\eta \right)dB\left(\theta \right),\hfill \\ \hfill x_2\left(t\right)& =x_2\left(0\right)+\frac{t}{2}{\int }_{-1}^1\left(-{\alpha }_2{x}_2\left(\eta \right)+{\beta }_2{x}_3\left(\eta \right)+{\gamma }_2{x}_1\left(\eta \right)+{\delta }_2{x}_2^{{\nu }_2}(\eta -{\tau }_2)e^{-{\rho }_2{x}_2(\eta -{\tau }_2)}\right)d\theta \hfill \\ \hfill & +\frac{t}{2}{\int }_{-1}^1{\sigma }_2{x}_2\left(\eta \right)dB\left(\theta \right),\hfill \\ \hfill x_3\left(t\right)& =x_3\left(0\right)+\frac{t}{2}{\int }_{-1}^1\left(-{\alpha }_3{x}_3\left(\eta \right)+{\beta }_3{x}_1\left(\eta \right)+{\gamma }_3{x}_2\left(\eta \right)+{\delta }_3{x}_3^{{\nu }_3}(\eta -{\tau }_3)e^{-{\rho }_3{x}_3(\eta -{\tau }_3)}\right)d\theta \hfill \\ \hfill & +\frac{t}{2}{\int }_{-1}^1{\sigma }_3{x}_3\left(\eta \right)dB\left(\theta \right),\hfill \end{array}$$

(4)

The spectral equations (semi-discretised) form of Equation (4) is given by

$$\begin{array}{cc}\hfill x_1\left(t\right)& =x_1\left(0\right)+\frac{t}{2}\sum _{k=0}^N\left(-{\alpha }_1{x}_1\left(\eta \right)+{\beta }_1{x}_2\left(\eta \right)+{\gamma }_1{x}_3\left(\eta \right)+{\delta }_1{x}_1^{{\nu }_1}(\eta -{\tau }_1)e^{-{\rho }_1{x}_1(\eta -{\tau }_1)}\right){\omega }_k\hfill \\ \hfill & +\frac{t}{2}\sum _{k=0}^N{\sigma }_1{x}_1\left(\eta \right){\omega }_k^{*},\hfill \\ \hfill x_2\left(t\right)& =x_2\left(0\right)+\frac{t}{2}\sum _{k=0}^N\left(-{\alpha }_2{x}_2\left(\eta \right)+{\beta }_2{x}_3\left(\eta \right)+{\gamma }_2{x}_1\left(\eta \right)+{\delta }_2{x}_2^{{\nu }_2}(\eta -{\tau }_2)e^{-{\rho }_2{x}_2(\eta -{\tau }_2)}\right){\omega }_k\hfill \\ \hfill & +\frac{t}{2}\sum _{k=0}^N{\sigma }_2{x}_2\left(\eta \right){\omega }_k^{*},\hfill \\ \hfill x_3\left(t\right)& =x_3\left(0\right)+\frac{t}{2}\sum _{k=0}^N\left(-{\alpha }_3{x}_3\left(\eta \right)+{\beta }_3{x}_1\left(\eta \right)+{\gamma }_3{x}_2\left(\eta \right)+{\delta }_3{x}_3^{{\nu }_3}(\eta -{\tau }_3)e^{-{\rho }_3{x}_3(\eta -{\tau }_3)}\right){\omega }_k\hfill \\ \hfill & +\frac{t}{2}\sum _{k=0}^N{\sigma }_3{x}_3\left(\eta \right){\omega }_k^{*},\hfill \end{array}$$

(5)

where the Legendre–Gauss quadrature with weights are

$$\begin{array}{c}\hfill {\omega }_k=\frac{2}{{\left[L_{N+1}^{{}^{\prime }}\left(s_k\right)\right]}^2(1-s_k^2)},0\le k\le N.\end{array}$$

Similarly, ${\omega }_k^{*}=\sqrt{{\omega }_k}\times randn(1,N)$ is the stochastic weight function.

To find the numerical solution for the proposed system, we used the Legendre polynomials of the following form:

$$\begin{array}{c}\hfill x_1\left(t\right)=\sum _{n=0}^N{a}_n{P}_n\left(t\right),x_2\left(t\right)=\sum _{n=0}^N{b}_n{P}_n\left(t\right),x_3\left(t\right)=\sum _{n=0}^N{c}_n{P}_n\left(t\right)\end{array}$$

(6)

In the above equation, $a_n,b_n,c_n$ are the Legendre coefficients for the classes $x_1,x_2,x_3$, respectively, where $P_n\left(t\right)$ are the Legendre polynomials. Incorporating Equation (6) into Equation (5), we get the following algebraic system

$$\begin{array}{cc}\hfill \sum _{n=0}^N{a}_n{P}_n\left(t\right)& =\sum _{n=0}^N{a}_n{P}_n\left(0\right)+\frac{t}{2}\sum _{k=0}^N(-{\alpha }_1\sum _{n=0}^N{a}_n{P}_n\left(\eta \right)+{\beta }_1\sum _{n=0}^N{b}_n{P}_n\left(\eta \right)+{\gamma }_1\sum _{n=0}^N{c}_n{P}_n\left(\eta \right)\hfill \\ \hfill & +{\delta }_1\sum _{n=0}^N{a}_n^{{\nu }_1}{P}_n^{{\nu }_1}(\eta -{\tau }_1)e^{-{\rho }_1{\sum }_{n=0}^N{a}_n{P}_n(\eta -{\tau }_1)}){\omega }_k+\frac{t}{2}\sum _{k=0}^N{\sigma }_1\sum _{n=0}^N{a}_n{P}_n\left(\eta \right){\omega }_k^{*},\hfill \\ \hfill \sum _{n=0}^N{b}_n{P}_n\left(t\right)& =\sum _{n=0}^N{b}_n{P}_n\left(0\right)+\frac{t}{2}\sum _{k=0}^N(-{\alpha }_2\sum _{n=0}^N{b}_n{P}_n\left(\eta \right)+{\beta }_2\sum _{n=0}^N{c}_n{P}_n\left(\eta \right)+{\gamma }_2\sum _{n=0}^N{a}_n{P}_n\left(\eta \right)\hfill \\ \hfill & +{\delta }_2\sum _{n=0}^N{b}_n^{{\nu }_2}{P}_n^{{\nu }_2}(\eta -{\tau }_2)e^{-{\rho }_2{\sum }_{n=0}^N{b}_n{P}_n(\eta -{\tau }_2)}){\omega }_k+\frac{t}{2}\sum _{k=0}^N{\sigma }_2\sum _{n=0}^N{b}_n{P}_n\left(\eta \right){\omega }_k^{*},\hfill \\ \hfill \sum _{n=0}^N{c}_n{P}_n\left(t\right)& =\sum _{n=0}^N{c}_n{P}_n\left(0\right)+\frac{t}{2}\sum _{k=0}^N(-{\alpha }_3\sum _{n=0}^N{c}_n{P}_n\left(\eta \right)+{\beta }_3\sum _{n=0}^N{a}_n{P}_n\left(t\right)\left(\eta \right)+{\gamma }_3\sum _{n=0}^N{b}_n{P}_n\left(\eta \right)\hfill \\ \hfill & +{\delta }_3\sum _{n=0}^N{c}_n^{{\nu }_3}{P}_n^{{\nu }_3}(\eta -{\tau }_3)e^{-{\rho }_3{\sum }_{n=0}^N{c}_n{P}_n(\eta -{\tau }_3)}){\omega }_k+\frac{t}{2}\sum _{k=0}^N{\sigma }_3\sum _{n=0}^N{c}_n{P}_n\left(\eta \right){\omega }_k^{*}.\hfill \end{array}$$

(7)

Thus there is $3N+3$ unknowns in the system given in Equation (7) with $3N$ nonlinear algebraic equations. After incorporating the initial conditions, we get

$$\begin{array}{c}\hfill \sum _{n=0}^N{a}_n{P}_n\left(0\right)=\sum _{n=0}^N{\left({\varrho }_1\right)}_n,\sum _{n=0}^N{b}_n{P}_n\left(0\right)=\sum _{n=0}^N{\left({\varrho }_2\right)}_n\sum _{n=0}^N{c}_n{P}_n\left(0\right)=\sum _{n=0}^N{\left({\varrho }_3\right)}_n.\end{array}$$

(8)

Now, using Equation (7) along with Equation (8) results in $3N+3$ nonlinear equations having $3N+3$ unknowns. We obtain the numerical solution to the proposed stochastic system given in Equation (1) by incorporating the values of these unknowns into Equation (6).

4. Preliminary Results

In the current section, we recommend a few fundamental lemmas and definitions, which might be useful for showing the continuation of the unique global positive solution of Equation (1).

Definition 1.

The proposed system in Equation (1) is bounded in the mean if for each positive $M>0$ free from the initial conditions of Equation (2) as

$$\begin{array}{c}\hfill \underset{t\to \infty }{lim}sup\mathbb{E}\left|x\left(t\right)\right|\le M\end{array}$$

(9)

Lemma 1.

Let $\nu ,\rho >0,$ and $f\left(x\right)=x^{\nu }{e}^{-\rho x},$ then $f\left(x\right)\le {\left(\frac{\nu }{\rho e}\right)}^{\nu }$ for $x\in {\mathbb{R}}_{+}.$

Proof. 

The proof is simple and is left for the reader. □

Lemma 2.

If $a_i\in \mathbb{R},b_i,c_i\in {\mathbb{R}}_{+},(i=1,2,3)$, then $\frac{a_1{x}^2+(b_1+c_1)x+a_2{y}^2+(b_2+c_2)y+a_3{z}^2+(b_3+c_3)z}{1+x^2+y^2+z^2}$

$\le D(a_1,a_2,a_3,b_1,b_2,b_3,c_1,c_2,c_3)$ where $D(a_1,a_2,a_3,b_1,b_2,b_3,c_1,c_2,c_3)=$

$$\begin{array}{c}\hfill \left\{\begin{array}{cc}\left(a_1+\sqrt{a_1^2+b_1^2+c_1^2}+a_2+\sqrt{a_2^2+b_2^2+c_2^2}+a_3+\sqrt{a_3^2+b_3^2+c_3^2}\right)/2,\hfill & a_1,a_2,a_3\ge 0,\hfill \\ -(b_1^2+c_1^2)/4a_1-(b_2^2+c_2^2)/4a_2-(b_3^2+c_3^2)/4a_3,\hfill & a_1,a_2,a_3<0,\hfill \\ \left(a_1+\sqrt{a_1^2+b_1^2+c_1^2}\right)/2-(b_2^2+c_2^2)/4a_2-(b_3^2+c_3^2)/4a_3,\hfill & a_1\ge 0;a_2,a_3<0,\hfill \\ \left(a_2+\sqrt{a_2^2+b_2^2+c_2^2}+a_3+\sqrt{a_3^2+b_3^2+c_3^2}\right)/2-(b_1^2+c_1^2)/4a_1,\hfill & a_1<0;a_2,a_3\ge 0.\hfill \end{array}\right.\end{array}$$

Proof. 

By using Lemma 1.2 of [38] for the two connected neoclassical models, we can obtain the result easily for three connected neoclassicals, so we discard the proof. □

Lemma 3.

For any given initial conditions of Equation (2), there exists a unique global positive solution $x\left(t\right)=\left(x_1\left(t\right),x_2\left(t\right),x_3\left(t\right)\right)$ of Equation (1) in a closed interval $[-\tau ,+\infty ]$, and each $x_i\left(t\right),(i=1,2,3)$ will be a positive with unit probability.

Proof. 

It is simple to see that for $t\in [0,\tau ]$, then the proposed system given in Equation (1) along with the initial conditions of Equation (2) reduces to the linear stochastic system, now by using Theorem 3.3.1 of [39], provided that there is a unique stable solution $x\left(t\right)$ in the interval $[0,\tau ]$: if solution $x\left(t\right)$ is in the interval $[0,\tau ]$ once it is known, then we can easily proceed such arguments in the intervals $[\tau ,2\tau ]$, $[2\tau ,3\tau ]$... Therefore, we will obtain the solution of the max interval $[-\tau ,{\mu }_e]$, where ${\mu }_e$ denotes the explosion time. Now, to prove ${\mu }_e=\infty$, we assume that $m_0\ge 1,$ is a sufficiently large number, such as:

$$\begin{array}{c}\hfill \frac{1}{m_0}<\underset{\tau \le t\le 0}{min}{\varrho }_i\left(t\right)\le \underset{\tau \le t\le 0}{max}{\varrho }_i\left(t\right)<m_0.\end{array}$$

Therefore, for each integer $m\ge m_0$, the stopping time is defined by:

$$\begin{array}{c}\hfill {\mu }_m=inf\left\{t\in [0,{\mu }_e):x_i\left(t\right)\in \left(\frac{1}{m},m\right),i=1,2,3\right\},\end{array}$$

where we assume that $\varphi$ is the empty set with the usual convention $inf\varphi =+\infty$. Obviously, ${\mu }_m$ is consistently increasing as $m\to \infty$. We set ${\mu }_{\infty }={lim}_{m\to \infty }{\mu }_m$, where ${\mu }_{\infty }\le {\mu }_e$. If ${\mu }_{\infty }=\infty$ can be proven, then ${\mu }_e=\infty$ where $x_i\left(t\right)\in {\mathbb{R}}_{+}i=1,2,3$ as $t\ge 0$. For this we need to prove that ${\mu }_{\infty }=\infty$. To do this, we must define $C^2$-function $V:{\mathbb{R}}_{+}\times {\mathbb{R}}_{+}\times {\mathbb{R}}_{+}\to {\mathbb{R}}_{+}$ by $V(x_1,x_2,x_3)={\sum }_{i=1}^3(x_i-1-ln{x}_i)$. For $t\in [0,{\mu }_m\wedge T)$ to show this, we use the It$\widehat{o}$ formula:

$$\begin{array}{cc}\hfill dV\left(x_1\left(t\right),x_2\left(t\right),x_3\left(t\right)\right)& =LV\left(x_1\left(t\right),x_2\left(t\right),x_3\left(t\right),x_1(t-{\tau }_1),x_2(t-{\tau }_2),x_3(t-{\tau }_3)\right)dt\hfill \\ \hfill & +\sum _{i=1}^3{\sigma }_i\left(x_i\left(t\right)-1\right)d{B}_i\left(t\right),\hfill \end{array}$$

(10)

where $m\ge m_0,T>0$ is arbitrary, and the operator’s LV is defined by

$$\begin{array}{cc}\hfill & LV\left(x_1\left(t\right),x_2\left(t\right),x_3\left(t\right),x_1(t-{\tau }_1),x_2(t-{\tau }_2),x_3(t-{\tau }_3)\right)\hfill \\ \hfill & =\sum _{i=1}^3\left({\alpha }_i+\frac{{\sigma }_i^2}{2}+{\delta }_i{x}_i^{{\nu }_i}(t-{\tau }_i)e^{-\rho x_i(t-{\tau }_i)}-\frac{{\delta }_i{x}_i^{{\nu }_i}(t-{\tau }_i)e^{-\rho x_i(t-{\tau }_i)}}{x_i\left(t\right)}\right)\hfill \\ \hfill & -\left({\alpha }_1-({\beta }_3+{\gamma }_2)\right)x_1\left(t\right)-\left({\alpha }_2-({\beta }_1+{\gamma }_3)\right)x_2\left(t\right)-\left({\alpha }_3-({\beta }_2+{\gamma }_1)\right)x_3\left(t\right)\hfill \\ \hfill & -\frac{({\beta }_1+{\gamma }_3)x_2^2\left(t\right)+({\beta }_2+{\gamma }_1)x_3^2\left(t\right)+({\beta }_3+{\gamma }_2)x_1^2\left(t\right)}{x_1\left(t\right)x_2\left(t\right)x_3\left(t\right)}.\hfill \end{array}$$

(11)

We use the inequality $y\le 3\left(y-1-lny\right)+3$ for all $y\in {\mathbb{R}}_{+}$, along with Lemma 1, then we can find Equation (11):

$$\begin{array}{cc}\hfill & LV\left(x_1\left(t\right),x_2\left(t\right),x_3\left(t\right),x_1(t-{\tau }_1),x_2(t-{\tau }_2),x_3(t-{\tau }_3)\right)\hfill \\ \hfill & \le \sum _{i=1}^3\left({\alpha }_i+\frac{{\sigma }_i^2}{2}+{\delta }_i{\left(\frac{{\nu }_i}{{\rho }_{i}e}\right)}^{{\nu }_i}\right)+6max\left\{|{\alpha }_1-({\gamma }_2+{\beta }_3)|,|{\alpha }_2-({\beta }_1+{\gamma }_3)|,|{\alpha }_3-({\gamma }_1+{\beta }_2)|\right\}\hfill \\ \hfill & +3max\left\{|{\alpha }_1-({\gamma }_2+{\beta }_3)|,|{\alpha }_2-({\beta }_1+{\gamma }_3)|,|{\alpha }_3-({\gamma }_1+{\beta }_2)|\right\}V\left(x_1\left(t\right),x_2\left(t\right),x_3\left(t\right)\right)\hfill \\ \hfill & =3max\left\{|{\alpha }_1-({\gamma }_2+{\beta }_3)|,|{\alpha }_2-({\beta }_1+{\gamma }_3)|,|{\alpha }_3-({\gamma }_1+{\beta }_2)|\right\}V\left(x_1\left(t\right),x_2\left(t\right),x_3\left(t\right)\right)+L,\hfill \end{array}$$

(12)

where $L={\sum }_{i=1}^3\left({\alpha }_i+\frac{{\sigma }_i^2}{2}+{\delta }_i{\left(\frac{{\nu }_i}{{\rho }_{i}e}\right)}^{{\nu }_i}\right)+6max\left\{|{\alpha }_1-({\gamma }_2+{\beta }_3)|,|{\alpha }_2-({\beta }_1+{\gamma }_3)|,|{\alpha }_3-({\gamma }_1+{\beta }_2)|\right\}.$ We assume that each $m\ge m_0$ applies integrals on both sides of Equation (10) from 0 to ${\mu }_m\wedge T$, then

$$\begin{array}{cc}\hfill & \mathbb{E}V\left(x_1({\mu }_m\wedge T),x_2({\mu }_m\wedge T),x_3({\mu }_m\wedge T)\right)\hfill \\ \hfill & \le L_1+3max\left\{|{\alpha }_1-({\gamma }_2+{\beta }_3)|,|{\alpha }_2-({\beta }_1+{\gamma }_3)|,|{\alpha }_3-({\gamma }_1+{\beta }_2)|\right\}\hfill \\ \hfill & \times \mathbb{E}{\int }_0^{{\mu }_m\wedge T}V\left(x_1\left(t\right),x_2\left(t\right),x_3\left(t\right)\right)dt\hfill \\ \hfill & \le L_1+3max\left\{|{\alpha }_1-({\gamma }_2+{\beta }_3)|,|{\alpha }_2-({\beta }_1+{\gamma }_3)|,|{\alpha }_3-({\gamma }_1+{\beta }_2)|\right\}\hfill \\ \hfill & \times {\int }_0^T\mathbb{E}V\left(x_1({\mu }_m\wedge t),x_2({\mu }_m\wedge t),x_3({\mu }_m\wedge t)\right)dt,\hfill \end{array}$$

(13)

where $L_1:=V\left(x_1\left(0\right),x_2\left(0\right),x_3\left(0\right)\right)+LT$. Using the Gronwall inequality, we obtain from Equation (13) that

$$\begin{array}{cc}\hfill & \mathbb{E}V\left(x_1({\mu }_m\wedge T),x_2({\mu }_m\wedge T),x_3({\mu }_m\wedge t)\right)\hfill \\ \hfill & \le L_1{e}^{3Tmax\left\{|{\alpha }_1-({\beta }_3+{\gamma }_2)|,|{\alpha }_2-({\beta }_1+{\gamma }_3)|,|{\alpha }_3-({\beta }_2+{\gamma }_1)|\right\}}.\hfill \end{array}$$

(14)

Since for each $\eta \in \{{\mu }_m\wedge T\}$ there certainly exists one of $x_1({\mu }_m,\eta )$ or $x_2({\mu }_m,\eta )$ or $x_3({\mu }_m,\eta )$, which are equal to m or $1/m$, therefore, $V\left(x_1({\mu }_m\wedge T),x_2({\mu }_m\wedge T),x_3({\mu }_m\wedge t)\right)\ge (m-1-lnm)\wedge \left(\frac{1}{m}+lnm-1\right)$. Then it follows from Equation (14) that

$$\begin{array}{cc}\hfill & L_1{e}^{3Tmax\left\{|{\alpha }_1-({\gamma }_2+{\beta }_3)|,|{\alpha }_2-({\beta }_1+{\gamma }_3)|,|{\alpha }_3-({\gamma }_1+{\beta }_2)|\right\}}\hfill \\ \hfill & \ge \mathbb{E}V\left(x_1({\mu }_m\wedge T),x_2({\mu }_m\wedge T),x_3({\mu }_m\wedge t)\right)\hfill \\ \hfill & \ge \mathbb{E}\left[I_{{\mu }_m\le T}\left(\eta \right)V\left(x_1({\mu }_m\wedge T),x_2({\mu }_m\wedge T),x_3({\mu }_m\wedge t)\right)\right]\hfill \\ \hfill & \ge P\left\{{\mu }_m\le T\right\}(m-lnm-1)\wedge (\frac{1}{m}+lnm-1),\hfill \end{array}$$

here, $I_{{\mu }_m\le T}$ should be the indicator function of $\left\{{\mu }_m\le T\right\}$. Since $m\to \infty$, there exists ${lim}_{m\to \infty }P\left\{{\mu }_m\le T\right\}=0$; therefore, $P\left\{{\mu }_{\infty }\le T\right\}=0$. Since T is an arbitrary positive, we must have $P\left\{{\mu }_{\infty }<\infty \right\}=0$. Therefore, $P\left\{{\mu }_{\infty }=\infty \right\}=1$ is the required result. □

Remark 1.

It is essential to the inspection whether or not the solution of Equation (1), along with initial values of Equation (2), will not collapse to infinity in a finite time (global existence). Indeed, we cannot obtain the global existence of the proposed solution only from the explicit expression of the given system. Although Lemma 3 is fundamental to the study of the global existence of the positive solution for the proposed system of Equation (1). It is worth declaring that by using Lemma 3, we can show the proposed stochastic delay Equation (1) in the sense that we have a positive solution that will not collapse to infinity in finite time.

5. Main Results

In the present section, we discuss the important properties of the proposed system given in Equation (1), which are the criteria for the alternate boundedness in the mean.

Theorem 1.

If $({\alpha }_1>{\beta }_3+{\gamma }_2)$, $({\alpha }_2>{\beta }_1+{\gamma }_3)$ and $({\alpha }_3>{\beta }_2+{\gamma }_1)$, then the global solution $x\left(t\right)=\left(x_1\left(t\right),x_2\left(t\right),x_3\left(t\right)\right)$ of Equation (1) with the initial values Equation (2) of $t\ge 0$ are positive almost surely and satisfy:

$$\begin{array}{c}\hfill \underset{t\to \infty }{lim}sup\mathbb{E}\left|x\left(t\right)\right|\le \frac{\delta }{\alpha }\end{array}$$

(15)

and

$$\begin{array}{c}\hfill \underset{t\to \infty }{lim}sup\frac{1}{t}{\int }_0^t\mathbb{E}\left(x_1^p\left(t\right)+x_2^p\left(t\right)+x_3^p\left(t\right)\right)ds\le Q_1+Q_2+Q_3,\end{array}$$

(16)

where $\alpha =min\left\{{\alpha }_1-({\beta }_3+{\gamma }_2),{\alpha }_2-({\beta }_1+{\gamma }_3),{\alpha }_3-({\beta }_2+{\gamma }_1)\right\}$, $\delta ={\sum }_{i=1}^3{\delta }_i{\left(\frac{{\nu }_i}{{\rho }_{i}e}\right)}^{{\nu }_i}$, and $p\ge 1$ such that $A_1:={\alpha }_1-({\beta }_3+{\gamma }_2)-\frac{p-1}{2}{\sigma }_1^2+\frac{p-1}{p}\left({\beta }_3+{\gamma }_2-({\beta }_1+{\gamma }_1)\right)>0$,$A_2:={\alpha }_2-({\beta }_1+{\gamma }_3)-\frac{p-1}{2}{\sigma }_2^2+\frac{p-1}{p}\left({\beta }_1+{\gamma }_3-({\beta }_2+{\gamma }_2)\right)>0$ and$A_3:={\alpha }_3-({\beta }_2+{\gamma }_1)-\frac{p-1}{2}{\sigma }_3^2+\frac{p-1}{p}\left({\beta }_2+{\gamma }_1-({\beta }_3+{\gamma }_3)\right)>0$,$Q_i={max}_{y\ge 0}\left\{-p{A}_i{y}^p+p{\delta }_i{\left(\frac{{\nu }_i}{{\rho }_{i}e}\right)}^{{\nu }_i}{y}^{p-1}\right\},i=1,2,3.$ Namely, Equation (1) is ultimately bounded in the mean.

Proof. 

In the highlights of Lemma 3, we can easily see that $x\left(t\right)>0$ for $t\ge 1$ almost surely. Moreover, by using Lemma 1, we get:

$$\begin{array}{c}\hfill d\left(x_1\left(t\right)+x_2\left(t\right)+x_3\left(t\right)\right)\le \left[-\alpha \left(x_1\left(t\right)+x_2\left(t\right)+x_3\left(t\right)\right)+\delta \right]dt+\sum _{i=1}^3{\sigma }_i{x}_i\left(t\right)d{B}_i\left(t\right),\end{array}$$

(17)

Now, applying It$\widehat{o}$ formula, Equation (17) takes the form:

$$\begin{array}{c}\hfill d\left[e^{\alpha t}\left(x_1\left(t\right)+x_2\left(t\right)+x_3\left(t\right)\right)\right]\le \delta e^{\alpha t}dt+\sum _{i=1}^3{\sigma }_i{e}^{\alpha t}{x}_i\left(t\right)d{B}_i\left(t\right),\end{array}$$

(18)

now integrating Equation (18) from $0,t$, we get:

$$\begin{array}{cc}\hfill & e^{\alpha t}\mathbb{E}\left(x_1\left(t\right)+x_2\left(t\right)+x_3\left(t\right)\right)\le x_1\left(0\right)+x_2\left(0\right)+x_3\left(0\right)+\frac{\delta }{\alpha }(e^{\alpha t}-1),\hfill \\ \hfill & \Rightarrow \underset{t\to \infty }{lim}sup\mathbb{E}\left(x_1\left(t\right)+x_2\left(t\right)+x_3\left(t\right)\right)\le \frac{\delta }{\alpha },\hfill \end{array}$$

In view of Lemma 1, Young’s inequality and the It$\widehat{o}$ formula follow from Equation (1), such that:

$$\begin{array}{cc}\hfill d\left(x_1^p\left(t\right)+x_2^p\left(t\right)+x_3^p\left(t\right)\right)& =p\{-\left({\alpha }_1-\frac{p-1}{2}{\sigma }_1^2\right)x_1^p\left(t\right)-\left({\alpha }_2-\frac{p-1}{2}{\sigma }_2^2\right)x_2^p\left(t\right)\hfill \\ \hfill & -\left({\alpha }_3-\frac{p-1}{2}{\sigma }_3^2\right)x_3^p\left(t\right)+({\beta }_3+{\gamma }_2)x_1^{p-1}\left(t\right)x_2\left(t\right)x_3\left(t\right)\hfill \\ \hfill & +({\beta }_1+{\gamma }_3)x_2^{p-1}\left(t\right)x_1\left(t\right)x_3\left(t\right)+({\beta }_2+{\gamma }_1)x_3^{p-1}\left(t\right)x_2\left(t\right)x_1\left(t\right)\hfill \\ \hfill & +\sum _{i=1}^3{\delta }_i{x}_i^{p-1}\left(t\right)x_i^{{\nu }_i}(t-{\tau }_i)e^{-{\rho }_i{x}_i(t-{\tau }_i}\}dt+\sum _{i=1}^{3}p{\sigma }_i{x}_i^p\left(t\right)d{B}_i\left(t\right)\hfill \\ \hfill & \le p\{-({\alpha }_1-({\beta }_3+{\gamma }_2)-\frac{p-1}{2}{\sigma }_1^2\hfill \\ \hfill & +\frac{p-1}{p}\left({\beta }_3+{\gamma }_2-({\beta }_1+{\gamma }_1)\right))x_1^p\left(t\right)\hfill \\ \hfill & -\left({\alpha }_2-({\beta }_1+{\gamma }_3)-\frac{p-1}{2}{\sigma }_2^2+\frac{p-1}{p}\left({\beta }_1+{\gamma }_3-({\beta }_2+{\gamma }_2)\right)\right)x_2^p\left(t\right)\hfill \\ \hfill & -\left({\alpha }_3-({\beta }_2+{\gamma }_1)-\frac{p-1}{2}{\sigma }_3^2+\frac{p-1}{p}\left({\beta }_2+{\gamma }_1-({\beta }_3+{\gamma }_3)\right)\right)x_3^p\left(t\right)\hfill \\ \hfill & +\sum _{i=1}^3{\delta }_i{\left(\frac{{\nu }_i}{{\rho }_{i}e}\right)}^{{\nu }_i}{x}_i^{p-1}\left(t\right)\}dt+\sum _{i=1}^{3}p{\sigma }_i{x}_i^p\left(t\right)d{B}_i\left(t\right)\hfill \\ \hfill & =\sum _{i=1}^3\left\{-p{A}_i{x}_i^p\left(t\right)+p{\delta }_i{\left(\frac{{\nu }_i}{{\rho }_{i}e}\right)}^{{\nu }_i}{x}_i^{p-1}\left(t\right)\right\}dt+\sum _{i=1}^{3}p{\sigma }_i{x}_i^p\left(t\right)d{B}_i\left(t\right)\hfill \\ \hfill & \le \sum _{i=1}^3{Q}_{i}dt+\sum _{i=1}^{3}p{\sigma }_i{x}_i^p\left(t\right)d{B}_i\left(t\right),\hfill \end{array}$$

which suggests

$$\begin{array}{c}\hfill \underset{t\to \infty }{lim}sup\frac{1}{t}{\int }_0^t\mathbb{E}\left\{x_1^p\left(t\right)+x_2^p\left(t\right)+x_3^p\left(t\right)\right\}ds\le Q_1+Q_2+Q_3.\end{array}$$

(19)

□

To define the asymptotic estimation for the solution of almost surely, Mao [39] defines the assumptions: ${lim}_{t\to \infty }sup\frac{1}{t}ln\left|x\left(t\right)\right|$, known as the sample Lyapunov exponent. Therefore, we will next estimate the Lyapunov exponent of Equation (1) along with with the initial conditions of Equation (2).

Theorem 2.

The sample Lyapunov exponent of the $x\left(t\right)=\left(x_1\left(t\right),x_2\left(t\right),x_3\left(t\right)\right)$ solution of Equation (1) with initial the conditions of Equation (2) satisfies:

$$\begin{array}{c}\hfill \underset{t\to \infty }{lim}sup\frac{lnx\left(t\right)}{t}\le \frac{G}{3}\end{array}$$

(20)

where $G=D\left\{-3{\alpha }_1+{\beta }_1+{\beta }_2+{\beta }_3+{\sigma }_1^2,-3{\alpha }_2+{\beta }_1+{\beta }_2+{\beta }_3+{\sigma }_2^2,-3{\alpha }_3+{\beta }_1+{\beta }_2+{\beta }_3+{\sigma }_3^2,-3{\alpha }_1+{\gamma }_1+{\gamma }_2+{\gamma }_3+{\sigma }_1^2,-3{\alpha }_2+{\gamma }_1+{\gamma }_2+{\gamma }_3+{\sigma }_2^2,-3{\alpha }_3+{\gamma }_1+{\gamma }_2+{\gamma }_3+{\sigma }_3^2,3{\delta }_1{\left(\frac{{\nu }_1}{{\rho }_{1}e}\right)}^{{\nu }_1},3{\delta }_2{\left(\frac{{\nu }_2}{{\rho }_{2}e}\right)}^{{\nu }_2},3{\delta }_3{\left(\frac{{\nu }_3}{{\rho }_{3}e}\right)}^{{\nu }_3}\right\}.$

Proof. 

Using the Young inequality, It$\widehat{o}$ formula and Lemma 1 along with Lemma 2, then from Equation (1) we get:

$$\begin{array}{cc}\hfill ln(1+x_1^2\left(t\right)& +x_2^2\left(t\right)+x_3^2\left(t\right))\hfill \\ \hfill & =ln\left(1+x_1^2\left(0\right)+x_2^2\left(0\right)+x_3^2\left(0\right)\right)+{\int }_0^t\frac{1}{1+x_1^2\left(s\right)+x_2^2\left(s\right)+x_3^2\left(s\right)}\hfill \\ \hfill & \times [(-3{\alpha }_1+{\sigma }_1^2)x_1^2\left(s\right)+(-3{\alpha }_2+{\sigma }_2^2)x_2^2\left(s\right)+(-3{\alpha }_3+{\sigma }_3^2)x_3^2\left(s\right)\hfill \\ \hfill & +3({\beta }_1+{\beta }_2+{\beta }_3+{\gamma }_1+{\gamma }_2+{\gamma }_3)x_1\left(s\right)x_2\left(s\right)x_3\left(s\right)\hfill \\ \hfill & +\sum _{i=1}^{3}3{\delta }_i{x}_i\left(s\right)x_i^{{\nu }_i}(s-{\tau }_i)e^{-{\rho }_i{x}_i(s-\tau )}]ds\hfill \\ \hfill & +\sum _{i=1}^3\left[M_i\left(t\right)-{\int }_0^t\frac{3{\sigma }_i^2{x}_i^4\left(s\right)}{{\left(1+x_1^2\left(s\right)+x_2^2\left(s\right)+x_3^2\left(s\right)\right)}^2}ds\right]\hfill \\ \hfill & \le ln\left(1+x_1^2\left(0\right)+x_2^2\left(0\right)+x_3^2\left(0\right)\right)+{\int }_0^t\frac{1}{1+x_1^2\left(s\right)+x_2^2\left(s\right)+x_3^2\left(s\right)}\hfill \\ \hfill & \times [(-3{\alpha }_1+{\beta }_1+{\gamma }_1+{\beta }_2+{\gamma }_2+{\beta }_3+{\gamma }_3+{\sigma }_1^2)x_1^2\left(s\right)\hfill \\ \hfill & +(-3{\alpha }_2+{\beta }_1+{\gamma }_1+{\beta }_2+{\gamma }_2+{\beta }_3+{\gamma }_3+{\sigma }_2^2)x_2^2\left(s\right)\hfill \\ \hfill & +(-3{\alpha }_3+{\beta }_1+{\gamma }_1+{\beta }_2+{\gamma }_2+{\beta }_3+{\gamma }_3+{\sigma }_3^2)x_3^2\left(s\right)\hfill \\ \hfill & +\sum _{i=1}^{3}3{\delta }_i{\left(\frac{{\nu }_i}{{\rho }_{i}e}\right)}^{{\nu }_i}{x}_i\left(s\right)]ds+\sum _{i=1}^3\left[M_i\left(t\right)-{\int }_0^t\frac{3{\sigma }_i^2{x}_i^4\left(s\right)}{{\left(1+x_1^2\left(s\right)+x_2^2\left(s\right)+x_3^2\left(s\right)\right)}^2}ds\right]\hfill \\ \hfill & \le ln\left(1+x_1^2\left(0\right)+x_2^2\left(0\right)+x_3^2\left(0\right)\right)+{\int }_0^{t}Gds\hfill \\ \hfill & +\sum _{i=1}^3\left[M_i\left(t\right)-{\int }_0^t\frac{3{\sigma }_i^2{x}_i^4\left(s\right)}{{\left(1+x_1^2\left(s\right)+x_2^2\left(s\right)+x_3^2\left(s\right)\right)}^2}ds\right]\hfill \end{array}$$

(21)

$M_i\left(t\right)={\int }_0^t\frac{3{\sigma }_i{x}_i^2\left(s\right)}{1+x_1^2\left(s\right)+x_2^2\left(s\right)+x_3^2\left(s\right)}d{B}_i\left(s\right),i=1,2,3.$ Now, for each positive n, applications of the exponential martingale inequality [39] yield to:

$$\begin{array}{c}\hfill p\left\{\underset{0\le t\le n}{sup}\left[M_i\left(t\right)-{\int }_0^t\frac{3{\sigma }_i^2{x}_i^4\left(s\right)}{{\left(1+x_1^2\left(s\right)+x_2^2\left(s\right)+x_3^2\left(s\right)\right)}^2}ds\right]>3lnn\right\}\le \frac{1}{n^2},i=1,2,3.\end{array}$$

By applying the lemma of Borel-Cantelli, for certainly all $\omega \in \mathsf{\Lambda }$ there are $n_i=n_i\left(\omega \right)\ge 1(i=1,2,3)$ random integers such as:

$$\begin{array}{c}\hfill \underset{0\le t\le n}{sup}\left[M_i\left(t\right)-{\int }_0^t\frac{3{\sigma }_i^2{x}_i^4\left(s\right)}{{\left(1+x_1^2\left(s\right)+x_2^2\left(s\right)+x_3^2\left(s\right)\right)}^2}ds\right]\le 3lnn,n\ge n_i.\end{array}$$

Therefore,

$$\begin{array}{c}\hfill M_i\left(t\right)\le {\int }_0^t\frac{3{\sigma }_i^2{x}_i^4\left(s\right)}{{\left(1+x_1^2\left(s\right)+x_2^2\left(s\right)+x_3^2\left(s\right)\right)}^2}ds+3lnn,(i=1,2,3).\end{array}$$

(22)

Then using Equation (21), together with Equation (22), implies

$$\begin{array}{c}\hfill ln\left(1+\sum _{i=1}^3{x}_i^2\left(t\right)\right)\le Gt+4lnn+ln\left(1+\sum _{i=1}^3{x}_i^2\left(0\right)\right),\end{array}$$

However, each $0\le t\le n,n\ge n_1\vee n_2\vee n_3$. Hence for each $\omega \in \mathsf{\Lambda }$, if $n\ge n_1\vee n_2\vee n_3,n-1\le t\le n$, certainly we have:

$$\begin{array}{c}\hfill \frac{ln\left(1+x_1^2\left(t\right)+x_2^2\left(t\right)+x_3^2\left(t\right)\right)}{t}\le \frac{\left[Gt+4lnn+ln\left(1+x_1^2\left(0\right)+x_2^2\left(0\right)+x_3^2\left(0\right)\right)\right]}{n-1}.\end{array}$$

When n tends to infinity, then we get:

$$\begin{array}{cc}\hfill \underset{n\to \infty }{lim}sup\frac{ln{x}_i\left(t\right)}{t}& \le \underset{n\to \infty }{lim}sup\frac{ln\left(1+x_1^2\left(t\right)+x_2^2\left(t\right)+x_3^2\left(t\right)\right)}{3t}\hfill \\ \hfill & \le \underset{n\to \infty }{lim}sup\frac{\left[Gn+\frac{4lnn}{n-1}+ln\left(1+x_1^2\left(0\right)+x_2^2\left(0\right)+x_3^2\left(0\right)\right)\right]}{3(n-1)}\hfill \\ \hfill & =\frac{G}{3},i=1,2,3.\hfill \end{array}$$

(23)

□

Remark 2.

For the existence of a positive solution, the conditions are not necessary from Lemma 3. Therefore, in this article, we have generalized the main results [29,32,40].

6. Results and Discussion

In the present section, we provide some test examples along with numerical simulations to confirm the theoretical justifications.

Consider the stochastic delay differential NGM system given in Equation (1), with the parameter values given by ${\alpha }_1=1.32,{\alpha }_2=1.9,{\alpha }_3=1.9,{\beta }_1=1,{\beta }_2=1,{\beta }_3=0.5,{\gamma }_1=1,{\gamma }_2=0.8,{\gamma }_3=0.6,{\delta }_1=3,{\delta }_2=2,{\delta }_3=2,{\nu }_i=2,{\tau }_i={\rho }_i=1(i=1,2,3)$, with initial values ${\varrho }_1={\varrho }_2={\varrho }_3=1$. From Theorem 2 with Lemma 3, it follows that the proposed three connected stochastic delay neoclassical growth systems, along with the initial conditions given in Equation (2), have a unique global positive solution, as shown in Figure 1. It also satisfies the sample Lyapunov exponent for the proposed parameter values ${lim}_{n\to \infty }sup\frac{1}{t}ln{x}_i\left(t\right)\le 24/e^2,(i=1,2,3).$ Although we choose $p=1.5$ then we have each $A_1:={\alpha }_1-({\beta }_3+{\gamma }_2)-\frac{p-1}{2}{\sigma }_1^2+\frac{p-1}{p}\left({\beta }_3+{\gamma }_2-({\beta }_1+{\gamma }_1)\right)>0$, $A_2:={\alpha }_2-({\beta }_1+{\gamma }_3)-\frac{p-1}{2}{\sigma }_2^2+\frac{p-1}{p}\left({\beta }_1+{\gamma }_3-({\beta }_2+{\gamma }_2)\right)>0$ and $A_3:={\alpha }_3-({\beta }_2+{\gamma }_1)-\frac{p-1}{2}{\sigma }_3^2+\frac{p-1}{p}\left({\beta }_2+{\gamma }_1-({\beta }_3+{\gamma }_3)\right)>0$, $Q_i={max}_{y\ge 0}\left\{-p{A}_i{y}^p+p{\delta }_i{\left(\frac{{\nu }_i}{{\rho }_{i}e}\right)}^{{\nu }_i}{y}^{p-1}\right\},i=1,2,3.$ Namely, Equation (1) is ultimately bounded in the mean by Theorem 1, as shown in Figure 1. Similarly, for the same parameter values given in Figure 1, we draw the comparison of the deterministic (to take in Equation (1) ${\sigma }_i=0,i=1,2,3$) with the stochastic one in Figure 2. We can clearly see that both solutions are in very good agreement. In Figure 3, we use the parameter values ${\alpha }_1=1.26,{\alpha }_2=1.8,{\alpha }_3=1.6,{\beta }_1=1,{\beta }_2=0.8,{\beta }_3=0.5,{\gamma }_1=1,{\gamma }_2=0.8,{\gamma }_3=0.7,{\delta }_1=3,{\delta }_2=2,{\delta }_3=2,{\nu }_i=2,{\tau }_i={\rho }_i=1,{\sigma }_i=1,(i=1,2,3)$. For the above parameter values, the proposed stochastic delay NGM system has an unstable positive solution, clearly seen in Figure 3. Again, using the same parameter values as given in Figure 3 above, we draw the comparisons of both the stochastic system with the deterministic one in Figure 4. Using the parameter values ${\alpha }_1=1.3,{\alpha }_2=1.9,{\alpha }_3=1.9,{\beta }_1=1,{\beta }_2=1,{\beta }_3=0.5,{\gamma }_1=1,{\gamma }_2=0.8,{\gamma }_3=0.6,{\delta }_1=5,{\delta }_2=5,{\delta }_3=4,{\nu }_i=0.7,{\tau }_i={\rho }_i=1,{\sigma }_i=1,(i=1,2,3)$. From the above parameter values, the system given in Equation (1) satisfies the sample Lyapunov exponent and each $A_i,i=1,2,3$ is greater then zero, along with $Q_i={max}_{y\ge 0}\left\{-p{A}_i{y}^p+p{\delta }_i{\left(\frac{{\nu }_i}{{\rho }_{i}e}\right)}^{{\nu }_i}{y}^{p-1}\right\},i=1,2,3.$ Therefore, from Theorem 1 and Theorem 2, the models are exponentially mean square stable and merge to zero, as shown in Figure 5. Similarly, for the same parameter values as given in Figure 6 above, we draw the comparison of the deterministic with the stochastic one.

7. Conclusions

In this article, we consider a novel approach for three connected delay differential NGMs under stochastic perturbations. It is observed that the nonlinearity and delay can be sources of continuous time chaos. Constant delay can generate complex dynamics involving chaos via period-doubling bifurcation. Stability conditions for positive and zero equilibria of the proposed model are obtained. For numerical simulations, we convert the proposed system to a nonlinear system using a polynomial with Legendre-Gauss quadrature and respective weight functions. We consider both deterministic and stochastic models. It is shown that the proposed stochastic delay NGM system given in Equation (1), along with initial conditions given in Equation (2), has a global positive solution that is conclusively bounded. The numerical results confirm the theoretical justifications.