Exploring the solutions of a financial bubble model via a new fractional derivative

1.   Introduction

Previous studies indicate that categorizing and analyzing the financial system's price crashes, bubbles, momentum, and liquidity have played a vital role in financial behavior over the last few years. Asset bubbles exist when the prices of commodities such as real estate, stocks, or gold shoot up abruptly within a short period, not necessarily due to their intrinsic value. Consulting Evans ^[1] provides additional understanding of important issues relating to new financial instruments. This economic model helped us to understand financial fundamentals, including the formation of a financial bubble and how to deal with it. One very good indicator of a bubble is when the market value is way over and above the intrinsic value of an asset, whereby people seemingly buy a particular asset for no good reason, known as irrational exuberance. In other words, a bubble occurs when market prices are greater than the asset's fundamental value.

For several years, great effort has been devoted to studying the financial bubbles. According to Barlevy ^[2], a financial bubble occurs when the market price of an asset is higher than its fair or main value. This deviation from true value can lead to concerns about market distortions. Furthermore, financial bubbles can develop even when an asset's essential value does not change. Thus, while trends in fashion may show price fluctuations that correspond to shifts in fair value, sudden price spikes in other assets may indicate a bubble. Conventional macroeconomic models, which assume perfect financial markets and ignore financial frictions, prove inadequate for understanding financial crises ^[3,4,5,6]. These models assume that financial frictions only affect non-financial firms and view financial intermediaries as a simple curtain. Nonetheless, groundbreaking research in the literature, such as that carried out by ^[7,8], has examined different types of systems and shed light on the intricate dynamics of financial crises. As demand for the asset increases, herd behavior causes prices to soar. This disrupts the normal supply and demand equilibrium and leads to market instability. During the growth phase, the bubble often experiences exponential expansion beyond the natural increase in market size. Anyhow, this unsustainable price direction eventually collapses, especially when triggered by even minor events. This has made the use of the two-scale economic theory, introduced by He ^[9,10], more effective in studying such bubbles, which will give more accurate and reliable results. In summary, an asset bubble occurs when market prices become disconnected from the fundamental value, and understanding the causes and signs can help protect investors when these bubbles burst. According to Sornette and Cauwels ^[11], a bubble usually begins when new information, such as the opening of a new market, captures the market's attention and raises expectations of future performance. Investors are drawn in as word spreads of the potential for high returns. First, the astute investors spot the real opportunity, then the inattentive ones.

The Riemann-Liouville and Caputo operators are examples of nonlocal fractional operators that contain a singular kernel that occasionally fail to explain complex dynamic systems. As a result, researchers have offered several options to improve the description of actual event models by using a novel strategy and an additional tool. In this regard, novel fractional operators including a nonsingular kernel ^[12] have emerged. In fact, researchers have made significant progress in applying fractional calculus to real-world problems. The Atangana-Baleanu-Caputo(ABC) operator is the most optimal emulation operator among nonsingular kernel operators, and it depends on the Mittag-Leffler function (MLF) ^[13]. Then, Abdeljawad and Baleanu ^[14], created a new generalized Atangana-Baleanu -Caputo (GABC) operator of singular and nonsingular kernels with a generalized MLF of three parameters. Recently, a generalized weighted $\phi$-fractional operator covering all definitions of nonsingular kernels have been presented by Thabet et al. ^[15]. These new fractional derivative operators have applications in science, engineering, and financial bubbles ^{[16,17,18,19]}.

Using fixed point theory, several scholars have investigated the qualitative properties of solutions to fractional problems ^{[20,21,22,23]}. In particular, Wang ^[24] employed the Guo-Krasnoselskii and Avery-Henderson fixed point theory. In addition, there are numerous simple ways to prove the stability of linear systems of fractional order. In ^[25,26,27] established the existence and Hyers–Ulam stability of solutions for their proposed equations. However, these approaches are not applicable to fractional-order nonlinear systems. Diethelm ^[28] showed that under certain conditions a fractional system is stable; however, this conclusion is only applicable to scalar fractional systems. Therefore, alternative methods must be used to establish the stability of nonlinear fractional systems in the vector situation. Li et al. ^[29] introduced a fractional-order extension of the direct Lyapunov approach as one such method. However, because the fractional-order scenario involves additional complexity in identifying a Lyapunov candidate function, the application of this technique is often very challenging. In order to demonstrate the stability of fractional systems, some authors have suggested the use of Lyapunov functionals. Although two well-known works^[30,31] can be consulted, there is no clear relationship between the fractional differential equation and the Lyapunov function. In addition, ^[32] suggests more Lyapunovs where their relationship to the fractional differential problem is more fundamental; however, these functionals are neither simple nor restricted to certain types of fractional systems.

