1. Introduction
The examination of economic growth models is one of the most frequently discussed issues in mathematical economics. Day [
1,
2] has investigated a neoclassical growth model, and a productivity and population growth model and showed the emergence of complex behavior even under simple economic structure. His models were based on discrete time scales and a mound-shaped production function that represented the negative effect of pollution resulting from increasing capital. It was demonstrated by numerical computations that these models could generate cyclic and even chaotic behavior. Following Day’s pioneering works, a lot of effort has been given to the understanding of complex economic dynamics. Day [
3], Puu [
4] and Bischi
et al. [
5] present the earlier contributions of this field. A large number of studies assumed discrete time scales. Li and Yorke [
6] have introduced the “period-three condition” to detect chaos, which has many applications in first-order nonlinear difference equations. The papers collected by Rosser [
7] offer many applications. Only a few studies are devoted to the case of continuous time scales, since there is no general criterion to detect chaos and the system must have at least three dimensions.
In this paper, we will examine an extension of the neoclassical growth model, which can be traced back to the early works of Solow [
8] and Swan [
9]. The neoclassical growth model is constructed with the two (usually implicitly mentioned) assumptions; one is the full-employment of labor and capital and the other is instantaneous adjustment in the output market. Thus, it is suitable for describing the long-run behavior of the economy. Due to the well-behaved production function, the steady state of the model is usually asymptotically stable. However, it is often observed in reality that growth path exhibits persistent fluctuations. The neoclassical model could be the good point of departure to show how such persistent behavior can emerge when nonlinearities and a production lag are present. Matsumoto and Szidarovszky [
10] attempt to fill the gap and have introduced a neoclassical model with a mound-shaped production function that was assumed to be a Cobb-Douglas type function of the form
F(
x)
=Axa(1 -
x)
b with
x being the capital per unit labor. Although they show emergence of erratic fluctuations in the capital accumulation process, the production function is restrictive in the sense that it is defined only in the unit interval. This paper modifies this drawback, considers another type of mound-shaped production function and will examine the stability of the steady state with and without time delays in the continuous-time framework. Two kinds of delays will be discussed, fixed and continuously distributed (continuously hereafter) delays. We keep the relatively simple model of Matsumoto and Szidarovszky [
10] in order to be able to compare the results and to illustrate that complex dynamics can be generated under simple economic assumptions with both function types.
This paper develops as follows. First the mathematical model is formulated without time delays, and complete steady state and stability analysis is presented. Then, models with fixed delays and then with continuous delays are introduced and complete stability analysis is given. The last section concludes the paper.
2. The Mathematical Model
Matsumoto and Szidarovszky [
10] have introduced a special growth model of the form
where
x is the capital per labor,
s and
α are positive parameters where
s ∈ (0,1) is the average propensity to save and
α = n+sμ with
μ being the depreciation ratio of capital and
n the growth rate of labor. In applying the Cobb-Douglas type function
F(
x) =
Axα (1 -
x)
b the value of
x has to be normalized into the (0,1) interval. In this paper we will assume that
which is a mound-shaped function in which
x can take any positive value. This function has zero value at
x = 0, converges to zero as
x → ∞, increases for
x <
γ /
δ and decreases for
x >
γ /
δ, so it has its maximum at
x =
γ /
δ. The drawbacks of the neoclassical production function that ignores natural resource or energy are partly remedied in our function.
δ of the exponential term reflecting a strength of a “negative effect” caused by increasing concentration of capital. The value of
δ is exogenously given, however, can be thought to be determined by a damaging degree of natural environment or energy resources. With the modified function, the mathematical model becomes
where
α, γ, δ and
β =
s ε are positive parameters. The number of steady states and their locations depend on the specific values of the model parameters. We will consider three different cases:
γ < 1,
γ = 1 and
γ > 1. We can give the following interpretation for the value of the parameter γ of the production function; γ can be thought as a proxy for measuring returns to scale of the production function. Indeed, when
x is small (
i.e., the exponential term is close to unity), output increases more than unity, exactly unity and less than unity if
γ > 1,
γ = 1 and
γ < 1, respectively. Let now
f(
x) denote the right hand side of equation (1). If
x(0) = 0, then the identically zero function is a solution which case is not interesting from the economic point of view and is eliminated from further considerations
Case I.
Assume first that
γ < 1. The steady states are the solutions of
f(
x) = 0. Notice that
f(0) = 0, so zero is a steady state.
f(
x) converges to - ∞ as
x → ∞ Since
f´(
x) converges to ∞ as
x tends to zero with positive values. Hence
f(
x) increases for small values of
x > 0. The steady state equation
f(
x) = 0 can be written as
so the positive steady state is the unique solution of equation
The left hand side is zero at
x = 0 and strictly increasing, furthermore, converges to ∞ as
x tends to infinity. The right hand side is
β > 0 at
x = 0, strictly decreases and converges to zero as
x → ∞. Hence there is a unique solution
> 0 of (3), and
f(
x) > 0 if
x <
and
f(
x) < 0 as
x >
These relations imply that if
x(0) <
, then
x (
t) increases and converges to
, and if
x (0) >
, then
x (
t) decreases and converges to
. If
x(0) =
, then
x (
t) remains
for all
x > 0. Thus
is globally asymptotically stable.
Case II.
Assume next that
γ = 1. Then the steady state equation has the form
so zero is a steady state and there is a unique root of the second factor,
If
β ≤
α, then the value of
f(
x) is negative for all
x > 0. Therefore
x(
t) is decreasing and converges to zero with arbitrary
x (0) > 0. If
β >
α, then
> 0, furthermore
f(
x) > 0 as
x <
, and
f(
x) < 0 as
x >
. If
x (0) <
, then
x (
t) increases and if
x (0) >
, then
x decreases and converges to
, and if
x (0) =
, then
x (
t) remains
for all
T > 0 Hence
is globally asymptotically stable.
Case III.
Consider finally the case of
γ > 1. The steady state equation has now the form
so zero is a steady state again, and any other steady state is the solution of equation
Notice that
and
Therefore
g(
x) as its global maximum at
increases for
x <
and decreases for
x >
. Now we have three sub-cases.
(i) if
g(
) < 0, then there is no positive steady state and with arbitrary
x(0)>0,
x(
t) decreases and converges to zero.
(ii) if
g(
) < 0, then
=
is the unique positive steady state and
f(
x) < 0 for all 0 <
x ≠
. If
x(0) <
, then
x (
t) decreases and converges to 0, and if
x(0) >
then
x (
t) decreases again and now converges to
. If
x(0) =
, then
x (
t) =
for all
t > 0.
(iii) If
g(
) > 0, then equation (7) has two positive solutions,
1 <
and
2 >
. Relation (6) implies that
f(
x) < 0 as
x <
1 or
x >
2, and
f(
x) > 0 if
1 <
x <
2. Therefore if
x(0) <
1, then
x (
t) decreases and converges to zero, if
1 <
x(0) <
2, then
x (
t) increases and converges to
2, if
x (0) >
2, then
x (
t) decreases and converges to
2. That is,
1 is locally unstable and
2 is locally asymptotically stable. If
x(0) =
1 or
x (0) =
2, then
x (
t) remains at that steady state level for all
t > 0.
3. Model with Fixed Delay
The fixed delay
T > 0 is assumed in the second term of the right hand side of equation (1), so we have the following equation:
Where
The local asymptotic behavior of the trajectory can be examined by linearization. Let
be a positive steady state. Then the linearized equation has the form
where
xδ(
t) is the deviation of
x (
t) from the steady state level. Looking for the solution in the usual form
xδ(
t) =
eλtu, we have
which gives the characteristic equation
or
Lemma 1 Assume that |
h´(
)| <
α. Then is locally asymptotically stable.
Proof. Assume that Re
λ ≥ 0. Then
and since
clearly
Therefore
implying that
λ cannot be an eigenvalue. Q.E.D.
Notice that
and at the steady state
implying that
Therefore
so the characteristic equation (11) can be rewritten as
We also mention that the condition of Lemma 1 can be rewritten as
or equivalently
In the special case of
γ = 1, this condition has the form
which is equivalent to relation
In order to give a complete stability analysis, we have to find the possible stability switches. Substituting any stability switch,
λ=
iω with
ω > 0 into equation (11) yields
Separating the real and imaginary parts gives two equations,
and
Adding the squares of these equations gives
so
In order to have solution we have to assume now that
that is, (14) is violated with strict inequalities. Concerning this assumption, we can give the following interpretation. Let
F(x) be
y. Then
So this means that the absolute value of the elasticity of output with respect to capital is larger than unity. From (16) we have that if
h΄(
)>0, then
T(n) =
T+(n) with
and if if
h´(
) < 0, then
T (n) =
T-(n) with
and by (13) and (19),
so
h´(
) cannot be zero.
By selecting
T as the bifurcation parameter and implicitly differentiating the characteristic equation with respect to
T, we have
implying that
If
λ =
iω, then
with real part
Therefore if a steady state is unstable with
T = 0, then it remains unstable for all
T > 0, and if a steady state is asymptotically stable at
T = 0, then this stability is lost at
T =
T(0) and cannot be regained later. That is, if the steady state is unstable without delay, then it remains unstable with any delay of positive length. If the steady state is asymptotically stable without delay, then it remains asymptotically stable until the delay reaches a certain threshold, and then becomes unstable and the stability cannot be regained later.
Taking,
α = 1,
β = 25,
γ = 1 and
δ = 1, we give an illustrative numerical example in
Figure 1. The critical value
γ -
δ is denoted by
Introducing the notation
z =
γ -
δ transforms the
T -(0) curve to
and then the corresponding critical value of the delay is
In
Figure 1(A), the steady state is locally asymptotically stable in the dark-gray region with
z > -1 due to Lemma 1. It is also locally asymptotically stable in the light-gray region, which is under the critical curve
T =
T-(0) and it is unstable in the white region above the curve. Setting
z =
zc and increasing
T along the vertical dotted line in
Figure 1(A), we can see that the steady state loses stability at
T =
Tc. Further increasing
T , as observed in
Figure 1(B), generates complex dynamics through a
quasi period-doubling bifurcation in which
T increases from
T c- 0.05 to 8.5 with an increment of 0.01 and the local maximum and minimum of the corresponding trajectory are plotted against each value of
T.
Figure 1.
Dynamics with α = 1, β = 25, γ = 1 and δ = 1.
Figure 1.
Dynamics with α = 1, β = 25, γ = 1 and δ = 1.
4. Model with Continuously Distributed Delay
Assuming continuously distributed delays in the second term of equation (1) gives the following Volterra-type integro-differential equation:
where
T > 0 is a positive parameter, the average delay and
m ≥ 0 is an integer. The kernel function has the form
If
m = 0, then the most current value has the highest weight, which is then decreasing exponentially. If
m ≥ 1, then the most current value has zero weight which is then increasing until
t -
s =
T, and is decreasing exponentially afterwards. As
m increases, the weighting function becomes more peaked around
t -
s =
T and as
m → ∞, it converges to the Dirac delta function centered at
t -
s =
T. If
T → 0, then it also converges to the Dirac delta function.
A similar model is investigated by Fanti and Manfredi [
11] where
m = 2 is selected and the stability of the system with a cubic characteristic polynomial is examined based on the Routh-Hurwitz criterion. Stamova and Stamov [
12] consider a generalized Solow model with endegenous labor growth and impulsive perturbations. Their stability analysis is based on the Lyapunov-Razumikhin sufficient stability conditions, which is a different approach than ours.
Linearizing equation (22), we have
where
xδ(
t) is the deviation of
x (
t) from the steady state level
. We are looking for the solution in the usual exponential form
then simple substitution shows that
Notice that by introducing the length of the delay as the new variable
S =
t -
s in the integral, we see that
and by letting
T → ∞, we have the characteristic equation
with
This equation can be rewritten as
Then similarly to the case of fixed delay we can prove the following result:
Lemma 2 Assume that |h´(
)| < α. Then
is locally asymptotically stable.
Proof. Assume that
Reλ≥ 0. Then
therefore
implying that
λ cannot be a solution of equation (23). Q.E.D.
It is well known that the Routh-Hurwitz stability theorem provides necessary and sufficient conditions for all roots of a polynomial equation with real coefficients to have negative real parts. It is also known that it is difficult to locate the eigenvalues with analytic methods in general. However in some special cases, analytic results are still possible to obtain, as it will be next demonstrated.
Case I. T = 0
Assume first that T = 0, which reduces equation (22) with delays to equation (1) without delays. The asymptotic properties of this equation were already discussed earlier.
Case II. T > 0 and m = 0
Assume next that
T > 0 and
m = 0, when the kernel function becomes exponentially declining. Then characteristic equation (23) becomes quadratic,
or
If
γ ≤ 1, then all coefficients are positive with a positive steady state, which is locally asymptotically stable. Assume next that
γ > 1. If
then the constant term is zero indicating that one eigenvalue is zero and the other is negative. So
is marginally stable in the linearized model, so no conclusion can be drawn about its asymptotical behavior in the nonlinear model. If
< (
γ - 1)/
δ, then
is unstable and if
> (
γ - 1)/
δ, then
is locally asymptotically stable.
Case III. T > 0 and m = 1
Assume now that
T > 0 and
m = 1, when the shape of the kernel function takes a bell-shaped form. Then we have a cubic characteristic polynomial:
or
If
γ ≤ 1, then all coefficients are positive at a positive steady state. If
γ > 1, then we can consider three cases. If
then zero is an eigenvalue and the other two eigenvalues have negative real parts implying that
in the linearized system is marginally stable. Therefore no conclusion can be drawn about the stability of
in the nonlinear system. If
< (
γ - 1)/
δ, then the constant term is negative, so
is unstable. If
> (
γ - 1)/
δ then all coefficient of (25) are positive. In this case and when
γ ≤ 1 the Routh-Hurwitz criterion implies that the real parts of the eigenvalues are negative if and only if
which can be reduced to a quadratic inequality in
T:
For the sake of simplicity, we re-introduce the notation =
γ -
δ . If
z ≥ 0, then this inequality holds implying the asymptotical stability of the steady state.So we can assume that
z < 0.The discriminant of the left hand side of inequality (26) is
If
z < -8, then
D > 0, so the left hand side of (26) has two roots
which are positive and
<
. Notice that
= 1/
α2 and (26) holds if and only if
T <
or
T >
when
is locally asymptotically stable. If
<
T < , then (26) is violated, so
is unstable. If
z = -8, then
D = 0 and there are equal roots
so if
T ≠ 1/
α, then
is locally asymptotically stable.
If -8 <
z < 0, then
D < 0, so (26) holds and
is asymptotically stable. The instability region is shown in
Figure 2(A) where
z is the horizontal axis and
T is the vertical axis. If we start with a very small value of
T with any given
z < -8, then
is asymptotically stable. If we gradually increase
T then
remains asymptotically stable until it reaches the critical value
, when the steady state becomes unstable. It remains unstable until
when stability is regained, and
remains asymptotically stable for all
T >
.
If the steady state is unstable without delay, then it remains unstable with continuous delay with any T and m. If it is asymptotically stable without delay, then either it remains asymptotically stable with all T and m, or loses stability at a certain value of the average delay T and stability is regained with an even larger value of T and the steady state remains asymptotically stable afterwards. In such cases small and large average delays lead to asymptotically stable steady states.
We will next show that at the critical values
and
, Hopf bifurcation occurs giving the possibility of the birth of limit cycles. We select
T as the bifurcation parameter. At the critical values (26) is satisfied with equality, so
and the characteristic equation (25) can be rewritten as
showing that there is a negative eigenvalue
and a pair of pure complex eigenvalues
Consider
λ as a function of the bifurcation parameter
T and differentiate implicitly equation (25) to have
By simple calculation we can see that at
λ = ±
iε,
with real part
Since
= 1/
α2, at
T =
the value of
d(
Reλ)/
dT changes from negative to positive showing the loss of stability, and if
T =
, then
d(
Reλ)/
dT changes from positive to negative indicating that stability is regained. Since at both critical values
d(
Reλ)/
dT ≠ 0, at both values Hopf bifurcation occurs giving the possibility of the birth of limit cycles.
We perform numerical simulations to illustrate the results obtained above. In
Figure 2(A), the steady state is locally asymptotically stable in the dark-gray region with
z > -1 due to Lemma 2. It is also locally asymptotically stable in the light gray region and unstable in the white region when
z < - 1.The appearance and disappearance of a limit cycle can be observed in
Figure 2(B) where we take
α = 1,
β =
e13,
γ = 1 and
δ = 1 implying
z = - 12,
Under these specifications, the Volterra-type integro-differential equation (22) can be written as a 3D system of differential equations,
where
and
When
T increases from
- 0.1 to
+ 0.3 with an increment of 0.01, the steady state loses stability at point
A and regains stability at point
B In
Figure 2(B), the local maximum and local minimum of a trajectory generated by the 3D system are depicted against each value of
T indicating the birth of a limit cycle for
<
T < .
Figure 2.
Dynamics with T > 0 and m = 1.
Figure 2.
Dynamics with T > 0 and m = 1.
The cases of
m ≥ 2 result in fourth or larger degree polynomial equations. The stability of the steady states can be examined similarly, but the mathematical details become much more complicated. It can be mathematically confirmed that as
m → ∞, equation (23) converges to the characteristic equation (11) of the model with fixed delay. In particular, if
m → ∞, then expression
converges to
eλT. For larger values of
m dynamics generated by the differential equation with continuously distributed time delay is similar to dynamics generated by the differential equation with fixed time delay.