Bifurcations in economic growth model with distributed time delay transformed to ODE

We consider the model of economic growth with time delayed investment function. Assuming the investment is time distributed we can use the linear chain trick technique to transform delay differential equation system to equivalent system of ordinary differential system (ODE). The time delay parameter is a mean time delay of gamma distribution. We reduce the system with distribution delay to both three and four-dimensional ODEs. We study the Hopf bifurcation in these systems with respect to two parameters: the time delay parameter and the rate of growth parameter. We derive the results from the analytical as well as numerical investigations. From the former we obtain the sufficient criteria on the existence and stability of a limit cycle solution through the Hopf bifurcation. In numerical studies with the Dana and Malgrange investment function we found two Hopf bifurcations with respect to the rate growth parameter and detect the existence of stable long-period cycles in the economy. We find that depending on the time delay and adjustment speed parameters the range of admissible values of the rate of growth parameter breaks down into three intervals. First we have stable focus, then the limit cycle and again the stable solution with two Hopf bifurcations. Such behaviour appears for some middle interval of admissible range of values of the rate of growth parameter.


Introduction
In economics many processes depends on past events, so it is natural to use the time delay differential equations to model economic phenomena. Two main areas of applications are business cycle and economic growth theories. In last decades the study of impact of investment delay was a subject of an detailed scrutiny as a mechanism for endogenous cycle to describe business cycles and growth cycle. One of the most influential models of business cycle with time delay was the Kaldor-Kalecki model [1] which was based on the Kaldor model, one earliest endogenous business cycle model [2,3]. The Kaldor is a prototype of a dynamical system with cyclic behaviour in which nonlinearity plays a crucial role to generate the endogenous cycles. The nonlinearities are common feature used to model the complexity of economic systems [4]. In turn, the investment delay was assumed to be the average time of making investment as it was proposed by Kalecki [5]. Apart from the pioneering interest in time delay in economics, the time delay systems were studied in many domains of science [6].
The Kaldor-Kalecki business cycle model was a subject of many studies as well as augmentations. This model was modified by incorporating the exponential trend to describe growth of an economy [7]. This new Kaldor-Kalecki growth model was formulated as similar way as the Kaldor growth model was obtained from the Kaldor business cycle model [8].
The models with distributed delays are a more realistic description of economic systems with time delay. There are some examples of such models in context of economic growth [9]. In this paper we modify the Kaldor-Kalecki growth model to allow for distributed time of investment. Instead of average time of investment completion the distributed time length of investment is considered. This approach provides more realistic characteristic of investment processes. This growth model has the form of the dynamical system with a distributed time delay. The unimodal distribution function for investment is assumed.
While the delay differential equations methods are developing rapidly, the mathematical methods for ordinary differential equations are superlative, especially when distributed delays are considered. Therefore, it is convenient to approximate a system with distributed delay with a system of ordinary differential equations. One of the method to reduce a system with the distributed delay to an ordinary differential equation system is the linear chain trick [10,11,12]. In consequence infinite-dimensional dynamical systems are approximated by finite-dimensional dynamical system, where the dimension of the system can be chosen.
The main aim of the paper is to study possible bifurcation due to change of the parameter values of the Kaldor-Kalecki growth model. We consider two simplest cases of three and four dimensional dynamical systems obtained through the linear chain trick from the Kaldor-Kalecki growth model with distributed delay. For both models we establish the conditions for existence of the Hopf bifurcation with respect to the time delay parameter and the rate of growth parameter. We show the both parameters play role in a scenario leading to the Hopf bifurcation and arising cyclic behaviour.
In the numerical part of the paper we determine in details the ranges of parameter values for which cyclical behaviour is possible. For these analysis we use the investment function obtained by Dana and Malgrange for French economy [8]. As in theoretical part of the paper we choose the time delay parameter and the rate of growth parameter. Additionally we consider the adjustment parameter. We show that the combination of these three parameters of model can lead to arising the cycles through Hopf. In the space of these three parameters of the model we obtain the surface (a section of a paraboloid) separating the regions with stable and cyclic solutions.