The fractional Adams–Bashforth method ^[33] is a powerful and excellent numerical approach that can produce a numerical solution that is closer to the precise answer; hence it is used to evaluate the approximate solution. This technique was created with classical differentiation, which takes the difference between two times, such as $\vartheta_{n+1}$ and $\vartheta_{n}$, and applies it to the fundamental theorem of calculus. We can observe, in the numerical simulation section, that the final formula of this technique comprises the fractional parameters, increasing the options for choosing that parameter and so providing more data about the dynamics behavior. Moreover, we show that the chaotic dynamic arising from the bubble's rupture may be effectively controlled.

Unlike previous models in the field of financial bubbles, this research paper contributes to studying the impact of behavioral contagion more accurately and more flexibly in the market by presenting a fractional-order mathematical model based on the generalized Mittag-Leffler function, providing a more comprehensive approach. Studying stability leads to a deeper understanding and opens new horizons and visions for those working and interested in the market, and provides numerical solutions to help in understanding how financial bubbles develop and how to address them.

1.1.   Model formulation

Very recently, the authors of ^[34] studied the following five-dimensional financial bubble model:

subject to the initial conditions:

where $\tilde{\mathcal{C}}$ refers to the population size within the economy, $\tilde{\mathcal{P}}$ is the asset price, $\tilde{\mathcal{V}}$ is the neutral sub-population, $\tilde{\mathcal{A}}$ is the optimist/bull group, $\tilde{\mathcal{U}}$ is the pessimist/bear group, and $\tilde{\mathcal{Q}}$ is refers actors before leaving the market permanently and becoming a quitter, and further parameter details are in Table 1.

Furthermore, most research examining financial bubble systems is limited to solving ordinary differential equations of integer order. It has been demonstrated that mathematical models that use integer-order ordinary differential equations are useful for comprehending bubbles dynamics. Motivated by the above work ^[34], in this article, we extend the abovementioned model (1.1) to the fractional differential model by using a new generalized fractional derivative in the Caputo sense as follows:

subject to the initial conditions:

where $^{c}D_{0, {\rm{ \mathsf{ σ}}} }^{{{\rm{ \mathsf{ α}}}, {\rm{ \mathsf{ β}}}, {\rm{ \mathsf{ μ}}}}}$ is a new generalized fractional derivative in the Caputo sense of order ${\rm{ \mathsf{ α}}} \in (0, 1), Re({\rm{ \mathsf{ μ}}}) > 0, {\rm{ \mathsf{ β}}} > 0,$ and ${\rm{ \mathsf{ σ}}} \in \mathbb{R},$ which will defined in Section 2.

It is worth declaring that the novelty and contributions of this study are the following:

(1) Exploring the financial bubble systems under a new generalized fractional derivative involving a modified MLF of three parameters (1.3)–(1.4).

(2) Investigating the E & U of the solution of the fractional financial bubble model (1.3)–(1.4) by utilizing the Banach and Schauder FPTs.

(3) Discovering the equilibrium point and identifying the nullcline points of the suggested model (1.3)–(1.4).

(4) Investigating the global stability of the discovered equilibrium point by using the Lyapunov function.

(5) Studying approximate solutions of the model (1.3)–(1.4) by utilizing the fractional Adams-Bashforth technique of order 3. We also compute the residual error to evaluate the correctness of the proposed method.

(6) Performing simulations with different parameter values and fractional orders to show the applicability of the method in different contexts.

(7) Our fractional model (1.3)–(1.4) returns to the framework of the GABC, ABC, and Caputo–Fabrizo-Caputo(CFC) operators for $({\rm{ \mathsf{ α}}} = {\rm{ \mathsf{ β}}}), ({\rm{ \mathsf{ α}}} = {\rm{ \mathsf{ β}}}, {\rm{ \mathsf{ σ}}} = {\rm{ \mathsf{ μ}}} = 1)$, and $({\rm{ \mathsf{ β}}} = {\rm{ \mathsf{ σ}}} = {\rm{ \mathsf{ μ}}} = 1)$, respectively.

This article is organized as follows: Several preliminary outcomes are presented in Section 2. Then, we investigate the E&U theorems in Section 3. A stability analysis of the model is performed in Section 4. Also, we determine a numerical solution of the above model by using the Adams-Bashforth technique in Section 5. In addition, we discuss an error in Section 6. Finally, we establish the numerical simulations of our results in Section 7.

2.   Preliminaries

