Understanding Dynamics and Bifurcation Control Mechanism for a Fractional-Order Delayed Duopoly Game Model in Insurance Market

Abstract

Recently, the insurance industry in China has been greatly developed. The number of domestic insurance companies and foreign investment insurance companies has greatly increased. Competition between different insurance companies is becoming increasingly fierce. Grasping the internal competition law of different insurance companies is a very meaningful work. In this present work, we set up a novel fractional-order delayed duopoly game model in insurance market and discuss the dynamics including existence and uniqueness, non-negativeness, and boundedness of solution for the established fractional-order delayed duopoly game model in insurance market. By selecting the delay as a bifurcation parameter, we build a new delay-independent condition ensuring the stability and creation of Hopf bifurcation of the built fractional-order delayed duopoly game model. Making use of a suitable definite function, we explore the globally asymptotic stability of the involved fractional-order delayed duopoly game model. By virtue of hybrid controller which includes state feedback and parameter perturbation, we can effectively control the stability and the time of creation of Hopf bifurcation for the involved fractional-order delayed duopoly game model. The research indicates that time delay plays an all-important role in stabilizing the system and controlling the time of onset of Hopf bifurcation of the involved fractional-order delayed duopoly game model. To check the rationality of derived primary conclusions, Matlab simulation plots are explicitly presented. The established results in this manuscript are wholly novel and own immense theoretical guiding significance in managing and operating insurance companies.

1. Introduction

With the rapid development of the insurance market in China, different kinds of domestic insurance companies and foreign investment companies come forward in large numbers. In order to survive and better serve the people, the competition among various insurance companies is very fierce. The level of monopoly of the insurance market in China has gradually declined, but it still remains at the state of oligopoly. The monopoly competition game between oligarchs has become a very important research topic. In order to reveal the inherent law of competition among different oligarchs, it is necessary for us to set up mathematical models on competition and explore the quantitative relation of competition models among different oligarchs. Plenty of excellent and meaningful works on this topic have been published. For instance, Elabbsy et al. [1] set up a nonlinear triopoly game model concerning heterogeneous players and discussed the stability of equilibrium points, and displayed the bifurcation and chaos of this model via computer simulations. Yu and Yu [2] explored the stability and duality of dynamic Cournot and Bertrand duopoly model involving comprehensive preference. Tu and Wang [3] studied the local stability, Hopf bifurcation, and chaos control issue for a dynamic R&D Bertrand triopoly game system. Ma and Wang [4] made a detailed analysis on the complex dynamical behavior of a Cournot–Bertrand mixed game model concerning delayed bounded rationality. Tu et al. [5] dealt with the complex dynamics and chaos control for a Bertrand duopoly game model involving R&D activities. Zhang et al. [6] investigated the stability and bifurcation phenomenon for the Bertrand model concerning bounded rationality. For more concrete literature, we refer the readers to [7,8,9,10,11,12,13,14].

The price plays a very crucial role in competition among different insurance companies. The price is a very sensitive aspect, which has important effect on people’s needs and the distribution of products [7]. How to adjust the price of products is naturally a focus problem for insurance companies. During the past several decades, many valuable works on this theme have been reported. For example, Ma and Li [7] considered the complex dynamical behavior for a Bertrand–Stackelberg pricing model. Wei et al. [15] explored the pricing decision issue for complementary products involving companies’ different market powers. Mukhopadhyay et al. [16] discussed a Stackelberg price model of complementary goods under information asymmetry. Ma and Sun [17] studied the price game model and its complex characteristics of triopoly in different decision-making rule. In 2012, Xu and Ma [18] studied bifurcation dynamics of the following duopoly game involving time delay in insurance market:

$$\left\{\begin{array}{c}{\displaystyle {\dot{y}}_1\left(t\right)={\rho }_1{y}_1\left(t\right)[\alpha -2{\beta }_1\mu y_1\left(t\right)-2{\beta }_1(1-\mu )y_1(t-\theta )+{\delta }_1{y}_2\left(t\right)+{\beta }_1{\gamma }_1],}\hfill \\ {\displaystyle {\dot{y}}_2\left(t\right)={\rho }_2{y}_2\left(t\right)[\alpha -2{\beta }_2{y}_2\left(t\right)+{\delta }_2\mu y_1\left(t\right)+{\delta }_2(1-\mu )y_1(t-\theta )+{\beta }_2{\gamma }_2],}\hfill \end{array}\right.$$

(1)

where $y_1\left(t\right)$ stands for the price of the first insurance company, $y_2\left(t\right)$ stands for the price of the second insurance company, ${\rho }_1$ stands for the speed of price adjustment of the first insurance company, ${\rho }_2$ stands for the speed of price adjustment of the second insurance company, $\alpha$ represents the possible largest demand, ${\beta }_1$ denotes the effect of which the price of product 1 has on its quantity, ${\beta }_2$ denotes the effect of which the price of product 2 has on its quantity, ${\delta }_1$ stands for substitution rate which the products of 1 shows to the products of 2, ${\delta }_2$ stands for substitution rate which the products of 2 shows to the products of 1, $\mu \in (0,1)$ denotes the weight of the current price at time t, $1-\mu$ denotes the weight of the current price at time $t-\theta$, $\theta$ is a time delay, and all the parameters ${\rho }_i,{\beta }_i,{\gamma }_i,{\delta }_i,\alpha ,\mu ,\theta ,(i=1,2)$ are positive constants. For more details, one can refer to Refs. [18,19]. Taking advantage of stability criterion and bifurcation theory of delayed differential equation, Xu and Ma [18] set up a sufficient condition which guarantees the stability and the onset of Hopf bifurcation for model (1).

It is noteworthy that all the involved literature above on game model in insurance market (see [1,2,3,4,5,6,7,8,9,10,11,12,13,14,15,16,17,18,19]) are basically concerned with the integer-order dynamical models. A large number of studies indicate that fractional-order differential equation has been deemed as a more valid tactics to describe the authentic natural phenomenon in the world than the conventional integer-order counterparts. At present, fractional dynamical systems have been applied in many areas such ascomplex networks, biological systems, artificial intelligence, various waves in physics, viscoelasticity, capacitor principle, biomedical treatment, electrical engineering, economics, and so on [20,21,22,23,24,25,26]. Its great application value comes from the powerful memory trait and hereditary superiority for various materials and evolutionary process [27,28]. In recent years, fractional dynamical systems have attracted great attention from many scientific circles and great achievements have been acquired. For example, Xu et al. [29] revealed the impact of delay on Hopf bifurcation of a class of fractional-order delayed bidirectional associate memory neural networks. Eshaghi et al. [30] explored the Hopf bifurcation, chaos control, and synchronization issues for a chaotic fractional-order dynamical model. Zhou et al. [31] probed into the Hopf bifurcation control problem of fractional-order prey–predator system involving delays via hybrid controller. For more concrete publications, one can see [32,33,34,35,36].

Considering that the fractional-order delayed duopoly game model can better reflect the memory trait and hereditary superiority in price of two insurance companies and is motivated by the investigation above and based on model (1), in the current work, we will set up the following fractional-order delayed duopoly game model in insurance market:

$$\left\{\begin{array}{c}{\displaystyle {\mathcal{D}}^{\eta }{y}_1\left(t\right)={\rho }_1{y}_1\left(t\right)[\alpha -2{\beta }_1\mu y_1\left(t\right)-2{\beta }_1(1-\mu )y_1(t-\theta )+{\delta }_1{y}_2\left(t\right)+{\beta }_1{\gamma }_1],}\hfill \\ {\displaystyle {\mathcal{D}}^{\eta }{y}_1\left(t\right)={\rho }_2{y}_2\left(t\right)[\alpha -2{\beta }_2{y}_2\left(t\right)+{\delta }_2\mu y_1\left(t\right)+{\delta }_2(1-\mu )y_1(t-\theta )+{\beta }_2{\gamma }_2],}\hfill \end{array}\right.$$

(2)

where $\eta \in (0,1]$. All other parameters and variables own the same economic meaning as those in model (1).

In this manuscript, we principally probe into the following four problems: (a) Investigate the existence and uniqueness, non-negativeness, and boundedness of the solution for system (2); (b) Set up the delay-independent condition guaranteeing the stability and the occurrence of Hopf bifurcation of model (2); (c) Build the sufficient condition to ensure the globally asymptotic stability of model (2); and (d) Control the time of onset of Hopf bifurcation of model (2) via hybrid controller.

The chief contributions of this manuscript are elaborated as follows: (1) Based on the earlier works, a novel fractional-order delayed duopoly game model in insurance market is proposed. (2) The sufficient condition ensuring the globally asymptotic stability of model (2) is set up via constructing an appropriate positive definite function. (3) Hopf bifurcation of model (2) is successfully dominated via hybrid control strategy. So far, very few scholars focus on the Hopf bifurcation control issue of fractional-order models by utilizing hybrid controller. (4) The influence of time delay on the stability behavior and the occurrence of Hopf bifurcation of model (2) and its controlled system is revealed. (5) The research approach can be applied to study the bifurcation control issue of lots of fractional dynamical models in numerous areas.

The novelty of this research lies in the design of hybrid controller for the fractional-order delayed duopoly game model in insurance market. By designing a suitable hybrid controller, we can successfully control the stability region and Hopf bifurcation of model (2). The obtained results play a vital role in controlling the price of two insurance companies.

The structure of this research is arranged as follows. Some necessary basic knowledge about fractional-order differential equation is given in Section 2. Section 3 proves the existence and uniqueness, non-negativeness, and boundedness of the solution for model (2). A new delay-independent sufficient criterion which ensures the stability and the creation of Hopf bifurcation for model (2) is set up in Section 4. Section 5 explores the globally asymptotic stability of model (2) via a definite function. Hybrid control tactics are executed to control the stability and creation of Hopf bifurcation of model (2) in Section 6. Software simulation results are distinctly displayed to support the established key conclusions in Section 7. Section 8 draws a simple conclusion to complete this research.

2. Prerequisite Knowledge

In this segment, some essential basic theories about fractional-order differential equation are presented.

Definition 1

([37]). The fractional integral of order η for the function $h\left(ϵ\right)$ is defined in the following form:

$${\mathcal{I}}^{\eta }h\left(ϵ\right)=\frac{1}{\mathrm{\Gamma }\left(\eta \right)}{\int }_{{ϵ}_0}^{ϵ}{(ϵ-\nu )}^{\eta -1}h\left(\nu \right)d\nu ,$$

where $ϵ>{ϵ}_0,\eta >0,\mathrm{\Gamma }\left(\nu \right)={\int }_0^{\infty }{s}^{\nu -1}{e}^{-s}ds$ stands for the Gamma function.

Definition 2

([38]). Define the Caputo-type fractional-order derivative of order η for the function $h\left(ϵ\right)\in ([{ϵ}_0,\infty ),R)$ as follows:

$${\mathcal{D}}^{\eta }h\left(ϵ\right)=\frac{1}{\mathrm{\Gamma }(\iota -\eta )}{\int }_{{ϵ}_0}^{ϵ}\frac{h^{\left(l\right)}\left(s\right)}{{(ϵ-s)}^{\eta -l+1}}ds,$$

where $ϵ\ge {ϵ}_0$ and l denotes a positive integer ($l-1\le \eta <l)$. Peculiarly, if $\eta \in (0,1)$, then

$${\mathcal{D}}^{\eta }h\left(ϵ\right)=\frac{1}{\mathrm{\Gamma }(1-\eta )}{\int }_{{ϵ}_0}^{ϵ}\frac{h^{{}^{\prime }}\left(s\right)}{{(ϵ-s)}^{\eta }}ds.$$

Lemma 1.

([39]). Consider the fractional order system: ${\mathcal{D}}^{\eta }y=\mathcal{Q}y,y\left(0\right)=y_0$ where $\eta \in (0,1),$$y\in R^m,\mathcal{Q}\in R^{m\times m}.$ Let ${\nu }_l(l=1,2,\cdots ,m)$ denote the root of the characteristic equation of ${\mathcal{D}}^{\eta }y=\mathcal{Q}y$, then system ${\mathcal{D}}^{\eta }y=\mathcal{Q}y$ is said to be locally asymptotically stable ⇔$|arg\left({\nu }_l\right)|>\frac{\eta \pi }{2}(l=1,2,\cdots ,m)$. The system is stable ⇔$|arg\left({\nu }_l\right)|>\frac{\eta \pi }{2}(l=1,2,\cdots ,m)$ and each critical eigenvalue, which satisfies $|arg\left({\nu }_l\right)|=\frac{\eta \pi }{2}(l=1,2,\cdots ,m)$, owns geometric multiplicity one.

Lemma 2.

([39]). Let $\varphi \left(t\right)\in C[t_0,\infty )$ and satisfy

$$\left\{\begin{array}{c}{\displaystyle D^{\eta }\varphi \left(t\right)\le -{\sigma }_1\varphi \left(t\right)+{\sigma }_2,}\hfill \\ {\displaystyle \varphi \left(t_0\right)={\varphi }_{t_0},}\hfill \end{array}\right.$$

where $\eta \in (0,1),{\sigma }_1,{\sigma }_2\in R,{\sigma }_1\ne 0,t_0\ge 0$, then

$$\varphi \left(t\right)\le \left(\varphi \left(t_0\right)-\frac{{\sigma }_2}{{\sigma }_1}\right)E_{\eta }[-{\sigma }_1{(t-t_0)}^{\eta }]+\frac{{\sigma }_2}{{\sigma }_1}.$$

3. Dynamics Investigation on the Solution

In this segment, we are to explore the existence and uniqueness, non-negativeness, and boundedness of the solution for model (2) by virtue of Lemma 2 and Banach fixed point theorem.

Theorem 1.

Set $\mathrm{\Psi }=\{y_1,y_2)\in R^2:max\left\{\right|y_1|,|y_2\left|\right\}\le \mathcal{Y}\}$, where $\mathcal{Y}$ stands for a positive constant. For each $(y_{10},y_{20})\in \mathrm{\Psi }$, system (2) concerning the initial value $(y_{10},y_{20})$ possesses a unique solution $Y=(y_1,y_2)\in \mathrm{\Psi }.$

Proof. 

Set up the following mapping:

$$h\left(Y\right)=(h_1\left(Y\right),h_2\left(Y\right)),$$

(3)

where

$$\left\{\begin{array}{c}{\displaystyle h_1\left(Y\right)={\rho }_1{y}_1\left(t\right)[\alpha -2{\beta }_1\mu y_1\left(t\right)-2{\beta }_1(1-\mu )y_1(t-\theta )+{\delta }_1{y}_2\left(t\right)+{\beta }_1{\gamma }_1],}\hfill \\ {\displaystyle h_2\left(Y\right)={\rho }_2{y}_2\left(t\right)[\alpha -2{\beta }_2{y}_2\left(t\right)+{\delta }_2\mu y_1\left(t\right)+{\delta }_2(1-\mu )y_1(t-\theta )+{\beta }_2{\gamma }_2].}\hfill \end{array}\right.$$

(4)

For each $Y,\tilde{Y}\in \mathrm{\Psi }$, we get

$$\begin{array}{ccc}\hfill & \hfill & \left|\right|h\left(Y\right)-h\left(\tilde{Y}\right)\left|\right|\hfill \\ \hfill & \hfill & =\left|{\rho }_1{y}_1\left(t\right)[\alpha -2{\beta }_1\mu y_1\left(t\right)-2{\beta }_1(1-\mu )y_1(t-\theta )+{\delta }_1{y}_2\left(t\right)+{\beta }_1{\gamma }_1]\right.\hfill \\ \hfill & \hfill & -\left.\left\{{\rho }_1{\tilde{y}}_1\left(t\right)[\alpha -2{\beta }_1\mu {\tilde{y}}_1\left(t\right)-2{\beta }_1(1-\mu ){\tilde{y}}_1(t-\theta )+{\delta }_1{\tilde{y}}_2\left(t\right)+{\beta }_1{\gamma }_1]\right\}\right|\hfill \\ \hfill & \hfill & +\left|{\rho }_2{y}_2\left(t\right)[\alpha -2{\beta }_2{y}_2\left(t\right)+{\delta }_2\mu y_1\left(t\right)+{\delta }_2(1-\mu )y_1(t-\theta )+{\beta }_2{\gamma }_2]\right.\hfill \\ \hfill & \hfill & -\left.\left\{{\rho }_2{\tilde{y}}_2\left(t\right)[\alpha -2{\beta }_2{\tilde{y}}_2\left(t\right)+{\delta }_2\mu {\tilde{y}}_1\left(t\right)+{\delta }_2(1-\mu ){\tilde{y}}_1(t-\theta )+{\beta }_2{\gamma }_2]\right\}\right|\hfill \\ \hfill & \hfill & \le {\rho }_1\alpha |y_1\left(t\right)-{\tilde{y}}_1\left(t\right)|+4{\rho }_1{\beta }_1\mu \mathcal{Y}|y_1\left(t\right)-{\tilde{y}}_1\left(t\right)|+4{\rho }_1{\beta }_1(1-\mu )|y_1\left(t\right)-{\tilde{y}}_1\left(t\right)|\hfill \\ \hfill & \hfill & +{\rho }_1{\delta }_1\mathcal{Y}|y_1\left(t\right)-{\tilde{y}}_1\left(t\right)|+{\rho }_1{\delta }_1\mathcal{Y}|y_2\left(t\right)-{\tilde{y}}_2\left(t\right)|+{\rho }_1{\beta }_1{\gamma }_1|y_1\left(t\right)-{\tilde{y}}_1\left(t\right)|\hfill \\ \hfill & \hfill & +{\rho }_2\alpha |y_2\left(t\right)-{\tilde{y}}_2\left(t\right)|+4{\rho }_2{\beta }_2\mu \mathcal{Y}|y_2\left(t\right)-{\tilde{y}}_2\left(t\right)|+{\rho }_2{\delta }_2\mu \mathcal{Y}|y_1\left(t\right)-{\tilde{y}}_1\left(t\right)|\hfill \\ \hfill & \hfill & +{\rho }_2{\delta }_2\mu \mathcal{Y}|y_2\left(t\right)-{\tilde{y}}_2\left(t\right)|+{\rho }_2{\delta }_2(1-\mu )\mathcal{Y}|y_1\left(t\right)-{\tilde{y}}_1\left(t\right)|\hfill \\ \hfill & \hfill & +{\rho }_2{\delta }_2(1-\mu )\mathcal{Y}|y_2\left(t\right)-{\tilde{y}}_2\left(t\right)|+{\rho }_2{\beta }_2{\gamma }_2|y_2\left(t\right)-{\tilde{y}}_2\left(t\right)|\hfill \\ \hfill & \hfill & =\left[{\rho }_1\alpha +4{\rho }_1{\beta }_1\mu \mathcal{Y}+4{\rho }_1{\beta }_1(1-\mu )+{\rho }_1{\delta }_1\mathcal{Y}+{\rho }_1{\beta }_1{\gamma }_1+{\rho }_2{\delta }_2\mathcal{Y}\right]|y_1\left(t\right)-{\tilde{y}}_1\left(t\right)|\hfill \\ \hfill & \hfill & +[{\rho }_1{\delta }_1\mathcal{Y}+{\rho }_2\alpha +4{\rho }_2{\beta }_2\mu \mathcal{Y}+{\rho }_2{\delta }_2\mathcal{Y}+{\rho }_2{\beta }_2{\gamma }_2]|y_2\left(t\right)-{\tilde{y}}_2\left(t\right)|\hfill \\ \hfill & \hfill & ={\mathcal{Y}}_1|y_1\left(t\right)-{\tilde{y}}_1\left(t\right)|+{\mathcal{Y}}_2|y_2\left(t\right)-{\tilde{y}}_2\left(t\right)|\hfill \\ \hfill & \hfill & \le \mathcal{Y}\left|\right|Y-\tilde{Y}\left|\right|,\hfill \end{array}$$

(5)

where

$$\left\{\begin{array}{c}{\displaystyle {\mathcal{Y}}_1={\rho }_1\alpha +4{\rho }_1{\beta }_1\mu \mathcal{Y}+4{\rho }_1{\beta }_1(1-\mu )+{\rho }_1{\delta }_1\mathcal{Y}+{\rho }_1{\beta }_1{\gamma }_1+{\rho }_2{\delta }_2\mathcal{Y},}\hfill \\ {\displaystyle {\mathcal{Y}}_2={\rho }_1{\delta }_1\mathcal{Y}+{\rho }_2\alpha +4{\rho }_2{\beta }_2\mu \mathcal{Y}+{\rho }_2{\delta }_2\mathcal{Y}+{\rho }_2{\beta }_2{\gamma }_2}\hfill \end{array}\right.$$

(6)

and

$$\mathcal{Y}=max\{{\mathcal{Y}}_1,{\mathcal{Y}}_2\}.$$

(7)

Then $h\left(Y\right)$ obeys Lipschitz condition with respect to Y (one can see [39]). Taking advantage of Banach fixed point theorem, one concludes that Theorem 1 is true. □

Theorem 2.

(a) Every solution to system (2) starting with $R_{+}^2$ is non-negative; (b) If the following inequality

$$min\{2{\rho }_1{\beta }_1\mu ,2{\rho }_2{\beta }_2\}>\frac{1}{2}{\rho }_1{\delta }_1+\frac{1}{2}{\rho }_2{\delta }_2(1-\mu )$$

holds, then every solution to system (2) starting with $R_{+}^2$ is uniformly bounded.

Proof. 

Let the initial value of system (2) be $Y\left(t_0\right)=(y_1\left(t_0\right),y_2\left(t_0\right))$. Assume that ∃ a constant $t_{*}$ satisfying $t\in (t_0,t_{*})$ such that

$$\left\{\begin{array}{c}{\displaystyle y_1\left(t\right)=0,t\in (t_0,t_{*}),}\hfill \\ {\displaystyle y_1\left(t_{*}\right)=0,}\hfill \\ {\displaystyle y_1\left(t_{*}^{+}\right)<0.}\hfill \end{array}\right.$$

(8)

According to system (2), we have

$${\mathcal{D}}^{\eta }{y}_1{\left(t\right)|}_{y_1\left(t_{*}\right)=0}=0.$$

(9)

By Lemma of [40], one gets $y_1\left(t_{*}^{+}\right)=0$. By (8), we find that it is contradiction. Therefore, $y_1\left(t\right)\ge 0$, ∀ $t\ge t_0.$ Similarly, one can easily check that $y_2\left(t\right)\ge 0$, ∀ $t\ge t_0.$ The proof of (a) ends. In the sequel, we shall prove uniformly boundedness of system (2). Set

$$\mathrm{\Phi }\left(t\right)=y_1\left(t\right)+y_2\left(t\right).$$

(10)

Then

$$\begin{array}{ccc}{\mathcal{D}}^{\eta }\mathrm{\Phi }\left(t\right)+\kappa \mathrm{\Pi }\left(t\right)\hfill & =\hfill & {\mathcal{D}}^{\eta }{y}_1\left(t\right)+{\mathcal{D}}^{\eta }{y}_2\left(t\right)+\kappa y_1\left(t\right)+\kappa y_2\left(t\right)\hfill \\ \hfill & =\hfill & {\rho }_1{y}_1\left(t\right)[\alpha -2{\beta }_1\mu y_1\left(t\right)-2{\beta }_1(1-\mu )y_1(t-\theta )+{\delta }_1{y}_2\left(t\right)+{\beta }_1{\gamma }_1]\hfill \\ \hfill & \hfill & +{\rho }_2{y}_2\left(t\right)[\alpha -2{\beta }_2{y}_2\left(t\right)+{\delta }_2\mu y_1\left(t\right)+{\delta }_2(1-\mu )y_1(t-\theta )+{\beta }_2{\gamma }_2]\hfill \\ \hfill & \hfill & +\kappa y_1\left(t\right)+\kappa y_2\left(t\right)\hfill \\ \hfill & \le \hfill & ({\rho }_1\alpha +{\rho }_1{\beta }_1{\gamma }_1+\kappa )y_1\left(t\right)-2{\rho }_1{\beta }_1\mu y_1^2\left(t\right)+{\rho }_1{\delta }_1{y}_1\left(t\right)y_2\left(t\right)\hfill \\ \hfill & \hfill & +({\rho }_2\alpha +{\rho }_2{\beta }_2{\gamma }_1+\kappa )y_2\left(t\right)-2{\rho }_2{\beta }_2{y}_2^2\left(t\right)+{\rho }_2{\delta }_2(1-\mu )y_1(t-\theta )y_2\left(t\right)\hfill \\ \hfill & \le \hfill & ({\rho }_1\alpha +{\rho }_1{\beta }_1{\gamma }_1+\kappa )y_1\left(t\right)-2{\rho }_1{\beta }_1\mu y_1^2\left(t\right)+\frac{1}{2}{\rho }_1{\delta }_1{y}_1^2\left(t\right)+\frac{1}{2}{\rho }_1{\delta }_1{y}_2^2\left(t\right)\hfill \\ \hfill & \hfill & +({\rho }_2\alpha +{\rho }_2{\beta }_2{\gamma }_1+\kappa )y_2\left(t\right)-2{\rho }_2{\beta }_2{y}_2^2\left(t\right)+\frac{1}{2}{\rho }_2{\delta }_2(1-\mu )y_1^2(t-\theta )\hfill \\ \hfill & \hfill & +\frac{1}{2}{\rho }_2{\delta }_2(1-\mu )y_2^2\left(t\right)\hfill \\ \hfill & =\hfill & ({\rho }_1\alpha +{\rho }_1{\beta }_1{\gamma }_1+\kappa )y_1\left(t\right)-\left[2{\rho }_1{\beta }_1\mu -\frac{1}{2}{\rho }_1{\delta }_1-\frac{1}{2}{\rho }_2{\delta }_2(1-\mu )\right]y_1^2\left(t\right)\hfill \\ \hfill & \hfill & +({\rho }_2\alpha +{\rho }_2{\beta }_2{\gamma }_1+\kappa )y_2\left(t\right)-\left[2{\rho }_2{\beta }_2-\frac{1}{2}{\rho }_1{\delta }_1-\frac{1}{2}{\rho }_2{\delta }_2(1-\mu )\right]y_2^2\left(t\right)\hfill \\ \hfill & \le \hfill & \mathcal{Q},\hfill \end{array}$$

(11)

where $\kappa >0$ is a constant and

$$\mathcal{Q}=\frac{{\left[2{\rho }_1{\beta }_1\mu -\frac{1}{2}{\rho }_1{\delta }_1-\frac{1}{2}{\rho }_2{\delta }_2(1-\mu )\right]}^2}{4({\rho }_1\alpha +{\rho }_1{\beta }_1{\gamma }_1+\kappa )}+\frac{{\left[2{\rho }_2{\beta }_2-\frac{1}{2}{\rho }_1{\delta }_1-\frac{1}{2}{\rho }_2{\delta }_2(1-\mu )\right]}^2}{4({\rho }_2\alpha +{\rho }_2{\beta }_2{\gamma }_1+\kappa )}.$$

(12)

Then

$${\mathcal{D}}^{\eta }\mathrm{\Phi }\left(t\right)\le -\kappa \mathrm{\Phi }\left(t\right)+\mathcal{Q}.$$

(13)

According to Lemma 2, we get

$$\mathrm{\Phi }\left(t\right)\le \left(\mathrm{\Phi }\left(t_0\right)-\frac{\mathcal{Q}}{\kappa }\right)E_{\eta }[-\kappa {(t-t_0)}^{\eta }]+\frac{\mathcal{Q}}{\kappa },$$

(14)

then

$$\mathrm{\Phi }\left(t\right)\to \frac{\mathcal{Q}}{\kappa },ast\to \infty .$$

(15)

The proof of Theorem 2 finishes. □

4. Bifurcation Study

It is not difficult to know that model (2) owns the following equilibrium points:

$$Y_1(0,0),Y_2\left(0,\frac{\alpha +{\beta }_2{\gamma }_2}{2{\beta }_2}\right),Y_3\left(\frac{\alpha +{\beta }_1{\gamma }_1}{2{\beta }_1},0\right),Y_4(y_{1*},y_{2*}),$$

where

$$\left\{\begin{array}{c}{\displaystyle y_{1*}=\frac{2\alpha {\beta }_2+2{\beta }_1{\beta }_2{\gamma }_1+\alpha {\delta }_1+{\beta }_2{\gamma }_2{\delta }_1}{4{\beta }_1{\beta }_2-{\delta }_1{\delta }_2},}\hfill \\ {\displaystyle y_{2*}=\frac{2\alpha {\beta }_1+2{\beta }_1{\beta }_2{\gamma }_2+\alpha {\delta }_2+{\beta }_1{\gamma }_1{\delta }_2}{4{\beta }_1{\beta }_2-{\delta }_1{\delta }_2}.}\hfill \end{array}\right.$$

(16)

Suppose that

$$\left({\mathcal{S}}_1\right)4{\beta }_1{\beta }_2>{\delta }_1{\delta }_2$$

holds, then the equilibrium point $Y_4(y_{1*},y_{2*})$ is a positive equilibrium point. Considering the actual meaning in insurance market, we are only concerned with the positive equilibrium point $Y_4(y_{1*},y_{2*})$ of model (2). The linear system of model (2) near $Y_4(y_{1*},y_{2*})$ owns the expression:

$${\mathcal{D}}^{\eta }y\left(t\right)={\mathcal{A}}_{1}y\left(t\right)+{\mathcal{A}}_{2}y(t-\theta ),$$

(17)

where

$$y\left(t\right)=\left[\begin{array}{c}{y}_1\left(t\right)\\ y_2\left(t\right)\end{array}\right],{\mathcal{A}}_1=\left[\begin{array}{cc}{a}_1& a_2\\ a_3& a_4\end{array}\right],{\mathcal{A}}_2=\left[\begin{array}{cc}{a}_5& 0\\ a_6& 0\end{array}\right],$$

(18)

where

$$\left\{\begin{array}{c}{\displaystyle a_1={\rho }_1\alpha -2{\rho }_1{\beta }_1\mu y_{1*}+2{\rho }_1{\beta }_1{y}_{1*}+{\rho }_1{\delta }_1{y}_{2*}+{\rho }_1{\beta }_1{\gamma }_1,}\hfill \\ {\displaystyle a_2={\rho }_1{\beta }_1{y}_{1*},}\hfill \\ {\displaystyle a_3={\rho }_2{\beta }_2\mu y_{2*},}\hfill \\ {\displaystyle a_4={\rho }_2\alpha -4{\rho }_2{\beta }_2{y}_{2*}+{\rho }_2{\delta }_2{y}_{1*}+{\rho }_2{\beta }_2{\gamma }_2,}\hfill \\ {\displaystyle a_5=-2{\rho }_1{\beta }_1(1-\mu )y_{1*},}\hfill \\ {\displaystyle a_6={\rho }_2{\delta }_2(1-\mu )y_{2*}.}\hfill \end{array}\right.$$

(19)

The characteristic equation of system (17) owns the following expression:

$$det\left[\begin{array}{cc}{s}^{\eta }-a_1-a_5{e}^{-s\theta }& -a_2\\ -a_3-a_6{e}^{-s\theta }& s^{\eta }-a_4\end{array}\right]=0,$$

(20)

which generates

$$s^{2\eta }+b_1{s}^{\eta }+b_2+(b_3{s}^{\eta }+b_4)e^{-s\theta }=0,$$

(21)

where

$$\left\{\begin{array}{c}{\displaystyle b_1=-(a_1+a_4),}\hfill \\ {\displaystyle b_2=a_1{a}_4-a_2{a}_3,}\hfill \\ {\displaystyle b_3=-a_5,}\hfill \\ {\displaystyle b_4=a_4{a}_5-a_1{a}_6.}\hfill \end{array}\right.$$

(22)

When $\theta =0,$ then Equation (21) becomes:

$${\lambda }^2+(b_1+b_3)\lambda +b_2+b_4=0,$$

(23)

If

$$\left({\mathcal{S}}_2\right)b_1+b_3>0,b_2+b_4>0$$

is fulfilled, then the two roots ${\lambda }_1,{\lambda }_2$ of Equation (23) obey $|arg\left({\lambda }_1\right)|>\frac{\eta \pi }{2},\left|arg\left({\lambda }_2\right)\right|>\frac{\eta \pi }{2}$. It follows from Lemma 1 that the positive equilibrium point $Y_4(y_{1*},y_{2*})$ of model (2) involving $\theta =0$ is locally asymptotically stable.

Denote $s=i\varsigma =\varsigma \left(cos\frac{\pi }{2}+isin\frac{\pi }{2}\right)$ the root of Equation (21). Then, Equation (21) takes the following expression:

$$\begin{array}{ccc}\hfill & \hfill & {\varsigma }^{2\eta }(cos\eta \pi +isin\eta \pi )+b_1{\varsigma }^{\eta }\left(cos\frac{\eta \pi }{2}+isin\frac{\eta \pi }{2}\right)+b_2\\ \hfill & \hfill & +\left[b_3{\varsigma }^{\eta }\left(cos\frac{\eta \pi }{2}+isin\frac{\eta \pi }{2}\right)+b_4\right](cos\varsigma \theta -isin\varsigma \theta )=0.\end{array}$$

(24)

Then

$$\left\{\begin{array}{c}{\displaystyle {\mathcal{C}}_{1}cos\varsigma \theta +{\mathcal{C}}_{2}sin\varsigma \theta ={\mathcal{C}}_3,}\hfill \\ {\displaystyle {\mathcal{C}}_{2}cos\varsigma \theta -{\mathcal{C}}_{1}sin\varsigma \theta ={\mathcal{C}}_4,}\hfill \end{array}\right.$$

(25)

where

$$\left\{\begin{array}{c}{\displaystyle {\mathcal{C}}_1=b_3{\varsigma }^{\eta }cos\frac{\eta \pi }{2}+b_4,}\hfill \\ {\displaystyle {\mathcal{C}}_2=b_3{\varsigma }^{\eta }sin\frac{\eta \pi }{2},}\hfill \\ {\displaystyle {\mathcal{C}}_3=-{\varsigma }^{2\eta }cos\eta \pi -b_1{\varsigma }^{\eta }cos\frac{\eta \pi }{2}-b_2,}\hfill \\ {\displaystyle {\mathcal{C}}_4=-{\varsigma }^{2\eta }sin\eta \pi -b_1{\varsigma }^{\eta }sin\frac{\eta \pi }{2}.}\hfill \end{array}\right.$$

(26)

It follows from (25) that

$${\mathcal{C}}_1^2+{\mathcal{C}}_2^2={\mathcal{C}}_3^2+{\mathcal{C}}_4^2,$$

(27)

which implies

$${\varsigma }^{4\eta }+{\tau }_1{\varsigma }^{3\eta }+{\tau }_2{\varsigma }^{2\eta }+{\tau }_3{\varsigma }^{\eta }+{\tau }_4=0,$$

(28)

where

$$\left\{\begin{array}{c}{\displaystyle {\tau }_1=2b_1\left(cos\eta \pi cos\frac{\eta \pi }{2}+sin\eta \pi sin\frac{\eta \pi }{2}\right),}\hfill \\ {\displaystyle {\tau }_2=2b_{2}cos\eta \pi +b_1^2-b_3^2,}\hfill \\ {\displaystyle {\tau }_3=2(b_1{b}_2-b_3{b}_4)cos\frac{\eta \pi }{2},}\hfill \\ {\displaystyle {\tau }_4=b_2^2-b_4^2.}\hfill \end{array}\right.$$

(29)

Let

$$\mathsf{{\rm Y}}\left(\varsigma \right)={\varsigma }^{4\eta }+{\tau }_1{\varsigma }^{3\eta }+{\tau }_2{\varsigma }^{2\eta }+{\tau }_3{\varsigma }^{\eta }+{\tau }_4.$$

(30)

Suppose that

$$\left({\mathcal{S}}_3\right)|b_2|<|b_4|$$

holds true, considering $\frac{d\mathsf{{\rm Y}}\left(\varsigma \right)}{d\varsigma }>0$,∀$\varsigma >0$, then one can easily know that Equation (28) has at least one positive real root. Therefore, Equation (21) owns at least one pair of purely roots. Making use of Sun et al. [41], the following assertion can be easily available.

Lemma 3.

(i) If ${\tau }_l>0(l=1,2,3,4)$ holds, then Equation (21) owns no root with zero real parts provided that $\theta \ge 0$. (ii) If $\left({\mathcal{S}}_3\right)$ holds and ${\tau }_l>0(l=1,2,3)$, then Equation (21) has a pair of purely imaginary roots $\pm i{\varsigma }_0$ for $\theta ={\theta }_h(h=1,2,\cdots ,)$ where

$${\theta }_0^{\left(h\right)}=\frac{1}{{\varsigma }_0}\left[arccos\left(\frac{{\mathcal{C}}_1{\mathcal{C}}_3+{\mathcal{C}}_2{\mathcal{C}}_4}{{\mathcal{C}}_1^2+{\mathcal{C}}_2^2}\right)+2h\pi \right],$$

(31)

where $h=0,1,2,\cdots ,$ and ${\varsigma }_0>0$ represents the unique zero of $\mathsf{{\rm Y}}\left(\varsigma \right).$

Here we omit the concrete proof of Lemma 1, one can consult [41] . Denote ${\theta }_0={\theta }_0^{\left(0\right)}$.

In the sequel, we make the following necessary assumption:

$$\left({\mathcal{S}}_4\right){\mathcal{I}}_{1R}{\mathcal{I}}_{2R}+{\mathcal{I}}_{1I}{\mathcal{I}}_{2I}>0,$$

where

$$\left\{\begin{array}{c}{\displaystyle {\mathcal{I}}_{1R}=2\eta {\varsigma }_0^{2\eta -1}cos\frac{(2\eta -1)\pi }{2}+\eta b_1{\varsigma }_0^{\eta -1}cos\frac{(\eta -1)\pi }{2}}\hfill \\ {\displaystyle +\eta b_3{\varsigma }_0^{\eta -1}\left[cos\frac{(\eta -1)\pi }{2}cos{\varsigma }_0{\theta }_0+sin\frac{(\eta -1)\pi }{2}sin{\varsigma }_0{\theta }_0\right],}\hfill \\ {\displaystyle {\mathcal{I}}_{1I}=2\eta {\varsigma }_0^{2\eta -1}sin\frac{(2\eta -1)\pi }{2}+\eta b_1{\varsigma }_0^{\eta -1}sin\frac{(\eta -1)\pi }{2}}\hfill \\ {\displaystyle -\eta b_3{\varsigma }_0^{\eta -1}\left[cos\frac{(\eta -1)\pi }{2}sin{\varsigma }_0{\theta }_0-sin\frac{(\eta -1)\pi }{2}cos{\varsigma }_0{\theta }_0\right],}\hfill \\ {\displaystyle {\mathcal{I}}_{2R}=\left(b_3{\varsigma }_0^{\eta }cos\frac{\eta \pi }{2}+b_4\right){\varsigma }_{0}sin{\varsigma }_0{\theta }_0+\left(b_3{\varsigma }_0^{\eta }cos\frac{\eta \pi }{2}\right){\varsigma }_{0}cos{\varsigma }_0{\theta }_0,}\hfill \\ {\displaystyle {\mathcal{I}}_{2I}=\left(b_3{\varsigma }_0^{\eta }cos\frac{\eta \pi }{2}+b_4\right){\varsigma }_{0}sin{\varsigma }_0{\theta }_0-\left(b_3{\varsigma }_0^{\eta }cos\frac{\eta \pi }{2}\right){\varsigma }_{0}cos{\varsigma }_0{\theta }_0.}\hfill \end{array}\right.$$

(32)

Lemma 4.

Let $s\left(\theta \right)={\omega }_1\left(\theta \right)+i{\omega }_2\left(\theta \right)$ be the root of Equation (21) near $\theta ={\theta }_0$ such that ${\omega }_1\left({\theta }_0\right)=0,{\omega }_2\left({\theta }_0\right)={\varsigma }_0$, then $Re\left(\frac{ds}{d\theta }\right){|}_{\theta ={\theta }_0,\varsigma ={\varsigma }_0}>0.$

Proof. 

Using Equation (21), one derives

$$\begin{array}{ccc}\hfill & \hfill & \left[2\eta s^{2\eta -1}+\eta b_1{s}^{\eta -1}\right]\frac{ds}{d\theta }+\eta b_3{s}^{\eta -1}{e}^{-s\theta }\frac{ds}{d\theta }\hfill \\ \hfill & \hfill & -e^{-s\theta }\left(\frac{ds}{d\theta }\theta +s\right)\left(b_3{s}^{\eta }+b_4\right)=0,\hfill \end{array}$$

(33)

which implies

$${\left(\frac{ds}{d\theta }\right)}^{-1}=\frac{{\mathcal{I}}_1\left(s\right)}{{\mathcal{I}}_2\left(s\right)}-\frac{\theta }{s},$$

(34)

where

$$\left\{\begin{array}{c}{\displaystyle {\mathcal{I}}_1\left(s\right)=2\eta s^{2\eta -1}+\eta b_1{s}^{\eta -1}+\eta b_3{s}^{\eta -1}{e}^{-s\theta },}\hfill \\ {\displaystyle {\mathcal{I}}_2\left(s\right)=s{e}^{-s\theta }[b_3{s}^{\eta }+b_4].}\hfill \end{array}\right.$$

(35)

Hence

$$Re{\left[{\left(\frac{ds}{d\theta }\right)}^{-1}\right]}_{\theta ={\theta }_0,\varsigma ={\varsigma }_0}=Re{\left[\frac{{\mathcal{I}}_1\left(s\right)}{{\mathcal{I}}_2\left(s\right)}\right]}_{\theta ={\theta }_0,\varsigma ={\varsigma }_0}=\frac{{\mathcal{I}}_{1R}{\mathcal{I}}_{2R}+{\mathcal{I}}_{1I}{\mathcal{I}}_{2I}}{{\mathcal{I}}_{2R}^2+{\mathcal{I}}_{2I}^2}.$$

(36)

Taking advantage of $\left({\mathcal{S}}_4\right)$, one gets

$$Re{\left[{\left(\frac{ds}{d\theta }\right)}^{-1}\right]}_{\theta ={\theta }_0,\varsigma ={\varsigma }_0}>0,$$

(37)

which completes the proof. □

Making use of Lemma 1, one gets the following result.

Theorem 3.

If $\left({\mathcal{S}}_1\right)$–$\left({\mathcal{S}}_4\right)$ holds, then $Y_4(y_{1*},y_{2*})$ of model (2) is locally asymptotically stable provided that $\theta \in [0,{\theta }_0)$ and a Hopf bifurcation of model (2) arises near $Y_4(y_{1*},y_{2*})$ for $\theta ={\theta }_0.$

5. Global Asymptotic Stability Exploration

In this part, we will explore the global stability issue of the positive equilibrium point $Y(y_{1*},y_{2*})$ of model (2). Firstly, we give the following assumption:

$$\left({\mathcal{S}}_5\right){\rho }_1{\rho }_2{\beta }_1{\beta }_2<\frac{{({\rho }_1{\delta }_1+{\rho }_2{\delta }_2)}^2}{4}.$$

Theorem 4.

If $\left({\mathcal{S}}_5\right)$ is fulfilled, then the positive equilibrium point $Y(y_{1*},y_{2*})$ of model (2) is globally asymptotically stable.

Proof. 

Setting up the following positive definite function:

$$\mathcal{V}\left(t\right)=\sum _{l=1}^2\left(y_l\left(t\right)-y_{l*}-y_{l*}ln\frac{y_l\left(t\right)}{y_{l*}}\right).$$

(38)

Then

$$\begin{array}{ccc}{\mathcal{D}}^{\eta }\mathcal{V}\left(t\right)\hfill & =\hfill & \frac{y_1\left(t\right)-y_{1*}}{y_1\left(t\right)}{\mathcal{D}}^{\eta }{w}_1\left(t\right)+\frac{y_2\left(t\right)-y_{2*}}{y_2\left(t\right)}{\mathcal{D}}^{\eta }{y}_2\left(t\right)\hfill \\ \hfill & \le \hfill & (y_1\left(t\right)-y_{1*})\left[{\rho }_1(\alpha -2{\beta }_1\mu y_1\left(t\right)-2{\beta }_1(1-\mu )y_1(t-\theta )+{\delta }_1{y}_2\left(t\right)+{\beta }_1{\gamma }_1)\right]\hfill \\ \hfill & \hfill & +(y_2\left(t\right)-y_{2*})\left[{\rho }_2(\alpha -2{\beta }_2{y}_2\left(t\right)+{\delta }_2\mu y_1\left(t\right)+{\delta }_2(1-\mu )y_1(t-\theta )+{\beta }_2{\gamma }_2)\right]\hfill \\ \hfill & =\hfill & ((y_1\left(t\right)-y_{1*})[-2{\rho }_1{\beta }_1\mu y_1\left(t\right)+2{\rho }_1{\beta }_1\mu y_{1*}-2{\rho }_1{\beta }_1(1-\mu )y_1(t-\theta )\hfill \\ \hfill & \hfill & +2{\rho }_1{\beta }_1(1-\mu )y_{1*}+{\rho }_1{\delta }_1{y}_2\left(t\right)-{\rho }_1{\delta }_1{y}_{2*}]\hfill \\ \hfill & \hfill & +(y_2\left(t\right)-y_{2*})[-2{\rho }_2{\beta }_2{y}_2\left(t\right)+2{\rho }_2{\beta }_2{y}_{2*}+{\rho }_2{\delta }_2\mu y_1\left(t\right)\hfill \\ \hfill & \hfill & -{\rho }_2{\delta }_2\mu y_{1*}+{\rho }_2{\delta }_2(1-\mu )y_1(t-\theta )-{\rho }_2{\delta }_2(1-\mu )y_{1*}]\hfill \\ \hfill & \le \hfill & -{\rho }_1{\beta }_1{(y_1\left(t\right)-y_{1*})}^2+{\rho }_1{\delta }_1(y_1\left(t\right)-y_{1*})(y_2\left(t\right)-y_{2*})\hfill \\ \hfill & \hfill & -{\rho }_2{\beta }_2{(y_2\left(t\right)-y_{2*})}^2+{\rho }_2{\delta }_2(y_1\left(t\right)-y_{1*})(y_2\left(t\right)-y_{2*})\hfill \\ \hfill & =\hfill & -{\rho }_1{\beta }_1{(y_1\left(t\right)-y_{1*})}^2-{\rho }_2{\beta }_2{(y_2\left(t\right)-y_{2*})}^2\hfill \\ \hfill & \hfill & +({\rho }_1{\delta }_1+{\rho }_2{\delta }_2)(y_1\left(t\right)-y_{1*})(y_2\left(t\right)-y_{2*}).\hfill \end{array}$$

(39)

By $\left({\mathcal{S}}_5\right)$, we can know that ${\mathcal{D}}^{\eta }\mathcal{V}\left(t\right)\le 0$, which completes the proof. □

6. Hybrid Control Technique for Bifurcation Control

In this part, we will make use of an appropriate hybrid controller which consists of state feedback and parameter perturbation to control the stability and Hopf bifurcation for model (2). By virtue of the research ideas in [42,43,44], we get the fractional-order controlled duopoly game model:

$$\left\{\begin{array}{c}{\displaystyle {\mathcal{D}}^{\eta }{y}_1\left(t\right)={\sigma }_1\left\{{\rho }_1{y}_1\left(t\right)[\alpha -2{\beta }_1\mu y_1\left(t\right)-2{\beta }_1(1-\mu )y_1(t-\theta )+{\delta }_1{y}_2\left(t\right)+{\beta }_1{\gamma }_1]\right\}}\hfill \\ {\displaystyle +{\sigma }_2(y_1\left(t\right)-y_{1*}),}\hfill \\ {\displaystyle {\mathcal{D}}^{\eta }{y}_1\left(t\right)={\sigma }_1\left\{{\rho }_2{y}_2\left(t\right)[\alpha -2{\beta }_2{y}_2\left(t\right)+{\delta }_2\mu y_1\left(t\right)+{\delta }_2(1-\mu )y_1(t-\theta )+{\beta }_2{\gamma }_2]\right\}}\hfill \\ {\displaystyle +{\sigma }_2(y_2\left(t\right)-y_{2*}),}\hfill \end{array}\right.$$

(40)

where ${\sigma }_1,{\sigma }_2$ represent feedback gain parameters. Models (40) and (2) own the same equilibrium points

$$Y_1(0,0),Y_2\left(0,\frac{\alpha +{\beta }_2{\gamma }_2}{2{\beta }_2}\right),Y_3\left(\frac{\alpha +{\beta }_1{\gamma }_1}{2{\beta }_1},0\right),Y_4(y_{1*},y_{2*}).$$

If $\left({\mathcal{S}}_1\right)$ holds, then $Y_4(y_{1*},y_{2*})$ is positive equilibrium point. The linear system of model (41) near $Y_4(y_{1*},y_{2*})$ owns the expression:

$${\mathcal{D}}^{\eta }y\left(t\right)={\mathcal{B}}_{1}y\left(t\right)+{\mathcal{B}}_{2}y(t-\theta ),$$

(41)

where

$$y\left(t\right)=\left[\begin{array}{c}{y}_1\left(t\right)\\ y_2\left(t\right)\end{array}\right],{\mathcal{B}}_1=\left[\begin{array}{cc}{g}_1& g_2\\ g_3& g_4\end{array}\right],{\mathcal{B}}_2=\left[\begin{array}{cc}{g}_5& 0\\ g_6& 0\end{array}\right],$$

(42)

where

$$\left\{\begin{array}{c}{\displaystyle g_1={\sigma }_1[{\rho }_1\alpha -2{\rho }_1{\beta }_1\mu y_{1*}+2{\rho }_1{\beta }_1{y}_{1*}+{\rho }_1{\delta }_1{y}_{2*}+{\rho }_1{\beta }_1{\gamma }_1]+{\sigma }_2,}\hfill \\ {\displaystyle g_2={\sigma }_1{\rho }_1{\beta }_1{y}_{1*},}\hfill \\ {\displaystyle g_3={\sigma }_1{\rho }_2{\beta }_2\mu y_{2*},}\hfill \\ {\displaystyle g_4={\sigma }_1[{\rho }_2\alpha -4{\rho }_2{\beta }_2{y}_{2*}+{\rho }_2{\delta }_2{y}_{1*}+{\rho }_2{\beta }_2{\gamma }_2]+{\sigma }_2,}\hfill \\ {\displaystyle g_5=-2{\sigma }_1{\rho }_1{\beta }_1(1-\mu )y_{1*},}\hfill \\ {\displaystyle g_6={\sigma }_1{\rho }_2{\delta }_2(1-\mu )y_{2*}.}\hfill \end{array}\right.$$

(43)

The characteristic equation of system (41) owns the following expression:

$$det\left[\begin{array}{cc}{s}^{\eta }-g_1-a_5{e}^{-s\theta }& -g_2\\ -g_3-g_6{e}^{-s\theta }& s^{\eta }-g_4\end{array}\right]=0,$$

(44)

which generates

$$s^{2\eta }+h_1{s}^{\eta }+h_2+(h_3{s}^{\eta }+h_4)e^{-s\theta }=0,$$

(45)

where

$$\left\{\begin{array}{c}{\displaystyle h_1=-(g_1+g_4),}\hfill \\ {\displaystyle h_2=g_1{g}_4-g_2{g}_3,}\hfill \\ {\displaystyle h_3=-g_5,}\hfill \\ {\displaystyle h_4=g_4{g}_5-g_1{g}_6.}\hfill \end{array}\right.$$

(46)

When $\theta =0,$ then Equation (45) becomes:

$${\lambda }^2+(h_1+h_3)\lambda +h_2+h_4=0,$$

(47)

If

$$\left({\mathcal{S}}_6\right)h_1+h_3>0,h_2+h_4>0$$

is fulfilled, then the two roots ${\lambda }_1,{\lambda }_2$ of Equation (47) obey $|arg\left({\lambda }_1\right)|>\frac{\eta \pi }{2},\left|arg\left({\lambda }_2\right)\right|>\frac{\eta \pi }{2}$. It follows from Lemma 1 that the positive equilibrium point $Y_4(y_{1*},y_{2*})$ of model (40) involving $\theta =0$ is locally asymptotically stable.

Denote $s=i\xi =\xi \left(cos\frac{\pi }{2}+isin\frac{\pi }{2}\right)$ the root of Equation (45). Then, Equation (45) takes the following expression:

$$\begin{array}{ccc}\hfill & \hfill & {\xi }^{2\eta }(cos\eta \pi +isin\eta \pi )+h_1{\xi }^{\eta }\left(cos\frac{\eta \pi }{2}+isin\frac{\eta \pi }{2}\right)+h_2\hfill \\ \hfill & \hfill & +\left[h_3{\xi }^{\eta }\left(cos\frac{\eta \pi }{2}+isin\frac{\eta \pi }{2}\right)+h_4\right](cos\xi \theta -isin\xi \theta )=0.\hfill \end{array}$$

(48)

Then

$$\left\{\begin{array}{c}{\displaystyle {\mathcal{G}}_{1}cos\xi \theta +{\mathcal{G}}_{2}sin\xi \theta ={\mathcal{G}}_3,}\hfill \\ {\displaystyle {\mathcal{G}}_{2}cos\xi \theta -{\mathcal{G}}_{1}sin\xi \theta ={\mathcal{G}}_4,}\hfill \end{array}\right.$$

(49)

where

$$\left\{\begin{array}{c}{\displaystyle {\mathcal{G}}_1=h_3{\xi }^{\eta }cos\frac{\eta \pi }{2}+h_4,}\hfill \\ {\displaystyle {\mathcal{G}}_2=h_3{\xi }^{\eta }sin\frac{\eta \pi }{2},}\hfill \\ {\displaystyle {\mathcal{G}}_3=-{\xi }^{2\eta }cos\eta \pi -h_1{\xi }^{\eta }cos\frac{\eta \pi }{2}-h_2,}\hfill \\ {\displaystyle {\mathcal{G}}_4=-{\xi }^{2\eta }sin\eta \pi -h_1{\xi }^{\eta }sin\frac{\eta \pi }{2}.}\hfill \end{array}\right.$$

(50)

It follows from (49) that

$${\mathcal{G}}_1^2+{\mathcal{G}}_2^2={\mathcal{G}}_3^2+{\mathcal{G}}_4^2,$$

(51)

which implies

$${\xi }^{4\eta }+{\iota }_1{\xi }^{3\eta }+{\iota }_2{\xi }^{2\eta }+{\iota }_3{\xi }^{\eta }+{\iota }_4=0,$$

(52)

where

$$\left\{\begin{array}{c}{\displaystyle {\iota }_1=2h_1\left(cos\eta \pi cos\frac{\eta \pi }{2}+sin\eta \pi sin\frac{\eta \pi }{2}\right),}\hfill \\ {\displaystyle {\iota }_2=2h_{2}cos\eta \pi +h_1^2-h_3^2,}\hfill \\ {\displaystyle {\iota }_3=2(h_1{h}_2-h_3{h}_4)cos\frac{\eta \pi }{2},}\hfill \\ {\displaystyle {\iota }_4=h_2^2-h_4^2.}\hfill \end{array}\right.$$

(53)

Let

$$\mathrm{\Pi }\left(\xi \right)={\xi }^{4\eta }+{\iota }_1{\xi }^{3\eta }+{\iota }_2{\xi }^{2\eta }+{\iota }_3{\xi }^{\eta }+{\iota }_4.$$

(54)

Suppose that

$$\left({\mathcal{S}}_7\right)|h_2|<|h_4|$$

holds true, considering $\frac{d\mathrm{\Pi }\left(\xi \right)}{d\xi }>0$,∀$\xi >0$, then one can easily know that Equation (52) has at least one positive real root. Therefore, Equation (45) owns at least one pair of purely roots. Making use of Sun et al. [41], the following assertion can be easily available.

Lemma 5.

(i) If ${\iota }_l>0(l=1,2,3,4)$ holds, then Equation (45) owns no root with zero real parts provided that $\theta \ge 0$. (ii) If $\left({\mathcal{S}}_7\right)$ holds and ${\iota }_l>0(l=1,2,3)$, then Equation (45) has a pair of purely imaginary roots $\pm i{\xi }_0$ for $\theta ={\theta }_h(h=1,2,\cdots ,)$ where

$${\theta }_0^{\left(h\right)}=\frac{1}{{\xi }_0}\left[arccos\left(\frac{{\mathcal{G}}_1{\mathcal{G}}_3+{\mathcal{G}}_2{\mathcal{G}}_4}{{\mathcal{G}}_1^2+{\mathcal{G}}_2^2}\right)+2h\pi \right],$$

(55)

where $h=0,1,2,\cdots ,$ and ${\xi }_0>0$ represents the unique zero of $\mathrm{\Pi }\left(\xi \right).$

Here we omit the concrete proof of Lemma 5; one can consult [41]. Denote ${\theta }_{0*}={\theta }_0^{\left(0\right)}$. In the sequel, we make the following necessary assumption:

$$\left({\mathcal{S}}_8\right){\mathcal{I}}_{1R}{\mathcal{I}}_{2R}+{\mathcal{I}}_{1I}{\mathcal{I}}_{2I}>0,$$

where

$$\left\{\begin{array}{c}{\displaystyle {\mathcal{I}}_{1R}=2\eta {\varsigma }_0^{2\eta -1}cos\frac{(2\eta -1)\pi }{2}+\eta h_1{\xi }_0^{\eta -1}cos\frac{(\eta -1)\pi }{2}}\hfill \\ {\displaystyle +\eta h_3{\xi }_0^{\eta -1}\left[cos\frac{(\eta -1)\pi }{2}cos{\xi }_0{\theta }_{0*}+sin\frac{(\eta -1)\pi }{2}sin{\xi }_0{\theta }_{0*}\right],}\hfill \\ {\displaystyle {\mathcal{I}}_{1I}=2\eta {\xi }_0^{2\eta -1}sin\frac{(2\eta -1)\pi }{2}+\eta h_1{\xi }_0^{\eta -1}sin\frac{(\eta -1)\pi }{2}}\hfill \\ {\displaystyle -\eta h_3{\xi }_0^{\eta -1}\left[cos\frac{(\eta -1)\pi }{2}sin{\xi }_0{\theta }_{0*}-sin\frac{(\eta -1)\pi }{2}cos{\xi }_0{\theta }_{0*}\right],}\hfill \\ {\displaystyle {\mathcal{I}}_{2R}=\left(h_3{\xi }_0^{\eta }cos\frac{\eta \pi }{2}+h_4\right){\xi }_{0}sin{\xi }_0{\theta }_{0*}+\left(h_3{\xi }_0^{\eta }cos\frac{\eta \pi }{2}\right){\xi }_{0}cos{\xi }_0{\theta }_{0*},}\hfill \\ {\displaystyle {\mathcal{I}}_{2I}=\left(h_3{\xi }_0^{\eta }cos\frac{\eta \pi }{2}+h_4\right){\xi }_{0}sin{\xi }_0{\theta }_{0*}-\left(h_3{\xi }_0^{\eta }cos\frac{\eta \pi }{2}\right){\xi }_{0}cos{\xi }_0{\theta }_{0*}.}\hfill \end{array}\right.$$

(56)

Lemma 6.

Let $s\left(\theta \right)={\varpi }_1\left(\theta \right)+i{\varpi }_2\left(\theta \right)$ be the root of Equation (45) near $\theta ={\theta }_{0*}$ such that ${\varpi }_1\left({\theta }_{0*}\right)=0,{\varpi }_2\left({\theta }_{0*}\right)={\xi }_0$, then $Re\left(\frac{ds}{d\theta }\right){|}_{\theta ={\theta }_{0*},\xi ={\xi }_0}>0.$

Proof. 

Using Equation (45), one derives

$$\begin{array}{ccc}\hfill & \hfill & \left[2\eta s^{2\eta -1}+\eta h_1{s}^{\eta -1}\right]\frac{ds}{d\theta }+\eta h_3{s}^{\eta -1}{e}^{-s\theta }\frac{ds}{d\theta }\hfill \\ \hfill & \hfill & -e^{-s\theta }\left(\frac{ds}{d\theta }\theta +s\right)\left(h_3{s}^{\eta }+h_4\right)=0,\hfill \end{array}$$

(57)

which implies

$${\left(\frac{ds}{d\theta }\right)}^{-1}=\frac{{\mathcal{I}}_1\left(s\right)}{{\mathcal{I}}_2\left(s\right)}-\frac{\theta }{s},$$

(58)

where

$$\left\{\begin{array}{c}{\displaystyle {\mathcal{I}}_1\left(s\right)=2\eta s^{2\eta -1}+\eta h_1{s}^{\eta -1}+\eta h_3{s}^{\eta -1}{e}^{-s\theta },}\hfill \\ {\displaystyle {\mathcal{I}}_2\left(s\right)=s{e}^{-s\theta }[h_3{s}^{\eta }+h_4].}\hfill \end{array}\right.$$

(59)

Hence,

$$Re{\left[{\left(\frac{ds}{d\theta }\right)}^{-1}\right]}_{\theta ={\theta }_{0*},\xi ={\xi }_0}=Re{\left[\frac{{\mathcal{I}}_1\left(s\right)}{{\mathcal{I}}_2\left(s\right)}\right]}_{\theta ={\theta }_{0*},\xi ={\xi }_0}=\frac{{\mathcal{I}}_{1R}{\mathcal{I}}_{2R}+{\mathcal{I}}_{1I}{\mathcal{I}}_{2I}}{{\mathcal{I}}_{2R}^2+{\mathcal{I}}_{2I}^2}.$$

(60)

Taking advantage of $\left({\mathcal{S}}_8\right)$, one gets

$$Re{\left[{\left(\frac{ds}{d\theta }\right)}^{-1}\right]}_{\theta ={\theta }_{0*},\xi ={\xi }_0}>0,$$

(61)

which completes the proof. □

Making use of Lemma 1, one gets the following result.

Theorem 5.

If $\left({\mathcal{S}}_1\right),\left({\mathcal{S}}_6\right),\left({\mathcal{S}}_7\right)$ and $\left({\mathcal{S}}_8\right)$ hold, then $Y_4(y_{1*},y_{2*})$ of model (40) is locally asymptotically stable provided that $\theta \in [0,{\theta }_{0*})$ and a Hopf bifurcation of model (40) arises near $Y_4(y_{1*},y_{2*})$ for $\theta ={\theta }_{0*}.$

Theorem 6.

In 2012, Xu and Ma [18] explored the local stability and the creation of Hopf bifurcation of integer-order (1). In this current work, we mainly explore the various dynamics including the existence and uniqueness, non-negativeness, boundedness of the solution, local stability, onset of Hopf bifurcation, and Hopf bifurcation control problem for the established fractional-order delayed duopoly game model (2), which comes from the modified version of integer-order delayed duopoly game model (1). All investigation approaches and ideas practically differ from those in Xu and Ma [18] . The exploration idea of Xu and Ma [18] can not be applied to study the dynamical characteristics of model (2) in this work. From this viewpoint, we hold that our works replenish the research of [18] and expedite the development of bifurcation principle of fractional differential system.

7. Software Simulations

Example 1.

Consider the following fractional-order delayed duopoly game model:

$$\left\{\begin{array}{c}{\displaystyle {\mathcal{D}}^{0.92}{y}_1\left(t\right)={\rho }_1{y}_1\left(t\right)[\alpha -2{\beta }_1\mu y_1\left(t\right)-2{\beta }_1(1-\mu )y_1(t-\theta )+{\delta }_1{y}_2\left(t\right)+{\beta }_1{\gamma }_1],}\hfill \\ {\displaystyle {\mathcal{D}}^{0.92}{y}_1\left(t\right)={\rho }_2{y}_2\left(t\right)[\alpha -2{\beta }_2{y}_2\left(t\right)+{\delta }_2\mu y_1\left(t\right)+{\delta }_2(1-\mu )y_1(t-\theta )+{\beta }_2{\gamma }_2],}\hfill \end{array}\right.$$

(62)

where ${\rho }_1=0.2,{\rho }_2=0.2,\alpha =5,{\beta }_1=4.5,{\beta }_2=5,\mu =0.2,{\delta }_1=0.7,{\delta }_2=0.6,{\gamma }_1=\frac{1}{90},$${\gamma }_2=0.0001.$ By algebraic operation, one can derivethatthe unique positive equilibrium point of system (62) takes the value $Y(0.5978,0.5359)$. It is checked that the hypotheses $\left({\mathcal{S}}_1\right)$–$\left({\mathcal{S}}_4\right)$ in Theorem 3 are all met. Utilizing Matlab software, one can determine that ${\varsigma }_0=5.3122,{\theta }_0=1.7$. To test the correctness of the key conclusions of Theorem 3, in the sequel, we will fix two delay values. Firstly, set $\theta =1.52$ which is less than ${\theta }_0=1.7$, namely, θ falls into the range of value $[0,{\theta }_0)$. For this case, the Matlab simulation plots are provided in Figure 1. Apparently, Figure 1 demonstrates that the price of the first insurance company $y_1$ will approach to $0.5978$ and the price of the second insurance company $y_2$ will approach to $0.5359$ with the increase of time t. Secondly, set $\theta =1.94$ which is greater than ${\theta }_0=1.7$, namely, θ exceeds the key value ${\theta }_0)$. For this case, the Matlab simulation plots are provided in Figure 2. Apparently, Figure 2 demonstrates that the price of the first insurance company $y_1$ will oscillate around the value $0.5978$ and the price of the second insurance company $y_2$ will oscillate around the value $0.5359$ with the increase of time t. That is to say, a Hopf bifurcation (a limit cycle) will take place near the equilibrium point $Y(0.5978,0.5359)$. In addition,in order to intuitively display the bifurcation value of delay, we give the bifurcation diagrams thatshowthe bifurcation point ${\theta }_0\approx 1.7$ (see Figure 3 and Figure 4).

Example 2.

Consider the following fractional-order controlled delayed duopoly game model:

$$\left\{\begin{array}{c}{\displaystyle {\mathcal{D}}^{0.92}{y}_1\left(t\right)={\sigma }_1\left\{{\rho }_1{y}_1\left(t\right)[\alpha -2{\beta }_1\mu y_1\left(t\right)-2{\beta }_1(1-\mu )y_1(t-\theta )+{\delta }_1{y}_2\left(t\right)+{\beta }_1{\gamma }_1]\right\}}\hfill \\ {\displaystyle +{\sigma }_2(y_1\left(t\right)-y_{1*}),}\hfill \\ {\displaystyle {\mathcal{D}}^{0.92}{y}_1\left(t\right)={\sigma }_1\left\{{\rho }_2{y}_2\left(t\right)[\alpha -2{\beta }_2{y}_2\left(t\right)+{\delta }_2\mu y_1\left(t\right)+{\delta }_2(1-\mu )y_1(t-\theta )+{\beta }_2{\gamma }_2]\right\}}\hfill \\ {\displaystyle +{\sigma }_2(y_2\left(t\right)-y_{2*}).}\hfill \end{array}\right.$$

(63)

where ${\rho }_1=0.2,{\rho }_2=0.2,\alpha =5,{\beta }_1=4.5,{\beta }_2=5,\mu =0.2,{\delta }_1=0.7,{\delta }_2=0.6,{\gamma }_1=\frac{1}{90},$${\gamma }_2=0.0001.$ Let ${\sigma }_1=0.2,{\sigma }_2=0.4$. By algebraic operation, one can derive the unique positive equilibrium point of system (63) takes the value $Y(0.5978,0.5359)$. It is checked that the hypotheses $\left({\mathcal{S}}_1\right)$–$\left({\mathcal{S}}_8\right)$ in Theorem 5 are all met. Utilizing Matlab software, one can determine that ${\xi }_0=4.007,{\theta }_{0*}=1.33$. To test the correctness of the key conclusions of Theorem 5, in the sequel, we will fix two delay values. Firstly, set $\theta =1.1$ which is less than ${\theta }_{0*}=1.33$, namely, θ falls into the range of value $[0,{\theta }_{0*})$. For this case, the Matlab simulation plots are provided in Figure 5. Apparently, Figure 5 demonstrates that the price of the first insurance company $y_1$ will approach to $0.5978$ and the price of the second insurance company $y_2$ will approach to $0.5359$ with the increase of time t. Secondly, set $\theta =1.45$ which is greater than ${\theta }_{0*}=1.33$, namely, θ exceeds the key value ${\theta }_0)$. For this case, the Matlab simulation plots are provided in Figure 6. Apparently, Figure 6 demonstrates that the price of the first insurance company $y_1$ will oscillate around the value $0.5978$ and the price of the second insurance company $y_2$ will oscillate around the value $0.5359$ with the increase of time t. Thai is to say, a Hopf bifurcation (a limit cycle) will take place near the equilibrium point $Y(0.5978,0.5359)$. In addition,in order to intuitively display the bifurcation value of delay, we give the bifurcation diagrams thatshowthe bifurcation point ${\theta }_{0*}\approx 1.7$ (see Figure 7 and Figure 8).

Remark 1.

By making use of a suitable hybrid controller, we can narrow the stability region and advance the onset of Hopf bifurcation of the fractional-order duopoly game model (2). From an economic point of view, we can advance the cyclic state of the price of the tow insurance companies via adjusting the delay and feedback gain parameters.

8. Conclusions

The price of insurance companies plays an important role in dominating the market and attracting the consumers. The price competition of insurance companies is a vital topic. In this current research, we propose a new fractional-order duopoly game model with delays in insurance market. The existence and uniqueness, non-negativeness, boundedness of solution, stability, Hopf bifurcation, globally asymptotic stability, and Hopf bifurcation control of the involved fractional-order delayed duopoly game model in insurance market have been systematically explored. A series of sufficient conditions which guarantee the existence and uniqueness, non-negativeness, boundedness of solution, stability of the positive equilibrium, onset of Hopf bifurcation, and globally asymptotic stability of the addressed fractional-order delayed duopoly game model in insurance market, are derived. By virtue of hybrid control technique, we successfully control the stability domain and the time of generation of Hopf bifurcation of the involved fractional-order delayed duopoly game model in insurance market. The obtained study results own great theory value and praxis function use for reference in administering and running insurance companies. In addition, the research approach is also used to probe into bifurcation dynamics and its control issue of a number of other fractional-order systems appearing in many areas. Here we must point out that although we can control the stability region and the time of onset of Hopf bifurcation of the fractional-order duopoly game model via hybrid controller, we may not be able to other fractional-order delayed models via the same hybrid controller. We must take some adequate measures to control the stability region and the time of onset of Hopf bifurcation according to different fractional-order delayed models. We will address this aspect in the near future.