Model
In Economic growth cycles driven by investment delay, Krawiec and Szyd lowski [7] formulated the model based on the Kaldor business cycle model with two modifications: exponential growth introduced by Dana and Malgrange [8] and Kaleckian investment time delay [1]. This model of economic growth is described by the following system of differential equations with time delay where I(y(t), k(t)) = k(t)Φ(y(t)/k(t)), with Φ(y/k) = c + d 1 + e −a(vy/k−1) , and α, γ, g, δ, G 0 , g and a, c, d, v are positive constants. It can be found that the system has a unique fixed point (y * , k * ) with positive coordinates, where , with x * the unique solution of the equation Because of the S-shape of function Φ(x), we have that x * always exists and the values of y * and k * depend only on x * (in our case, c < g + δ < c + d). Notice that, for economic considerations, the investment function I(y(t), k(t)) is such that and In this paper, we generalize their model by replacing the time delay in (2) with a distributed delay as follows, where κ(·) is a gamma distribution, i.e.
with m a positive integer that determines the shape of the weighting function. T ≥ 0 is a parameter associated with the mean time delay of the distribution. Notice that as T → 0 the distribution function approaches the Dirac distribution, and, thus, one recovers the time delay case. Henceforth, we will consider only the cases m = 1 (weak delay kernel) and m = 2 (strong delay kernel). Using the so-called linear chain trick technique [12], system (5)-(6) can be transformed into equivalent systems of ODEs. More precisely, defining the new variable one has the system (case m = 1) while defining the new variables We will now analyse the stability and Hopf bifurcation of systems (7)-(9) and (10)-(13) by determining eigenvalues of linearised systems around the critical point (y * , y * , k * ) and (y * , y * , y * , k * ), respectively.

The time delay bifurcation analysis
Case m = 1 The characteristic equation of the linearised system (7)-(9) at the critical point (y * , u * , k * ), where u * = y * , is given by where λ denotes a characteristic root. A direct calculation implies that (14) leads to λ 3 + a 1 (T )λ 2 + a 2 (T )λ + a 3 (T ) = 0, where The necessary and sufficient condition for the local stability of the equilibrium point is that all characteristic roots of (15) have negative real parts, which, from the Routh-Hurwitz condition, is equivalent to a 1 (T ) > 0, a 3 (T ) > 0 and a 1 (T )a 2 (T ) > a 3 (T ). Then, a 2 (T ) > 0 is necessarily satisfied. Let us examine whether these inequalities hold. First, we notice that A < 0. In fact, by contradiction, if A = 0, then a 2 (T ) = − x * I * y 2 < 0. On the other hand, if A > 0, then B > 0, and so a 2 (T ) < 0. The fact A < 0 implies that a 1 (T ) > 0 holds always true, while the inequality a 3 (T ) > 0 is valid if and only if B + αI * k I * y < 0. Thus, a 3 (T ) > 0 is always satisfied when B ≤ 0, and it is verified for g + δ − (g + αγ)x * < 0 when B > 0. Finally, let us consider a 1 (T )a 2 (T ) > a 3 (T ). Since the sign of a 1 (T )a 2 (T ) − a 3 (T ) depends on the sign of (AB)T 2 + (A 2 + αI * k I * y )T − A, which is a quadratic polynomial in T. We have now several cases.
ii) If B > 0, then AB < 0 and −A > 0. By Descartes' rule of signs we find that the polynomial (AB) iii) If B < 0, then AB > 0 and −A > 0. Applying again the Descartes' rule of signs we see that (AB)T 2 +(A 2 +αI * k I * y )T −A has (two) sign changes only if A 2 + αI * k I * y < 0, meaning that this polynomial may have two positive roots T * 2 < T * 3 . If this happens, then a 1 (T )a 2 (T ) − a 3 (T ) > 0 if 0 < T < T * 2 and T > T * 3 .
Let T = T * such that a 1 (T * )a 2 (T * ) − a 3 (T * ) = 0, namely T * = T * j (j = 0, 1, 2, 3). The curve T = T * divides the parameter space into stable and unstable parts. Choosing T as a bifurcation parameter, we apply the Hopf bifurcation theorem to establish the existence of a cyclical movement. This theorem asserts the existence of the closed orbit, if the characteristic equation (15) has a pair of pure imaginary roots and a non-zero real root, and if the real part of the imaginary roots is not stationary with respect to the changes of the parameter T . At the critical value T = T * , Eq. (15) factors as so we have the following three roots λ 1,2 = ±i a 2 (T * ) = ±iω * and λ 3 = −a 1 (T * ) < 0. Next, let us investigate the sign of the real parts of this roots as T varies. A differentiation of (15) with respect to T yields where Then, from (16), we get The previous analysis leads to the following conclusions.
2) If B < 0, then there exists 0 < T * 2 < T * 3 such that the equilibrium point (y * , y * , k * ) of (7)-(9) is locally asymptotically stable for all T < T * 2 and T < T * 3 , and unstable for all T * 2 < T < T * 3 . A comparison of 1/ √ −B with T * 2 and T * 3 yields that system (7)-(9) undergoes a Hopf bifurcation at (y * , y * , k * ) when T = T * 2 or T = T * 3 or T = T * 2 and T = T * 3 . For system (7)-(9) figure 1 presents the bifurcation diagram for the time delay parameter T . Case m = 2 The characteristic equation of the linearised system (10)-(13) at the critical point (y * , p * , w * , k * ), where p * = w * = y * , takes the form where I * y and I * k are defined as in (3) and (4), which leads to the following fourth order algebraic equation in λ where and (19) According to the Routh-Hurwitz conditions for stable roots, the equilibrium point (y * , y * , y * , k * ) of system (14) is locally asymptotically stable if Taking in mind that N < 0 and P > 0, we derive that conditions (3) is such that ϕ(0) < 0 and ϕ(+∞) = +∞. Hence, there exists at least a positive value of T, say T * 0 , such that ϕ(T ) = 0. We need now to recall Descartes' rule of signs and its corollary, that state "the number of positive roots of the polynomial ϕ(T ) is either equal to the number of sign differences between consecutive nonzero coefficients, or is less than it by an even number" and "the number of negative roots is the number of sign changes after multiplying the coefficients of odd-power terms by −1, or fewer than it by an even number", respectively. Applying these rules to the polynomial ϕ(T ), we get that ϕ(T ) has one positive zero and the number of negative zeros must be either 2 or 0. Therefore, ϕ(T ) < 0 if T < T * 0 . When M = 0, the equilibrium point (y * , y * , y * , k * ) of (14) is locally asymptotically stable for T < T * 0 . Assume there exists T * > 0 such that ϕ(T * ) = 0, i.e.
In this case, we can rewrite the characteristic equation (18) as so that we have two purely imaginary roots and two other roots, which have real parts different from zero since Differentiating the characteristic equation (18) with respect to T , we have where Letting λ = iω * in (21), a direct calculation yields Let us notice that sign Re (dλ/dT ) T =T * = sign [−ϕ (T * )], and recall that sign Re (dλ/dT ) T =T * > 0 and sign Re (dλ/dT ) T =T * < 0 correspond to crossings of the imaginary axis from right to left, and from left to right, respectively. Summarizing all the previous analysis, we have the following results.
Theorem 2 Let M be defined as in (19).
1) Let M = 0. There exists T * 0 > 0 such that the equilibrium point (y * , y * , y * , k * ) of (14) is locally asymptotically stable for T < T * 0 , unstable for T > T * 0 , and bifurcates to a limit cycle through a Hopf bifurcation at the equilibrium point when T = T * 0 .