In 2023, Thabet et al. ^[15] introduced new generalized fractional operator with respect to a function $\phi$ involving the weighted function $\omega$. One of its special cases represents a generalized to ABC and CFC operators for $\phi(\vartheta) = \vartheta,$ $\omega = 1$, and $\gamma = {\rm{ \mathsf{ β}}}$ as follows:

Definition 2.1. ^[15] The generalized fractional derivative in the framework of Caputo is defined as

where ${\rm{ \mathsf{ α}}}\in (0, 1), \, Re({\rm{ \mathsf{ μ}}}) > 0, \, {\rm{ \mathsf{ β}}} > 0, {\rm{ \mathsf{ σ}}}\in \mathbb{R}$, and a function ${f}\in H^{1}({a}, {b}).$ Furthermore, $E^{{\rm{ \mathsf{ σ}}}}_{{\rm{ \mathsf{ α}}}, {\rm{ \mathsf{ μ}}}}({\rm{ \mathsf{ λ}}}, z)$ is the modified Mittag-Leffler function (MLF) of three parameters defined as $E^{{\rm{ \mathsf{ σ}}}}_{{\rm{ \mathsf{ α}}}, {\rm{ \mathsf{ μ}}}}({\rm{ \mathsf{ λ}}}, z) = \sum\limits_{j = 0}^{\infty}{\rm{ \mathsf{ λ}}}^{j}\frac{({\rm{ \mathsf{ σ}}})_{j}\; z^{j{\rm{ \mathsf{ α}}}+{\rm{ \mathsf{ μ}}}-1}}{j!\; \Gamma(j{\rm{ \mathsf{ α}}}+{\rm{ \mathsf{ μ}}})}$, $\Gamma(\cdot)$ is the gamma function, ${\rm{ \mathsf{ λ}}} _{{\rm{ \mathsf{ α}}}} = \frac{-{\rm{ \mathsf{ α}}}}{1-{\rm{ \mathsf{ α}}}}$, $\Lambda (\cdot)$ is the normalization function such that $\Lambda (0) = \Lambda (1) = 1,$ and $({\rm{ \mathsf{ σ}}})_j = {\rm{ \mathsf{ σ}}}({\rm{ \mathsf{ σ}}}+1)\dots ({\rm{ \mathsf{ σ}}}+j-1).$

Remark 2.1. Be aware that the following established definitions can be found in the literature for every given value of the parameters ${\rm{ \mathsf{ α}}}, {\rm{ \mathsf{ β}}}, {\rm{ \mathsf{ μ}}},$ and ${\rm{ \mathsf{ σ}}}$.

(1) If ${\rm{ \mathsf{ α}}} = {\rm{ \mathsf{ β}}}$, in the formula (2.1), one finds the GABC fractional derivative that was recently defined in ^[35].

(2) If ${\rm{ \mathsf{ α}}} = {\rm{ \mathsf{ β}}},$ and ${\rm{ \mathsf{ σ}}} = {\rm{ \mathsf{ μ}}} = 1$, in the formula (2.1), we find the ABC fractional derivative that was defined in ^[13].

(3) If ${\rm{ \mathsf{ β}}} = {\rm{ \mathsf{ σ}}} = {\rm{ \mathsf{ μ}}} = 1$, in the formula (2.1), we get the CFC fractional derivative ^[36].

(4) The formula (2.1) reduces to the first derivative for ${\rm{ \mathsf{ α}}}, {\rm{ \mathsf{ σ}}}, {\rm{ \mathsf{ β}}}, {\rm{ \mathsf{ μ}}}\to 1$. Furthermore, it has a singular kernel for ${\rm{ \mathsf{ μ}}} \in(0, 1)$.

Moreover, the fractional derivative (2.1) is associated with the following fractional integral definition:

where $^{R}I_{a}^{i{\rm{ \mathsf{ β}}} -{\rm{ \mathsf{ μ}}}+1}{f}(\vartheta)$ is the Riemann-Liouville fractional integral defined in ^[37].

Lemma 2.1. ^[15] Let us consider ${\rm{ \mathsf{ α}}} \in (0, 1), Re({\rm{ \mathsf{ μ}}}) > 0, {\rm{ \mathsf{ β}}} > 0,$ and ${\rm{ \mathsf{ σ}}} \in \mathbb{R}$. In this case