The rate of growth bifurcation analysis
Let us consider the dynamics of the system (7)-(9) with respect to the change of the parameter g (the rate of economic growth).

Proposition 1
The critical point of system (7)-(9) (and equivalently system (1)-(2) always exists for the rate of growth parameter g in the interval c−δ < g < c + d − δ.
It is shown earlier that the system (1)-(2) has unique fixed point for c < g + δ < c + d. The economy with the investment function I(y, k) = kΦ(y, k) has a fixed point or a limit cycle solution only for some rates of growth within the interval (c − δ, c + d − δ). The parameters c and d of the investment function (and also capital stock depreciation) put some limit on minimal and maximal rates of growth. Let us consider the characteristic equation of the linearised system (10)-(13) at the critical point (y * , u * , k * ), in the form or λ 3 + a 1 (g)λ 2 + a 2 (g)λ + a 3 (g) = 0 (23) where a 1 (g) = 1 T − α(I * y (g) − γ) + g + x * (g)I * y (g) where λ is a root of the characteristic equation. The discriminant of the characteristic equation is Proposition 2 If the expression (24) is positive than the all eigenvalues are real and if it is negative there is one real, one pair of conjugate complex eigenvalues. For the zero value of the expression (24) the critical point is non-hyperbolic.
For real eigenvalues we have the following proposition Proposition 3 In the interval c−δ < g < c+d−δ, there are two subintervals with the positive values of the discriminant (24) there two cases for the values of rate of growth parameter g. In these subintervals there are three negative real eigenvalues.
The subintervals of the parameter g with two negative and one positive eigenvalues are non-physical regions as the critical point (y * , u * , k * ) does not lie in a positive quadrant.
For complex eigenvalues we have the following proposition Proposition 4 In the interval c − δ < g < c + d − δ and negative values of the discriminant (24) for the increasing value of the rate of growth parameter g there are two supercritical Hopf bifurcations. For the value g = g 1,Hopf the limit cycle is created, and for the value g = g 2,Hopf the limit cycle is destroyed (g 1,Hopf < g 2,Hopf ).
Therefore, as the the rate of growth parameter is increasing in the interval g min = c − δ < g < c + d − δ = g max the eigenvalues change as follows. In the first subinterval (g min ; g 1 ) there are three real eigenvalues (two negative, one positive). In the second subinterval (g 1 ; g 1,Hopf ) there are three real eigenvalues (three negative). In the third subinterval (g 1,Hopf ; g 2,Hopf ) there are one real eigenvalue (negative) and one conjugate complex eigenvalue (positive real parts). In the fourth subinterval (g 2,Hopf ; g 2 ) there are three real eigenvalues (three negative). And finally, in the fifth subinterval g 2 , (g max ) there are three real eigenvalues (two negative, one positive).
We some example values of parameters we can determine the values of the rate of growth parameter for which the eigenvalues change their character or sign. We assume the values of investment function parameters obtained by Dana and Malgrange, namely, c = 0.01, d = 0.026, a = 9, v = 4.23. We fix also the following model parameters α = 1, γ = 0.15, δ = 0.007, G 0 = 2 and T = 1. The rate of growth parameter g is taken within the interval g min = c − δ < g < c + d − δ = g max . The results are presented in table 1.
For system (7)-(9) figure 2 presents the bifurcation diagram for the rate of growth parameter g.

Numerical analysis of the Hopf bifurcation
The original Kaldor model exhibited the limit cycle behaviour due to the Hopf bifurcation caused by the increase of the parameter α value [3]. Later, it has been augmented both by introducing the investment lag T and exogenous growth trend g. The increase of the investment time delay parameter value also generates the limit cycle [1]. However, the dependence the Hopf bifurcation on the rate of growth parameter was not elaborated so far. Both Chang and Smyth in the Kaldor model [3] as well as Dana and Malgrange in the Kaldor model with exogenous growth trend investigated the parameter α as the bifurcation parameter [8].

Case m = 1
In this case we consider the three-dimensional model (7)-(9) for the state variables (y, u, k). In this model we study numerically the stability of the critical point (y * = u * , k * ) to find the values of parameter T for which the critical point loses the stability and the limit cycle is created through the Hopf bifurcation mechanism. We study in details the dependence of the bifurcation value of T on the model parameters α and g as well as the dependence the bifurcation value of g on the model parameter α and T . The bifurcation surface in the parameter space (a, g, T ) is presented in figure 3. The region below the surface corresponds to the asymptotic stability of the critical point (y * , u * , k * ). The region inside corresponds to parameter values for which system (4) has an unstable critical point with a limit cycle around it.
Let conduct more detailed analysis and consider relations between two parameters with a third parameter fixed. First, the relation of T on α for g = 0.016 is shown in Fig. 4. We find that the asymptotic stability region exist only if α < 0.7644 (with g = 0.016). In the interval of α ∈ (0.6, 0.764) we find the relation T bi (α) is Second, we analyze the dependence of the parameter T on the parameter g with the fixed value of α. We consider the three values of parameter α. The stability regions on the plane (g, T ) are shown for α = 0.6 and α = 0.9 in fig. 5. As the value of the parameter g increases the region B is growing.
The bifurcation line separating the region A and B is described by a quadratic equation T bi = a 2 g 2 + a 1 g + a 0 .  Figure 4: The plane of parameters (α, T ) for system m = 1 and g = 0.016 with investment function (25). Region I is the region of asymptotic stability, while region II is the region of parameters value for which a limit cycle solution exists.
an oscillating manner. It is the region of asymptotic stability of the model. In fig. 6 there are the two cases for α = 0.6, 0.9 with three solutions y(t) obtained for given T = 0.5, 1.5, 3 and the same initial function y(t) = 1, k(t) = 100 where t ∈ (−T, 0). We can see that the dumping of oscillations is weaker as the time delay T increases.

Case m = 2
In this case we consider the four-dimensional model (10)-(13) for the state variables (y, p, w, k). In this model the critical point values of (y * = p * = w * , k * ) are the same as the critical point values of (y * = u * , k * ) of three dimensional model presented in the previous section.
In a similar way as in the previous section we analyse the occurrence of the Hopf bifurcation for the parameter T depending on parameters α and g. The relation of T on α for g = 0.016 is shown in Fig. 7. The asymptotic stability region exist only if α < 0.7644 (with g = 0.016). It is the same result as in the case of model m = 1. Let us study the cycle characteristics for some values of the parameter α with different values of the delay parameter T . Figure 8 shows the solutions of y for two cases of α = 0.6, 0.9 and three values of the parameter T = 0.5, 1.5, 3.0. The amplitude and period of cycles are decreasing as the parameter T increases for α = 0.6.
Comparing m = 1 and m = 2 The models considered can be treated as the approximation of the delay Kaldor-Kalecki growth model. It is important to find how good the subsequent approximations are. Therefore, we compare the bifurcation values of the parameter T in models m = 1 and m = 2.
First, we consider the bifurcation diagram in the parameter plane (α, T ) presented in Fig. 9. We find that for the given value of parameter α the Hopf bifurcation value of the parameter T bi is lower for the model m = 2. This difference is zero for T = 0 and then increases as the value T bi increases for We compare trajectories of y(t) for systems m = 1 and m = 2 with the same initial conditions. In Figs. 10, 11 there are trajectories y(t) for assumed the parameter g = 0.016 and combinations of parameters α and T . We find that for the same parameters α, g and T the period of cycles is smaller for model m = 2 and amplitude is also smaller, although the difference is very small. For example, for α = 0.9, g = 0.016 and T = 3, the period of cycle in model m = 1 is 114.85 and in model m = 2 is 116.45 while amplitudes are 12.9555 and 12.966, respectively. Taking models with greater m we should obtain the cycles with longer periods and amplitudes.
We can explore further approximations of the model and compare values of the bifurcation parameter g for successive m in table 2.

Conclusions
We study the Kaldor-Kalecki growth model with distributed delay transform to the ordinary differential equation system using the linear chain trick. It allows to choose the arbitrary dimension of the system. We consider two cases where the system is transformed to three-dimensional and four-dimensional dynamical systems. The model possesses the equilibrium for some interval of values of rate of growth parameter. The interval depends on the values of investment function which is taken to model analysis.
We study a bifurcation to a limit cycle (the Hopf bifurcation) due to the change of two parameters: the time delay T and the rate of growth g. Additionally we take take into consideration the speed of adjustment parameter α which is the bifurcation parameter of original Kaldor model. For the better insight of dynamics we analyzed the bifurcations in the the three dimensional space of these parameters.
When we consider the dynamics of the model under the change of the growth rate parameter, we discover numerically two bifurcation values of the rate of growth parameter when the Hopf bifurcations occur. Increasing the value of the rate of growth parameter, for a smaller value of this parameter the limit cycle emerges then for a larger value of this parameter the limit cycle is destroyed to a stable focus through the Hopf bifurcations. Therefore, the cyclic behaviour takes place in some interval of the rate of growth parameter values. Outside of the interval the systems through damping oscillation goes to a stable stationary solution. The similar dynamic behaviour with two Hopf bifurcations separating stable, unstable and stable regions was also found in the macroeconomic model extending the Calvo and Obstfeld framework [13].  All numerical analyses have been done with Dana and Malgrange's investment function for the French macroeconomic data [8].
• The Kaldor-Kalecki model with distributed delay is reduced to the ordinary differential system using the linear chain trick technique.
• Depending on the value of the parameter m of the Γ distribution function the reduced system is (m + 2)-dimensional ordinary differential equation system. • For the increasing time delay parameter there is the supercritical Hopf bifurcation.
• For the increasing rate of growth parameter, first the limit cycle emerges and then the limit cycle disappears. Therefore, there are two supercritical Hopf bifurcations with two bifurcation values of the rate of growth parameter.
• For some values of parameters α and T , in the allowed range of the rate of growth parameter values, both for lower and higher values of the rate growth parameter the model has the stable stationary point while for the middle range of parameter values there is the limit cycle.
• The period of cycle increases and decreases as the rate of growth parameter increases in the range of unstable solution.
• Comparing the models with different m the stable region in the parameter space is slightly diminished as m is greater.