Theorem 2.1. (^[38], Banach's FPT) Assume that ${\digamma}$ is a Banach space, and let ${\varkappa} :{G} \to{G}$ be a contraction operator, in which case; ${\varkappa}$ owns an exactly one fixed point in ${G}$.

Theorem 2.2. (^[39], Schauder's FPT) Assume that ${\digamma}$ is a Banach space and let ${G}\in {\digamma}$ be a convex, bounded, and closed set. If ${\varkappa} :{G} \to{G}$ is a continuous operator such that ${\varkappa}{G} \in {\digamma}$ and ${\varkappa} {G}$ is relatively compact, then ${\varkappa}$ possesses at least one fixed point in ${G}$.

3.   Qualitative theorems

This part studies the E & U properties for the model (1.3)–(1.4). This model can be expressed in the following form:

where

and we can write the model (1.3)–(1.4) as follows:

such that

Theorem 3.1. Let ${\rm{ \mathsf{ α}}} \in (0, 1), {\rm{ \mathsf{ β}}} > 0, Re({\rm{ \mathsf{ μ}}}) > 0,$ and ${\rm{ \mathsf{ σ}}} \in \mathbb{R}$. A function ${\rm{ \mathsf{ χ}}}$ is a solution of the model (3.3) if and only if (iff) ${\rm{ \mathsf{ χ}}}$ satisfies the following integral equation:

provided that ${\rm{ \mathsf{ μ}}} = 1$ requires ${f}({0}, {\rm{ \mathsf{ χ}}}(0)) = 0.$

Proof. By using Definition 2.2, and Lemma 2.1, we obtain

Conversely, if we take the derivative defined in Definition 2.1 on both sides of Eq (3.4), and apply Lemma 2.1, we get the model (3.3). □

Now, we consider the Banach space of the continuous function $\varXi = {C}[0, T]$ under the norm $||{\rm{ \mathsf{ χ}}}|| = \sup_{t\in [0, T]}|{\rm{ \mathsf{ χ}}}(\vartheta)|.$ Moreover, we define the Banach space $\Sigma = (\varXi^5, \|{\rm{ \mathsf{ χ}}}\|)$ with the norm

Next, we use Schauder's FPT to demonstrate the E & U of the solution of model (3.3).

Theorem 3.2. Let a function ${f}\in\Sigma$ be continuous, and there is a constant $\kappa_1 > 0, \ni |{f}(\vartheta, {\rm{ \mathsf{ χ}}}(\vartheta))| \leq \kappa_1 (1+|{\rm{ \mathsf{ χ}}}(\vartheta)|)$, and $\forall \vartheta \in [0, {T}]$ for each ${\rm{ \mathsf{ χ}}} \in \Sigma$. Then at least one solution for the model (3.3) exists, provided that

and ${\rm{ \mathsf{ μ}}} = 1$ requires ${f}({0}, {\rm{ \mathsf{ χ}}}(0)) = 0.$

Proof. The solution of model (3.3) is an analogous to the solution of the fractional integral equation (3.4). Let us define the operator $\mathcal{J}:\Sigma \to\Sigma$ as follows:

We now take the bounded closed convex ball defined as $B_{\varsigma} = \left\lbrace {\rm{ \mathsf{ χ}}} \in \Sigma : ||{\rm{ \mathsf{ χ}}}|| \leq \varsigma, \varsigma > 0\right\rbrace$, where $\varsigma \geq \dfrac{\Psi_{2}}{1-\Psi_{1}}$, such that $\Psi_{2} = |{\rm{ \mathsf{ χ}}}_0|+\Big(\sum\limits_{i = 0}^{\infty}\binom{{\rm{ \mathsf{ σ}}}}{i}\dfrac{(1-{\rm{ \mathsf{ α}}})^{1-i}{\rm{ \mathsf{ α}}}^{i}}{|\Lambda({\rm{ \mathsf{ α}}})|} \dfrac{T^{i{\rm{ \mathsf{ β}}}-{\rm{ \mathsf{ μ}}}+1}}{\Gamma(i{\rm{ \mathsf{ β}}}-{\rm{ \mathsf{ μ}}}+2)}\kappa_1 \Big).$

Firstly, we show that $(\mathcal{J}B_{\varsigma})\subset B_{\varsigma}.$ Thus $\forall t \in [0, T],$ and one has

For ${\rm{ \mathsf{ χ}}}\in B_{\varsigma}$, we get

This means that $(\mathcal{J}B_{\varsigma})\subset B_{\varsigma}$.

Now, we prove that the operator $\mathcal{J}$ is continuous. For this, we take the sequence $\left\lbrace {\rm{ \mathsf{ χ}}}_{n}\right\rbrace$ such that ${\rm{ \mathsf{ χ}}}_{n}\to{\rm{ \mathsf{ χ}}}$ in $B_{\varsigma}$ as $n\to\infty$. Then, for each $\vartheta\in[0, T]$, we find

By the continuity of a function ${f}$, we get $\|\mathcal{J}{\rm{ \mathsf{ χ}}}_{n}-\mathcal{J} {\rm{ \mathsf{ χ}}}\|\to 0$ as $n\to\infty$. Thus, the operator $\mathcal{J}$ is continuous.

Next, let us prove that $(\mathcal{J}B_{\varsigma})$ is a relatively compact operator. On account of the fact that $(\mathcal{J}B_{\varsigma})\subset B_{\varsigma}$, $(\mathcal{J}B_{\varsigma})$ is uniformly bounded.

For indicating $\mathcal{J}$ is equicontinuous on $B_{\varsigma}$, let ${\rm{ \mathsf{ χ}}}\in B_{\varsigma}$ and $\vartheta_{1}, \vartheta_{2}\in [0, T],$ and $\vartheta_{2} < \vartheta_{1}$. We obtain

We find that $({\vartheta_{1}^{i{\rm{ \mathsf{ β}}}-{\rm{ \mathsf{ μ}}}+1}-\vartheta_{2}^{i{\rm{ \mathsf{ β}}}-{\rm{ \mathsf{ μ}}}+1}})\to 0,$ as $\vartheta_{1}\to \vartheta_{2}.$ Hence, we deduce that $\mathcal{J}$ is equicontinuous. From the conclusions above, we infer that $\mathcal{J}$ is completely continuous. As a consequence of Theorem 2.2, we came to conclusion that the model (3.3) possesses at least one solution. □

In what follows, we study the uniqueness result of the solution of the model (3.3) by the Banach FPT.

Theorem 3.3. The model (3.3) possesses a unique solution if there is a positive real number $\kappa_2 > 0$, such that $\forall {\rm{ \mathsf{ χ}}}, {\bar{\rm{ \mathsf{ χ}}}}\in \Sigma.$ One has $|{f}(\vartheta, {\rm{ \mathsf{ χ}}}(\vartheta))-{{f}(\vartheta, {\bar{\rm{ \mathsf{ χ}}}}(\vartheta))}| \leq \kappa_2\, |{\rm{ \mathsf{ χ}}}(\vartheta)-{\bar{\rm{ \mathsf{ χ}}}}(\vartheta)|$, provided that

and ${\rm{ \mathsf{ μ}}} = 1$ requires ${f}({0}, {\rm{ \mathsf{ χ}}}(0)) = 0.$

Proof. Let us define the operator $\mathcal{J}:\Sigma \to\Sigma$ as given in (3.7). Then, for each ${\rm{ \mathsf{ χ}}}, {\bar{\rm{ \mathsf{ χ}}}}\in \Sigma,$ and $\vartheta\in[0, {T}]$, one finds

Hence,

Since $\Psi_{3} < 1$, $\mathcal{J}$ is a contraction operator. Based on Theorem 2.1, we conclude that the model (3.3) has a unique solution. □

4.   Stability analysis

This section aims to investigate the equilibrium point, nullclines, and the Lyapunov stability of the model (1.3)–(1.4) at an equilibrium point.

4.1.   Equilibrium point

The equilibrium point's are the point where each state variable does not change in value, marked by the derivative of each variable equal to zero. In order to compute it for the system (1.3), let us set

Thus,

Therefore, the equilibrium point is

Remark 4.1. The first equation in the model (1.3) does not correspond to the social contagion component's dynamics, as described by the final four equations. Thus, we must first determine the steady state $(\tilde{\mathcal{V}}, \tilde{\mathcal{A}}, \tilde{\mathcal{U}}, \tilde{\mathcal{Q}})$ to perform the stability analysis.

4.2.   Nullclines

The goal of this section is to identify the nullcline points from the system of Eq (1.3). Nullcline points contrast with equilibrium points as they are locations where the vectors are vertical, either ascending or descending. By definition, the ${\rm{ \mathsf{ χ}}}_{i}-nullcline$ is defined as

In our case, we find that a $\tilde{\mathcal{V}}-$nullcline point will be obtained with $\big(-\delta\tilde{\mathcal{V}}(\vartheta)\tilde{\mathcal{A}}(\vartheta)\big) = 0$; an $\tilde{\mathcal{A}}-$nullcline point will be determined by imposing with $\big(\delta\tilde{\mathcal{V}}(\vartheta)\tilde{\mathcal{A}}(\vartheta) - \big(\gamma \tilde{\mathcal{U}}(\vartheta) +\rho)\tilde{\mathcal{A}}(\vartheta)\big) = 0$; an $\tilde{\mathcal{U}}-$nullcline will be found due to $\big(\big(\gamma \tilde{\mathcal{U}}(\vartheta) +\rho)\tilde{\mathcal{A}}(\vartheta)-\frac{1}{\eta} \tilde{\mathcal{U}}(\vartheta)\big) = 0;$ and a $\mathcal Q-$nullcline will be obtained by $\big(\frac{1}{\eta} \tilde{\mathcal{U}}(\vartheta)\big) = 0$.

Thus, for the $\tilde{\mathcal{V}}-$nullcline point, we have the following set $(\tilde{\mathcal{V}}, 0, 0, 0), \; (0, \tilde{\mathcal{A}}, 0, 0), \;$and $\forall \tilde{\mathcal{V}}, \tilde{\mathcal{A}}\in \mathbb{R}^{+}$; for the $\tilde{\mathcal{A}}-$nullcline point, we have $(\tilde{\mathcal{V}}, 0, \tilde{\mathcal{U}}, 0), \;$ $(0, 0, \tilde{\mathcal{U}}, 0),$ $(0, \tilde{\mathcal{A}}, \frac{-\rho}{\gamma}, 0), \;$and $\forall \tilde{\mathcal{V}}, \tilde{\mathcal{U}}, \tilde{\mathcal{A}}\in \mathbb{R}^{+}$; for the $\tilde{\mathcal{U}}-$nullcline and $\tilde{\mathcal{Q}}-$nullcline, we have $(0, 0, 0, 0).$

4.3.   Lyapunov stability

Before we investigate the Lyapunov stability of the model (1.3)–(1.4) at the equilibrium point, we need to prove the following auxiliary result.

Theorem 4.1. Let ${\rm{ \mathsf{ χ}}}(\vartheta)$ be a derivable function and continuous. Then, at any instant in time $t\geq0$, for $0 < ({\rm{ \mathsf{ β}}} {i}+{\rm{ \mathsf{ μ}}}) < 1$, we have

Proof. The expression above is true if

By using Definition (2.1), it can be written as

Multiplying Eq (4.3), by ${\rm{ \mathsf{ χ}}}(\vartheta),$ we get

and

Therefore,

By subtracting Eq (4.5) from Eq (4.4), one finds

By putting $\mathcal{Y}(\mathtt{s}) = \big({\rm{ \mathsf{ χ}}}(\vartheta)-{\rm{ \mathsf{ χ}}}(\mathtt{s})\big)$, which yields ${\mathcal{Y}'}(\mathtt{s}) = -{{\rm{ \mathsf{ χ}}}'}(\mathtt{s})$ in (4.6), then we have

Integrating by part for Eq (4.7), implies

In that way (4.7), can written as:

Since there is an indeterminate at $\mathtt{s} = \vartheta$ in the first term of the expression (4.8), we examine the appropriate limit as follows:

Therefore, the expression (4.8) reduces to

since $0\leq \mathtt{s}\leq t$, and $0 < ({\rm{ \mathsf{ β}}} {{i}}+{\rm{ \mathsf{ μ}}}) < 1.$ This completes the proof. □

Theorem 4.2. The following inequality holds for $0 < {i}{\rm{ \mathsf{ β}}}+{\rm{ \mathsf{ μ}}} < 1$:

Then, the model (1.3)–(1.4) is globally stable at the equilibrium point.

Proof. First, we introduce the Lyapunov function as follows:

We note that $L(\tilde{\mathcal{V}}, \tilde{\mathcal{A}}, \tilde{\mathcal{U}}, \tilde{\mathcal{Q}}) > 0,$ for $(\tilde{\mathcal{V}}, \tilde{\mathcal{A}}, \tilde{\mathcal{U}}, \tilde{\mathcal{Q}})\neq (0, 0, 0, 0),$ and $L(\tilde{\mathcal{V}}, \tilde{\mathcal{A}}, \tilde{\mathcal{U}}, \tilde{\mathcal{Q}}) = 0$, at the equilibrium point $(\tilde{\mathcal{V}}^{0}, \tilde{\mathcal{A}}^{0}, \tilde{\mathcal{U}}^{0}, \tilde{\mathcal{Q}}^{0}) = (0, 0, 0, 0),$ then, by using Theorem 4.1, one has

Therefore, we deduce that ${}^{c}D_{0, {\rm{ \mathsf{ σ}}} }^{{{\rm{ \mathsf{ α}}}, {\rm{ \mathsf{ β}}}, {\rm{ \mathsf{ μ}}}}}L(\tilde{\mathcal{V}}, \tilde{\mathcal{A}}, \tilde{\mathcal{U}}, \tilde{\mathcal{Q}}) < 0,$ $\forall\, \, (\tilde{\mathcal{V}}, \tilde{\mathcal{A}}, \tilde{\mathcal{U}}, \tilde{\mathcal{Q}})\neq(0, 0, 0, 0)$, if the following condition holds:

Hence, ${}^{c}D_{0, {\rm{ \mathsf{ σ}}} }^{{{\rm{ \mathsf{ α}}}, {\rm{ \mathsf{ β}}}, {\rm{ \mathsf{ μ}}}}}L(\tilde{\mathcal{V}}, \tilde{\mathcal{A}}, \tilde{\mathcal{U}}, \tilde{\mathcal{Q}})$ is negative definitely. Therefore, according to (Lemma 5, ^[40]), we deduce that the origin of the fractional model (1.3)–(1.4) is asymptotically stable. □

5.   Numerical method

In this section, we introduce a numerical method which is used to solve the fractional model (1.3)–(1.4). This technique is the fractional Adams-Bashforth method of order 3 (FABM3). Thus, by using the Theorem 3.1, we have

Now, we provide a rough solution that entails approximation at the points $\vartheta = \vartheta_{n}$, and $\vartheta = \vartheta_{n+1},$ such that ${h} = \vartheta_{n+1}-\vartheta_{n}$, where ${h}$ is the size step, $n = 0, 1, \dots, p, \; \vartheta_0 = 0, \; and \; \vartheta_{p+1} = T.$ Therefore, if $\vartheta = \vartheta_{n}$, one has

and if $\vartheta = \vartheta_{n+1}$, we find

Subtracting Eq (5.1) from Eq (5.2), yield

Using the Lagrange interpolation and the interpolating functions ${F}_{n}$, ${F}_{n-1}$, ${F}_{n-2},$ the proposed technique is of order 3.

where $F_{n} = {f}(\vartheta_{n}, {\rm{ \mathsf{ χ}}}(\vartheta_{n})), F_{n-1} = {f}(\vartheta_{n-1}, {\rm{ \mathsf{ χ}}}(\vartheta_{n-1})), \; and \; F_{n-2} = {f}(\vartheta_{n-2}, {\rm{ \mathsf{ χ}}}(\vartheta_{n-2})).$

In what follows, we can rewrite Eq (5.3) based on the Newton polynomial as:

where $E$ is related to an error term and is defined by ^[41] as follows:

Next, let

Thus, by putting Eqs (5.4) and (5.6) into Eq (5.3), one gets

We can also write Eq (5.7), as follows:

where

and

The fractional integral (5.9) can be evaluated as:

Now, by choosing ${\rm{ \mathsf{ χ}}} = \big(\vartheta_{n+1}-\vartheta\big)$, which implies that $d{\rm{ \mathsf{ χ}}} = -d\vartheta,$ and then

As a result, the first fractional integral computation is provided by

The second fractional integral (5.10) can be computed in the same manner as the previous procedures. We found that

Substitution ${\rm{ \mathsf{ χ}}} = \big(\vartheta_{n}-\vartheta\big)$, which implies that $d{\rm{ \mathsf{ χ}}} = -d\vartheta$, and then

Hence,

Now, putting Eqs (5.11) and (5.12) into Eq (5.8) gives the following numerical strategy for the FABM3:

6.   Error analysis

The error is typically analyzed by comparing the numerical solution with the exact solution ^[42,43]. In this section, we locate the error analysis for the suggested numerical method that includes the generalized fractional derivative in the framework of Caputo. To calculate it, we go back to (5.5):

where

Therefore, the error term is

Now, by applying the absolute value on both sides for $s\in[\vartheta_{p}, \vartheta_{p+1}]$, we get

where ${\varUpsilon} = \sup\limits_{s\in [t_{p}, t_{p+1}]}|{f^{(3)}{(\mathtt{s}, {\rm{ \mathsf{ χ}}}(\mathtt{s})}}|$.

As a result, in what follows, we drive the error estimate for Eq (6.1), which is

7.   Numerical simulations

This part uses the earlier approximate solutions to present the numerical approach for all categories for the fractional model (1.3)–(1.4). The numerical solutions are plotted for numerous values of fractional orders. Here, we consider the following parameter values that are given in ^[34]: the primary constraints $\tilde{\mathcal{V}}(0) = 995, \tilde{\mathcal{A}}(0) = 5,$ and $\tilde{\mathcal{U}}(0) = \tilde{\mathcal{Q}}(0) = 0,$ the population size in economy $\tilde{\mathcal{C}} = 1000$; and the parameter values ${k}_1 = 2, \zeta = 0.0002, {v} = 0.0018$, ${d}_1^{*} = 0.02, {d}_2^{*} = 0.001, {\varepsilon} = 0.0002, \mathtt{s}_1^{*} = 0.03, \mathtt{s}_2^{*} = 0.004$, $\delta = 0.0001, \gamma = 0.0002, \rho = 0.0002,$ and $\eta = 30.$ We simulate our numerical results for various fractional-orders values in Figures 1–5 for a set of fractional order values $(\alpha, \beta, \mu, \sigma)$.

In Figure 1, the curves correspond to the situations one expects to observe during the formation of a financial bubble. We note that there is a volatile increase in the price of assets which would result from an increased demand for the asset, which is reflected in a higher price for the asset. This is accompanied by a rapid decrease in value due to the bubble's implosion. After the bubble's implosion, the asset's price seems to settle around its price just before the bubble breaks. The reason for this is said to be adjustment back to equilibrium after the huge variability created by financial bubbles. It can be noted that the fractional orders are directly related to the rate of increase as well as the rate of decrease in the price levels. In Figure 2, we observe a gradual decline in the number of individuals in the neutral group. This indicates that individuals tend to shift to other groups, either the optimists group or the pessimists group. This leads to the system reaching a state of equilibrium. Figure 3 can be divided into two phases. During the growth phase as prices rise, optimism increases as people believe the market will continue to rise, in the decline phase, after the peak is reached, excessive optimism gradually diminishes and the number of optimists begins to decline. This behavior is the typical pattern of financial bubbles. On the other hand, in Figure 4, we see a time lag in the appearance of pessimists, which indicates that the market cycle is sequential, with optimists appearing first, followed by pessimists. We see an increase in the number of pessimists, followed by a period of contraction and their transition to the group that exits the market. In Figure 5, we observe a gradual increase in the number of individuals who decide to exit the market, indicating a loss of confidence in the market.

On the other hand, the considered derivative is a more general operator and includes some already defined fractional derivatives like $GABC$, $ABC$, and $CFC$ as special cases. The numerical simulations can be performed easily over a long range. Figures 6–10 show the graphics of the components of the fractional model (1.3)–(1.4) for a set of fractional-order values $(\alpha, \beta, \mu, \sigma)$ via four types of fractional derivatives which are a new fractional derivative(FD) (2.1), and its special cases, $GABC$, $ABC$, and $CFC$. This proves that the newly used FD is more flexible and practical, and the results can be easily compared with some of the existing fractional derivatives in the literature.

8.   Conclusions

This paper provides a comprehensive, deeper, and more accurate understanding of a model related to financial bubbles and analyzes their stages (growth, collapse, and stabilization) using a new fractional-order derivative. The model identifies equilibrium points, and nullclines, and verifies the uniqueness and existence of the solution using the Schauder and Banach fixed-point theorems. Stability was studied using various theories, and a powerful technique namely the third-order Adams-Bashforth method, was applied to compute approximate solutions for the proposed model. Finally, a numerical simulation was conducted to confirm the accuracy of the employed methodology in handling various scenarios by using a step size of $h = 1$. The results show that the new fractional-order derivative is more accurate in capturing the dynamics of time changes compared with other derivatives. Additionally, our results can be applied to various fractional derivatives as special cases, such as the GABC, ABC, and CFC derivatives at certain values for the given parameters (see Figures 6–10). The results confirm that the technique can produce accurate answers in many settings.

Moreover, in the future, we aim to extract valuable information from massive financial data to improve the parameter estimation process with the greater flexibility of the model to align with varying market environments. It is possible to broaden the model by considering other elements that contribute to the development of financial bubbles. Combining these aspects will improve the model's realism as well as its ability to help understand and explain financial market phenomena and predict future events. This analysis can also be extended to other models in the future.

Authors contributions

Sabri T. M. Thabet: Conceptualization, data curation, formal analysis, methodology, software, writing-original draft, writing-review and editing; Reem M. Alraimy: Conceptualization, data curation, formal analysis, investigation, methodology, software, writing-original draft; Imed Kedim: Formal analysis, investigation, methodology, writing-review and editing; Aiman Mukheimer: Formal analysis, investigation, methodology, writing-review and editing; Thabet Abdeljawad: Formal analysis, investigation, methodology, writing-review and editing. All authors have read and agreed to the published version of the article.

Use of Generative-AI tools declaration

The authors declare they have not used artificial intelligence (AI) tools in the creation of this article.

Acknowledgments

This study is supported via funding from Prince Sattam bin Abdulaziz University project number (PSAU/2025/R/1446). The authors A. Mukheimer and T. Abdeljawad would like to thank Prince Sultan University for paying the APC and for the support through TAS research lab.

Conflict of interest

The authors declare that they have no competing interests.