Note on Goodwin’s 1951 Nonlinear
Accelerator Model with an Investment Lag∗
Akio Matsumoto
Department of Economics, Chuo University742-1, Higashinakano, Hachioji
Tokyo, 211-0035, JapanMay 2, 2007

Abstract
It has been well-known that nonlinearity, time delay and local instability are significant sources for a birth of cyclical dynamics since
the pioneering work of Goodwin (1951). He has constructed a delay
nonlinear business cycle model with the nonlninear acceleration principle which gives rise to cyclic oscillations when a stationary state islocally unstable. However very little is known about time delay eﬀects
caused by investment lags in Goodwin’s cyclic dynamics, furthermore
global dynamics in the locally stable case has not been considered yet.
This study draws attentions to these unexplored aspects of Goodwin’s
business model. It is shown that increasing investment lag makes thelength of business cycle longer and its amplitude larger. It is further
demonstrated that multiple cycles coexist when the stationary state islocally stable. These results imply two issues: Goodwin’s model is not
only robust in cyclic oscillations regardless of local dynamic properties
but also corridor stable when the stationary state is locally stable.
Key words: fixed time delay, continuously distributed time delay,
corridor stability, coexistence of multiple limit cycles.
∗

The author appreciates financial supports from the Japan Ministry of Education,
Culture, Sports, Science and Technology (Grand-in-Aid for Scientifit Reseach A, 16203019)
and from Chuo University (Joint Reseach Project 0382). The paper is benefited fromcomments by Ferenc Szidarovszky and Mami Suzuki. The usual disclaimer applies.

11

Introduction

The contributions of Goodwin (1951) are reconsidered and further developed in this study. Goodwin introduced a nonlinear accelerator businesscycle model with an investment lag, numerically specified it and graphicallyshowed that it could generate a stable limit cycle when a stationary point islocally unstable. Since Goodwin’s work, it is expected that instability, nonlinearity and delay could be significant sources for the birth of cyclic behavior.
In view of the fact that it is diﬃcult to analytically solve delay nonlinear
models, it is a natural way to perform numerical studies or to convert themodel to a tractable one by using approximation. Indeed, considerable effort has been devoted to investigate the nonlinear structure of the ordinary
diﬀerential version of the unstable Goodwin’s model. Recently, Sasakura(1996) gives an elegant proof of the stability and the uniqueness of Goodwin’s cycle. More recently Lorenz and Nusse (2002), based on Lorenz (1987),
reconstructs Goodwin’s model as a forced oscillator system and demonstrates
the emergence of chaos when nonlinearities become stronger. In the existing
literature, however, there have been only limited analytical works on the delay diﬀerential version,1 and, furthermore, very little has yet been revealedwith respect to the circumstances under which the stationary point is locally asymptotically stable. The main purpose of this study is to provide aninvestigation of these unexplored aspects of Goodwin’s business cycle model.
We add two new observations to the existing results. First, we reformate
the model in terms of a nonlinear diﬀerential equation with the explicit treatment of time delay (i.e., fixed time delay and distributed time delay) and findout how the time delay aﬀects the emergence and characteristics of a limit
cycle. Second, we demonstrate that the nonlinear delay Goodwin’s modelhas the corridor stability when a stationary point is locally stable. Thatis, the model is stable and the trajectory returns to the stationary state forsmaller disturbances but is unstable and exhibits persistent fluctuations forlarger disturbances.
In what follows, Section 2 overviews the basic structure of Goodwin’s
nonlinear accelerator model and introduces time delay to see its eﬀect on
cyclic dynamics. Section 3 shows a coexistence of a stable stationary point,
an unstable limit cycle and a stable limit cycle. Section 4 concludes thepaper.
1

Yoshida and Asada (2007) investigates the impact of delayed government stabilization
policy on the dynamic behavior of a Keynes-Goodwin model. Also see Bischi, Chiarella,
Kopel and Szidarovszky (2007) for applications of the delay diﬀerential method to oligopoly
models.

2

2

Goodwin’s Business Cycle Model

This section is divided into three parts. In Section 2.1, we recapitulate thebasic elements of Goodwin’s model and perform simulations to see the birthof limit cycles when a stationary state is locally unstable and the investment
function is nonlinear. We, then, adopt an explicit treatment of the investment
lag into the model. In particular, we will examine fixed time delay in Section
2.2 and continuously distributed time delay in Section 2.3, and we will show
see how such a delay aﬀects the characteristics of cyclic dynamics.

2.1

Basic Model

Goodwin (1951) presents five diﬀerent versions of the nonlinear accelerator
model. The first version assumes a piecewise linear function with three levels
of investment, which can be thought as the crudest or simplest version of
the non-linear accelerator. This is a text-book model that can give a simple
exhibition on how nonlinearities give rise to endogenous cycles without relying on structurally unstable parameters, exogenous shocks, etc. The second
version replaces the piecewise linear investment function with a smooth nonlinear investment function. Although persistent cyclical oscillations of output
are shown to exist, the second version includes a unfavorable phenomenon,
namely, discontinuous investment jump , which is not realistic in the real
economic world. "In order to come close to reality" (p.11, Goodwin (1951)),
the investment lag is introduced in the third version. However, no analyticalconsiderations are given to this third version. The existence of a businesscycle is confirmed in the fourth version, which is a linear approximation of
the third version with respect to the investment lag. Finally alternation ofautonomous expenditure over time is taken into account in the fifth version.
To find out how nonlinearity works to generate endogenous cycles, we
review the second version of Goodwin’s model, which we call the basic model,⎧⎨ εẏ(t) = k̇(t) − (1 − α)y(t),(1)⎩k̇(t) = ϕ(ẏ(t)).

Here k is capital stock, y national income, α the marginal propensity to consume, which is positive and less than unity, and the reciprocal of ε a positive
adjustment coeﬃcient. The dot over variables stands for time diﬀerentiation. The first equation of (1) defines an adjustment process of the national
income. Accordingly, national income rises or falls if investment is larger
or smaller than savings. The second equation, in which ϕ(ẏ(t)) denotes theinduced investment, describes an accumulation process of capital stock based3

on the acceleration principle. According to the principle, investment depends
on the rate of changes in the national income. A distinctive feature of Goodwin’ model is to introduce a nonlinearity into the investment function in sucha way that the investment is proportional to the change in the national income in the neighborhood of the equilibrium income but becomes inflexible(i.e., less elastic) for extremely larger or smaller values of the income. Thisnonlinear acceleration principle is crucial in obtaining endogenous cycles in
Goodwin’s model. We will next retain this nonlinear assumption and specify
its explicit form. On the other hand, we depart from Goodwin’s non-essential
assumption of positive autonomous expenditure and will work with zero autonomous expenditure for the sake of simplicity. A direct consequence of thisassumption is that an equilibrium solution or a stationary point of the basicmodel is y(t) = ẏ(t) = 0 for all t.
Inserting the second equation of (1) into the first one and movining theterms on the left hand side to the right gives a single dynamics equation for
the national income y,εẏ(t) − ϕ(ẏ(t)) + (1 − α)y(t) = 0.

(2)This is a nonlinear diﬀerential equation. Although it is one-dimensional,
its nonlinearity prevents deriving an explicit form of the solution. In spiteof this simple form, it is possible to detect local dynamics by examining
its linearized version in a neighborhood of the stationary point and global
dynamics by performing numerical simulations.
The linear version of the income dynamic equation (2) isεẏ(t) − ν ẏ(t) + (1 − α)y(t) = 0,(3)

where ν = ϕ0 (0) is the slope of the investment function at the stationary
point. This is a first-order ordinary diﬀerential equation. Applying separation of variables gives a complete solution,y(t) = y0 eλt with λ =

1−α
,
ν −ε

(4)

where y0 is an initial condition. The stationary point is stable or unstableaccording to whether the eigenvalue λ is negative or positive. Since 1 − α isthe positive marginal propensity to save, the sign of the eigenvalue depends
on whether the numerator is positive or negative. Thus the stationary point
is locally stable if ν < ε and unstable if ν > ε.
We turn our attentions to global dynamics under the assumption of local
instability, namely, where ν > ε . We first specify the investment function
as well as the values of the coeﬃcients of (2) and then perform simulations4

to see what dynamics of y can be generated by the nonlinear equation (2).
Although Goodwin (1951) assumed piecewise linear investment function, we,for the sake of analytical convenience, adopt a smooth nonlinear investment
function of the form of an arctangent,©
ªϕ(ẏ(t)) = δ tan−1 (ẏ(t) − a) − tan−1 (−a) , δ > 0 and a > 0.
(5)This function has endogenous "ceiling" and "floor" and is asymmetric whenforthe parameter, a, is non-zero. In what follows, we set a = 1 and δ = 12πwhich the investment function (5) passes through the origin, and its ceilingis three time higher than its floor as it was the case in Goodwin’s model.
In numerical examples below, the initial point is set at y(0) = 10, andwe set of ε = 0.5 as the adjustment coeﬃcient and α = 0.6 as the marginal
propensity to consume as in Goodwin’s model. The numerical results are
given in Figure 1 in which an endogenous cycle is illustrated on the right andthe corresponding time path of output on the left. Along the upside-down
N-shaped locus on the right, the initial point, denoted by I0 , is displacedslightly upward to point A so that the output is increasing to yH , the highest level. Investment immediately switches discontinuously from positive tonegative. Consequently, the orbit jumps from point A to point B. With negative ẏ(t) at point B, the national income gradually declines from point B
to point C so that the output is decreasing to yL , the lowest level. Once
point C is reached, investment switches again discontinuously from negative
to positive. In other words, the orbit jumps again to point D from point C,from which the national income glides toward point A, and then the processrepeats itself. Thus we have a closed orbit constituting a self-sustaining cycle.
The points A and C are critical points at which one of the variables makes a
discontinuous jump. The same dynamics are depicted as a function of time ton the left in Figure 1. The kinked points of the time-trajectory correspondto the discontinuous jumps of the cyclic behavior. This is a simple exhibition of generating an endogenous cycle of output. We can summarize thesenumerical results as follows.
Result 1 Given ν > ε (i.e., locally unstable stationary point), a slow-rapid
limit cycle is brought about in the basic model due to the nonlinear
investment accelerator (5).

Insert Figure 1 Here.

5

2.2

Delay Model with Fixed Time Lags

Due to the fact that in real economy, plans and their realizations need time totake eﬀects, Goodwin (1951) introduces the investment lag, θ, between decisions to invest and the corresponding outlays. Inserting θ into the investment
function of the basic model yields the third version,
εẏ(t) − ϕ(ẏ(t − θ)) + (1 − α)y(t) = 0.

(6)This is a neutral delayed nonlinear diﬀerential equation, which we call thefixed delay model. Goodwin does not analyze dynamics generated by this
fixed delay model. Furthermore, to the best of our knowledge, no analyticalsolutions of the delayed model are available yet. However, it is possible,
again, to investigate dynamics of the delayed model by using linearization
for local dynamics and numerical simulations for global dynamics. Since a
cyclic oscillation has been shown to exist in the basic model, our main concern
is to see how the presence of the investment lag aﬀects characteristics of sucha slow-rapid cycle. To this end, we analytically investigate the stability of thecycle generated in the linearized model and numerically detect what eﬀects
are caused by the lag on cyclical dynamics.
The fixed delay model is autonomous and its special solution is constant(i.e., y(t) = 0) so that its linearized version takes the form of a linear neutral
autonomous delay diﬀerential equation,εẏ(t) − vẏ(t − θ) + (1 − α)y(t) = 0.

(7)

It is well known that if the characteristic polynomial of a linear neutral
equation has roots only with negative real parts, the stationary point islocally asymptotically stable. The normal procedure for solving this equation
is to try an exponential form of the solution. Substituting y(t) = y0 eλt
into (7) and rearranging terms, we obtain the corresponding characteristic
equation:
ελ − vλe−λθ + (1 − α) = 0.
To check stability, we determine conditions under which all roots of thischaracteristic equation lie in the left or right of the complex plane. Dividing both sides of the characteristic equation by ε and introducing the newvariables A = 1−αand B = − νε , we rewrite the characteristic equation asε
λ + A + Bλe−λθ = 0.

(8)Kuang (1993) derives explicit conditions for stability/instability of the n-thorder linear real scalar neutral diﬀerential diﬀerence equation with a single
6

delay. Since (8) is a special case of the n-th order equation, applying the resultof Kuang (1993, Theorem 1.2) implies that the real parts of the solutions ofequation (8) are positive for all θ if |B| > 1. Hence we have the following
result.
Theorem 1 If ν > ε, then the stationary point of (7) is unstable for all
θ > 0.
If v < ε (i.e., |B| < 1), (8) has at most finitely many eigenvalues withpositive real part. The roots of the characteristic equation are functions of thedelay. As the lengths of the delay change, the roots may change their signs
from positive to negative or vise versa so that the stability of the solution
may also change. Such phenomena are often referred to as stability switches.
We will next show that such stability switchings cannot take place in the
fixed delayed model.
For the following discussion we assume that v < ε. The case v = ε will betreated later as a critical case. It can be checked that λ = 0 is not a solutionof (8) because substituting λ = 0 yields A = 0 that contradicts A > 0. In thecase of v < ε, Kuang (1993, Theorem 1.4) shows that if the stability switchesat θ = θ̄, then (8) must have a pair of pure conjugate imaginary roots withθ = θ̄. Thus to find the critical value of θ̄, we assume that λ = iω, with ω > 0
is a root of (8) for θ = θ̄, θ̄ ≥ 0. Substituting λ = iω into (8), we have
A + Bω sin ωθ = 0,and
ω + Bω cos ωθ = 0.
Moving A and ω to the right hand side and adding the squares of the resultant
equations, we obtain
A2 + (1 − B 2 )ω 2 = 0.

Since A > 0 and 1 −B 2 > 0 as |B| < 1 is assumed, there is no ω that satisfies
the above equation. In other words, there are no roots of (8) crossing theimaginary axis when θ increases. Therefore, there are no stability switches
for any θ.
In case ε = ν in which |B| = 1, the characteristic equation becomesλ(1 − e−λθ ) + A = 0.

(9)

It is clear that λ = 0 is not a solution of (9) since A > 0. Thus we canassume that a root of (9) have non-negative real part, λ = u + iv with u ≥ 0
for some θ > 0. From (9), we have(u + A)2 + v2 = e−2uθ (u2 + v 2 ) ≤ (u2 + v 2 ),7

where the last inequality is due to e−2uθ ≤ 1 for u ≥ 0 and θ > 0. Hence
2uA + A2 ≤ 0,
where the direction of inequality contradicts the assumption that u ≥ 0 and
A > 0 in case ε = ν. Hence it is impossible that the characteristic equation
has roots with nonnegative real parts. Therefore, all roots of (9) must havenegative real parts for all θ > 0. Summarizing the above discussions gives
the following theorem.
Theorem 2 In case of v ≤ ε, the the stationary point of (7) is asymptoticallystable for all θ > 0.
Result 1 and Theorems 1 and 2 show that the stability condition of the
fixed delay model is the same as the one of the basic model. Thus it canbe seen that if the basic model is stable (resp. unstable), then the fixed
delay model is also stable (resp. unstable). In other words, introducing fixed
time lag does not change the stability condition of the basic model. It is,however, numerically confirmed that the fixed production lag has distinctive
eﬀects on the global dynamics as shown in Figure 2 in which we illustrate
four diﬀerent limit cycles for four diﬀerent values of θ. It can be observedthat the limit cycles change their shapes from parallelorgram-type cycle tovertically elongated parallelogram with rounded corners as the investment
lag increases from 0.125 to unity by doubling θ, that is, the width of thecycles becomes smaller and the height becomes larger. Since the peak andthe bottom of the output-cycle are reached at the point where the time
derivative is zero, the amplitude of the cycle is equal to the distance between
the two points at which the limit cycles cross the vertical axis. It is also
observed that the amplitude of the cycles increases and the rate of output
change becomes smaller as the production lag increases. These numericalresults are summarized in:
Result 2 The fixed time delay gets rid of discontinuous jumps and increasing
investment lag makes the length of cycles longer and its amplitude larger
in the case of ν > ε.

Insert Figure 2 about here.

8

2.3

Delay Model with Continuously Distributed Lags

Continuously distributed time delay is an alternative approach to deal with atime lag in investment. If the expected change of national income is denotedby ẏ e (t) at time t and is based on the entire history of the actual changes ofnational income from zero to t, the dynamic system can be written as thesystem of integro-diﬀerence equations,εẏ(t) − ϕ(ẏ e (t)) + (1 − α)y(t) = 0,eẏ (t) =

Z t
0

(10)
w(t − s, θ, m)ẏ(s)ds,

where the weighting function is
⎧
1 t−s
⎪
⎪
if m = 0,⎪ e− θ
⎨θw(t − s, θ, m) =³ ´⎪⎪
⎪ 1 m m+1 (t − s)m e− m(t−s)⎩θ
if m ≥ 1.
m! θHere m is a nonnegative integer and θ is a positive real parameter, whichis associated with the length of the delay. We call this dynamic system the
distributed delay model.
To examine local dynamics of the above system in the neighborhood of
the stationary point, we consider the linearized version,Z tw(t − s, θ, m)ẏ(s)ds + (1 − α)y(t) = 0.
εẏ(t) − ν
0Looking for the solution in the usual exponential formy(t) = y0 eλt and ẏ(t) = λy0 eλt ,we substitute y(t) and ẏ(t) into the linearized version to obtainZ tw(t − s, θ, m)e−λ(t−s) ds + (1 − α) = 0.
ελ − νλ0

As shown in Bischi et al, (2007), introducing the new variable z = t − ssimplifies the integral as
Z t
Z t−λ(t−s)w(t − s, θ, m)eds =w(z, θ, m)e−λz dz.
009By letting t → ∞ and assuming that Re(λ) + mθ > 0, we haveZ ∞1 − z −λze θ e dz = (1 + λθ)−1 if m = 0,
θ
0
and

Z ∞
0

That is,λθ
1 ³ m ´m+1 m − mz −λz
z e θ e dz = (1 + )−(m+1) if m > 1.
m! θ
m
Z ∞

−λz

w(z, θ, m)e

0

with

q=

⎧
⎨ 1
⎩

¶−(m+1)µλθds = 1 +q
if m = 0,m if m ≥ 1.

Then the characteristic equation becomes
µ
¶m+1
λθ
(ελ + (1 − α)) 1 +− νλ = 0.
q(11)

If there are no time delays, θ = 0, then the above equation is reduced to thesame characteristic equation as the one we have already derived above. We
will next show some simple cases in which analytical results can be obtained.
Since the case of m = 0 will be rigorously discussed in the next section,
we examine stability in cases with m ≥ 1. We expand the characteristic
equation (11) by using the binomial theorem to obtain,
a0 λm+2 + a1 λm+1 + ... + am+1 λ + am+2 = 0,(12)where the coeﬃcients ai are defined as
a0 = εθm+1 > 0,
ak =

½µ

m+1
k+1

¶

mε +

µ

m+1
k

¶

(1 − α)θ mk θm−k > 0 for k = 1, 2, ...m,am+1 = mm {mε + (m + 1)(1 − α)θ − mν} R 0,
am+2 = mm+1 (1 − α) > 0.

10¾

According to the Routh-Hurwitz criterion, the necessary and suﬃcient conditions that all roots of the characteristic equation (12) have negative realparts are the following:(1) the coeﬃcients are positive, ak > 0 for k = 1, 2, ...2m + 1,
(2) the principle minors of the Routh-Hurwitz determinant
are positive, ¯¯¯ a1 a0 0 0 ¯
¯
¯
¯
¯
¯
¯
¯ a1 a0 0 ¯¯ a3 a2 a1 a0 ¯¯
¯ a1 a0 ¯¯mmm¯¯ > 0, D3 = ¯ a3 a2 a0 ¯ > 0, D4 = ¯
D2 = ¯¯¯ a5 a4 a3 a2 ¯ > 0, ...
¯¯
a3 a2 ¯¯¯¯ a5 a4 a3 ¯¯ a7 a6 a5 a4 ¯Case 1. m = 1

Substituting m = 1 into (12) yields
a0 λ3 + a1 λ2 + a2 λ + a3 = 0,

(13)

where
a0
a1
a2
a3

=
=
=
=

εθ3 > 0,(2ε + (1 − α)θ)θ > 0,ε + (1 − α)2θ − ν R 0,
1 − α > 0.

It can be seen that the sign of a2 is not determined. In addition to a2 > 0, theRouth-Hurwitz criterion requires that the following second- and third-order
Routh-Hurwitz determinants are positive,¯¯¯¯¯ a1 a0 0 ¯
¯
¯
¯
¯
a a
D2m=1 = ¯¯ 1 0 ¯¯ > 0 and D3m=1 = ¯¯ a3 a2 a1 ¯¯ > 0.
a3 a2
¯ 0 0 a3 ¯Since D3m=1 = a3 D2m=1 and a3 = 1 − α > 0, we have¡¢¡¢sign D3m=1 = sign D2m=1 ,

where

©
ª
D2m=1 = θ 2(ε + (1 − α)θ)2 − (2ε + (1 − α)θ)ν .

(14)

Consider an a2 = 0 curve and a D2m=1 = 0 curve in the positive quadrantof the (θ, ν)-plane. First we assert that both curves are upward sloping andintersect only once for θ = 0. Moreover we can assert that a2 < 0 to the left
11

of the a2 = 0 curve and a2 > 0 to the right and that D2m=1 < 0 to the leftof the D2m=1 = 0 curve and D2m=1 > 0 to the right. Substituting a2 = 0 into
D2m=1 gives D2m=1 = −(1 − α)εθ2 < 0. This implies that the a2 = 0 curve islocated in the region where D2m=1 < 0. Thus, a2 > 0 and D3m=1 > 0 in theregion where D2m=1 > 0. Therefore the D2m=1 = 0 curve is the partition linethat divides the (θ, ν)-plane into two regions: one region below the line inwhich the stationary state is stable as the Routh-Hurwitz criterion is satisfied
and the other region above the line in which the stationary state is unstable.
Solving D2m=1 = 0 for ν yields the explicit expression of the partition line,ν=2(ε + (1 − α)θ)2.
2ε + (1 − α)θ

(15)

We now return to equation (13) to show the existence of a limit cycle
even with continuously distributed delay by applying the Hopf bifurcation
theorem. According to the theorem, we can establish the existence if the
cubic characteristic equation has a pair of pure imaginary roots and the realpart of these roots vary with a bifurcation parameter. We select ν as thebifurcation parameter and then calculate its value at the point for whichloss of stability just occurs. Substituting (15) into (13), we can obtain afactorized expression of the characteristic equation along the partition line,(2ε + (1 − α)θ + εθλ)(1 − α + (2εθ + (1 − α)θ2 )λ2 ) = 0,

which can be explicitly solved for λ. One of the characteristic roots is realand negative and the other two are pure imaginary:
2ε + (1 − α)θ
< 0,
λ1 = −εθs1−αλ2,3 = ±i= ±iω.
2εθ + (1 − α)θ2
In order to apply the Hopf bifurcation theorem, we need to check whether
the real part of the conjugate complex roots change its sign as the bifurcation
parameter passes through its critical value. Suppose that λ depends on ν,λ(ν), and then implicit-diﬀerentiation of (13) shows that
¡
¢
dλ
= λ.
3εθ3 λ2 + 2εθ + (1 − α)θ2 λ + ε + (1 − α)2θ − νdνThus( µ ¶ )¾½
−1dλd(Re λ)= sign Re
sign
dν
dν
λ=iω
λ=iω

©
ª
= sign (2εθ + (1 − α)θ2 )12

where we used the facts that the terms with λ are imaginary and the constantterms are real. Therefore we have¯d(Re λ) ¯¯> 0.
dν ¯λ=iωThis implies that the roots cross the imaginary axis at iω from left to rightas ν increases. Therefore the Hopf bifurcation theorem applies, and thus abirth of limit cycles is assured around the stationary point. The left part of
Figure 3 illustrates a limit cycle in a 3D space when the stationary state is
unstable, and the right side shown an orbit approaching the stationary state
when stable.2
Insert Figure 3 about here.

Case 2. m = 2
The characteristic equation is quartic in λ and its coeﬃcients are all
positive except a3 whose sign is not determined,
a3 = 4(2ε + 3(1 − α)θ) − 2ν).
The Routh-Hurwitz determinants can be defined in the same way as before,¯¯¯¯ a1 a0 0 0 ¯¯¯ a1 a0 0 ¯
¯
¯
¯
¯
¯
¯ a3 a2 a1 a0 ¯¯
¯ a1 a0 ¯
m=2
m=2
m=2
¯
¯
¯,
¯
¯
¯
, D3
=¯
= ¯ a3 a2 a1 ¯ and D4=¯
D2¯
a3 a2 ¯
0
a
a
a
4
3
2
¯ 0 a4 a3 ¯¯¯
¯ 0 0 0 a4 ¯
2z(t) and w(t) are defined as follows:z(t) =

Z t

1
− t−s
θ ẏ(s)ds
2 (t − s) e0 θ

andw(t) =
ν is selected as

and

Z t

1 − t−se θ ẏ(s)ds.
0 θ2(ε + (1 − α)θ)2+ 0.05 in the unstable case,2ε + (1 − α)θ
2(ε + (1 − α)θ)2− 0.05 in the stable case.
2ε + (1 − α)θ

13whereª©
D2m=2 = 2θ3 32ε2 + 3(1 − α)2 θ2 + ε(18(1 − α) + 4v) > 0,©
ª
D3m=2 = 8θ3 8(2ε + (1 − α)θ)3 − 2(28ε2 + 12(1 − α)εθ) + 3(1 − α)2 θ2 )ν − 8εν 2 ,D4m=2 = 8(1 − α)D3m=2 .

By the same way as given in the case of m = 1, we assert that the a3 = 0
curve and the D3m=2 = 0 curve are upward sloping, intersect once for θ = 0
and the former curve is located above the latter curve. Substituting a3 = 0
into D3m=2 gives
D3m=2 = −64(1 − α)2 θ4 (6ε + (1 − α)θ)2 < 0.
The direction of this inequality implies that the a3 = 0 curve is located in
the region in which D3m=2 < 0. Thus we have a3 > 0 and D4m=2 > 0 in theregion in which D3m=2 > 0. Hence the D3m=2 = 0 locus is the partition linethat divides the (θ, ν)-plane into two regions: one region below the line inwhich the stationary state is stable and the other region above the line inwhich the stationary state is unstable. From D2m=1 = 0 we obtain the explicit
expression of the partition line,8(2ε + (1 − α)θ)3 − 2(28ε2 + 12(1 − α)εθ) + 3(1 − α)2 θ2 )ν − 8εν 2 = 0. (16)
It can be confirmed that the D2m=1 = 0 curve is steeper than the D3m=2 = 0
curve, and both curves have the same vertical intercept for θ = 0. This meansthat the distributed delay model with m = 1 has a larger stable region thanthe model with m = 2. By the same procedure as in the case of m = 1 above,
we can show the birth of a limit cycle in the case of m = 2 as well.
After we repeat the above procedure for all values of m, it is then possible
to show that Dm+1 > 0 implies am+1 > 0 and Di > 0 for all i. Hence them
Dm+1= 0 locus is the partition line that divides the (θ, ν)-plane into two
regions, the stable region and the unstable region. The five partition lines
with m from 1 to 5 are depicted in Figure 4. It can be seen that all lines crossthe vertical axis for ν = ε and their slopes become smaller as m increases.
Notice that the dotted horizontal line is the partition line in the case of fixed
time delay. This implies that the stable region becomes smaller as the valueof m increases and converges to the region defined with the fixed time delay
when m tends to infinity.
Insert Figure 4 about here.
14

3

Coexistence of Multiple Cycles

In this section, we demonstrate a coexistence of multiple cycles when thestationary point is locally stable. This finding implies the robustness ofcyclic properties of Goodwin’s model regardless of whether the stationary
point is locally stable or unstable. In what follows, we select the slope of
the investment function evaluated at the stationary point as the bifurcation
parameter and investigate the possibility of the supercritical bifurcation in
Section 3.1. After proving the existence of an unstable limit cycle, we construct an invariant set in the state space and apply the Poincaré-Bendixson
theorem to find a stable limit cycle that encloses the unstable limit cycle in
Section 3.2. The coexistence of multiple cycles has been already shown for a
multiplier-accelerator model in Puu (1986), for Kaldor’s business cycle modelin Grasman and Wentzel (1994) and for a Metzlerian inventory cycle modelin Matsumoto (1996) using an approach that will be further explored.

3.1

Distributed Delay Model with m = 0

We consider the distributed delay model with m = 0 in which the weighting
function becomes exponentially declining and thus gives the most weight tothe most recent income change. The dynamic system of the Volterra-type
integro-diﬀerential equation isεẏ(t) − ϕ(ẏ e (t)) + (1 − α)y(t) = 0,eẏ (t) =Zt

1 − t−se θ ẏ(s)ds.
θ

(17)

0

As shown below, we have rewritten the continuously distributed delay equation as a system of ordinary diﬀerential equations. By doing so, we can use
all tools known from the stability theory of ordinary diﬀerential equations toanalyze the asymptotic behavior of the stationary point.
The time-diﬀerentiation of the second equation of (17) and introducing
z = ẏ e reduce it to1(18)ż(t) = (ẏ(t) − z(t)) .
θSolving the first equation for ẏ, replacing ẏ e with z, replacing ẏ in (18) with
the new expression of ẏ and then adding the new dynamic equation of z
will transform the integro-diﬀerential equation to the following 2D system of15

ordinary diﬀerential equations,1−α1y(t) + ϕ(z(t)),
ε
ε
µ
¶
1
1
1−α
ż(t) =−
y(t) + ϕ(z(t)) − z(t) .
θεεẏ(t) = −(19)

It can be seen that the stationary state of (19) is y = z = 0. By linearizing
the system at the stationary state, we get the Jacobian matrix,
⎞
⎛
1−α
ν
⎟
⎜ − ε
ε
⎟
⎜
J =⎜
⎟.
´ ⎠
⎝ 1−α
1 ³ν
−1
−
εθ
θ ε
The corresponding characteristic equation is quadratic in λ,3
λ2 +

1−α
ε + (1 − α)θ − ν
λ+
= 0.
εθ
εθ

Setting k = ν − [ε + (1 − α)θ] gives,
⎧
⎫
sµ ¶
2
⎨
1 k
4(1 − α) ⎬k
λ1,2 =±.
−
2 ⎩ εθεθεθ ⎭The product of the eigenvalues is positive,
λ1 λ2 =

1−α>0εθ

due to the assumptions imposed on parameters, 0 < α < 1 and (ε, θ) > 0.
These parametric restrictions ensure that the stationary point is not a saddle
point. It also follows that the sum of the eigenvalues is either positive or
negative according to whether k is negative or positive,
λ1 + λ2 =k
R 0 ⇔ k R 0.
εθ3
Notice that this characteristic equation is identical with the characteristic equationthat can be derived from Goodwin’s approximated model (i.e., equation (5f) of Goodwin(1951)). This means that both equations generates exactly the same dynamics in theneighborhood of θ = 0. See Szidarovszky and Matsumoto (2007) for the similarities anddissimilarities between the two dynamic equations.

16

To confirm the dependency of the stability on the parameters θ and ν, we
define the parameter region by Ω = {(θ, ν) |θ > 0 and ν > 0}, considering
the values of the other two parameters α and ε given. Since the stability ofthe stationary point depends of the sign of k so that a k = 0 curve becomesthe partition line that divides Ω into two regions, one stable and the otherunstable. Solving k = 0 gives the partition line
ν = [ε + (1 − α)θ],(20)

which is positive-sloping in Ω. Thus the stationary state is stable in theregion below this line in which k < 0 and unstable in the region above. It
also depends on the value of the discriminant of the characteristic equation
whether the local dynamics is oscillatory or monotonic. The curve along
which the discriminant is zero is determined bypν = [ε + (1 − α)θ] ± 2 (1 − α)εθ,(21)which distinguishes the parameter region for real roots from that for complex
roots.
These two curves divide the parameter region Ω as shown in Figure 5. For
combinations of θ and ν in a region marked as either [MS] (in the bottom
left and bottom right corners) or [MU], the stationary state is monotonic
stable and unstable while for those in [OS] or [OU], the stationary state isoscillatory stable and unstable.4 In Figure 5, we plot θ on the horizontal axis
and v on the vertical axis. For certain combinations of the parameters, θand ν, in either the light-gray region or the dark-gray region, the characteristic roots are complex and thus the system produces oscillations. Moreover,
the oscillations are explosive for the combinations in the light-gray region
and damped for those in the dark-gray region. On the other hand, for the
combinations of θ and ν in the white regions, the characteristic roots arereal and thus the system produces monotonic dynamics that is convergent ordivergent according to whether the combination is in the lower-white region
or the upper-white region. Figure 5 also implies that the distributed delay
model may produce qualitatively the same dynamics as ν increases regardlessof the value of the production lag. Thus, for a given θ, the linearized model
first generates monotonic stable dynamics, then oscillatory stable dynamics,
oscillatory unstable dynamics and finally monotonic unstable dynamics as ν
increases from zero.
4

The partition line crosses the vertical axis for ν = ε and is positive sloping while
the zero-discriminant locus is tangent to the vertical line for ν = ε and to the horizontal
εline for θ = 1−α. Thus the qualitatively similar divisions are obtained regardless of thespecified values of parameters, although we specify α = 0.6 and ε = 0.5 in Figure 5.

17Insert Figure 5 about here.

3.2

Hopf Cycle

By applying the Hopf bifurcation theorem, we investigate whether there isa limit cycle in the stable distributed delay model with m = 0. According
to the theorem, Hopf bifurcation occurs if the complex conjugate roots asfunctions of the bifurcation parameters cross the imaginary axis. Obviously,
the characteristic roots are complex conjugate with zero real part if k = 0.
As there are no other roots in the two-dimensional system, a limit cycle exists if the eigenvalues cross the imaginary axis with non-zero speed at thebifurcation point. Though there may exist several possibilities to parametrize the distributed delay model, it seems interesting to choose the slope of
the investment function evaluated at the stationary point as the bifurcation
parameter.
Since we consider the case of local stability (i.e., k < 0 or ν < ε+(1−α)θ)in this section, it can be easily seen that there exists a value ν 0 for whichν 0 − [ε + (1 − α)θ] = 0,
implying that the complex conjugate roots cross the imaginary axis. For
ν > ν 0 , respectively ν < ν 0 , the real part becomes positive, respectively
negative. Hence, v0 is indeed a bifurcation value of the distributed delay
model. Due to the Hopf bifurcation theorem, this establishes the existenceof closed orbits in a neighborhood of the stationary point (0, 0) at ν = ν 0 .5
The Hopf theorem, however, has no indication about the nature of thelimit cycle. There are two possibilities, one is that orbits spiral outward
from the stationary point toward a stable limit cycle, called the supercriticalbifurcation, and the other is that all orbits starting inside the cycle spiral in
toward the stationary point and becomes explosive outside the cycle, called
the subcritical bifurcation. To make the distinction between the sub- andsuper-critical Hopf bifurcation, we calculate the stability index.
The distributed delay model can be written asµ ¶ µ 1
µ ¶
¶
g (z)
y
ẏ
+
=J
z
ż
g 2 (z)5

See Lorenz (1993) for the Hopf bifurcation theorem, the stability index and the following coordinate transformations.

18

where J is the Jacobian matrix defined above, and gi (z) for i = 1, 2 arenonlinear terms that can be derived as1
g 1 (z) =(ϕ(z) − ν 0 z) ,θ
1
(ϕ(z) − ν 0 z) .
εθIn order to transform the Jacobian matrix of the distributed delay model
into the normal form, we introduce the coordinate transformation
µ
µ ¶
¶
µ ¶
0
1
u
y
with D ==D,
d21 d22vz
g 2 (z) =wherer
ν0 − εε
ν0 − εd21 =and d22 =
.
ν 0θ
ε
ν 0θ
Since matrix D transforms the coordinate system (y, z) into a new coordinate
system (u, v), the distributed delay model becomes
r
⎛
⎞
1−α
¶
µ ¶ µ
µ ¶
0
−
⎜
⎟ u
ag(u,
v)
u̇
εθ
⎟
r
=⎜
⎝ 1−α
⎠ v + bg(u, v)
v̇
0
εθ

where

g(u, v) = ϕ(d21 u + d22 v) − ν 0 (d21 u + d22 v).
1
1
and b = .
a = pεε(ν 0 − ε)It is well-known that the stability of the emerging cycle depends on up
to the third-order derivatives of the nonlinear function g(u, v). The stability
index is1I =[a (guuu + guvv ) + b (guuv + gvvv )]
16 r
¡
¢
¢¤
1
εθ £¡ 2+a − b2 guv (guu + gvv ) + ab (guu )2 − (gvv )2
16 1 − αwhere the partial derivatives are

guu = ϕ00 (d21 )2 , gvv = ϕ00 (d22 )2 , guv = ϕ00 d21 d22 ,guuu = ϕ000 (d21 )3 , gvvv = ϕ000 (d22 )3 , guvv = ϕ000 d21 (d22 )2 , guuv = ϕ000 (d21 )2 d22 .
19Then arranging terms gives
I=

ϕ000 (0)ν 0 (ν 0 − ε)
> 0.
16εθ3

This inequality indicates that the emerging cycle, enclosing the stable stationary point, is repelling (i.e., unstable). A subcritical Hopf bifurcation
occurs for ν < ν 0 in the distributed delay model with m = 0. We cansummarize our results as follows.
Theorem 3 Given the stability of the stationary point, there is an unstable
limit cycle that encloses a stable equilibrium in the distributed delay modelwith m = 0.

3.3

Stable and Unstable Limit Cycles

Assuming that ν < ν 0 , we investigate a number of a limit cycle in thissubsection. As shown in Theorem 3, an unstable limit cycle exists for v < v0
due to the Hopf theorem. We will also show that a stable limit cycle exists
and it encloses the unstable Hopf cycle. We will finally construct an invariant
set in such a way that once an orbit enters the set, it cannot escape fromit at any future time and then will apply the Poincaré-Bendixon theorem to
examine whether a stable cycle can arise in the set. A typical application of
the theorem is in the case in which there is a single unstable equilibrium insome invariant set. Remove the neighborhood of the stationary point in which
all orbits move away from it. Then the theorem ensures that the remaining
set must contain a limit cycle.6 Since the stationary point is now assumed
to be stable in the distributed delay model, we need another approach forsearching for such an invariant set.
In Figure 6 below, we may find a point A located so high on the verticalaxis that an orbit starting there comes back to cross again the vertical axis
at a lower point E after crossing not only the horizontal axis at points B
and D but also the vertical axis at point C. Consider a region bounded by
ABCDEA and an open region bounded by the unstable limit cycle that isshown to exist in Theorem 3. Then the invariant set to be considered can beobtained by deleting the latter region from the former region. The result isshown in Figure 6 in which the shaded region represents this invariant set.
The boundary of the inner white region corresponds to the unstable limit
cycle that surrounds the stable stationary point. The shaded region thus has
no stationary point, and any orbit starting inside the region stays within this6

See Chaper 2.2 of Lorenze (1993) for example.

20

region. Formally, the Poincaré-Bendixson theorem guarantees the existenceof one stable limit cycle in the shaded region as illustrated as the bold cyclein Figure 6. Combining the discussion above with Theorem 3 yields thefollowing:
Theorem 4 Given the stability of the stationary point, a stable limit cycle
coexists with an unstable limit cycle that encloses a stable equilibrium point
for the distributed delay model with m = 0.

Insert Figure 6 HereIn Figure 7, we present a bifurcation diagram in which the amplitude of
the cycle is on the vertical axis and the bifurcation parameter ν on the horizontal axis. ν u is the critical value at which the distributed delay model loses
its stability, and a limit cycle is born. For ν > ν u , the model is destabilized
for small perturbations so that any orbit moves away from the stationary
point. Nonlinearity of the model prevents it from diverging globally but
leads to a unique stable limit cycle. This is essentially the same cycle as theone that Goodwin (1951) demonstrates by applying the Lienard method. On
the other hand, for ν < ν u , it is locally stabilized but generates an unstable
cycle as well as a stable cycle for ν in the interval [ν s , ν u ]. It can be seen thatas ν decreases from ν u , the amplitude of the inner unstable cycle increases
and that of the outer stable cycle decreases. vs is the other critical value forwhich the two cycles coincide. For a lower ν < ν s , limit cycles no longer exit
since the invariant set vanishes.
The coexistence of multiple cycles and a stable stationary point reminds
us to "corridor stability," the notion which was introduced by Leijonhufvud(1973). It implies that a dynamic system is stable for small perturbations
but unstable for large perturbations. The interior of the white elliptic regionin Figure 7 is the corridor in which the stationary point is stable. Thus,
if perturbations around the stationary point are small enough not to takethe trajectory out of the corridor, the restoring eﬀect of the model forces the
trajectory to return to the stationary point. To the contrary, if perturbations
are large enough to take orbits outside the corridor, the lasting eﬀect of themodel leads to persistent cyclic oscillations. The foregoing numerical analysis
indicates that the nonlinearity of the investment function and the time delay
can be the sources of corridor stability in the distributed delay model withm = 0.
Insert Figure 7 Here21

4Concluding Remarks

We reconsidered Goodwin’s 1951 nonlinear accelerator model of business cycle and demonstrated two new features, the eﬀects caused by investment lagon the characteristic of Goodwin’s cycle when the stationary state was locallyunstable and the coexistence of multiple business cycles when locally stable.
Concerning the first feature, it is numerically confirmed that increasing investment lag makes the length of the business cycle longer and its amplitudelarger. Moreover, it is analytically confirmed that the fixed investment lag
has the stronger destabilizing eﬀect in comparison to the continuously distributed investment lag.
Concerning the second feature, we showed, by combining the result obtained from the Hopf bifurcation theorem with the one due to PoincareBendixon theorem, that two limit cycles can coexist with the stable stationary state: one cycle is unstable and surrounds the stationary state, and theother is stable and encloses the unstable limit cycle. This finding indicatesthe corridor stability of Goodwin’s model in which a damping force dominatesand makes trajectories approach the stationary state for small disturbances
but an anti-damping force dominates and makes trajectories converge to theouter stable limit cycle for larger disturbances. The results imply global
stability of Goodwin’s model regardless of the local dynamic properties.

22References[1] Bischi, G., C. Chiarella, M. Kopel and F. Szidarovszky, Nonlinear
Oligopolies: Stability and Bifurcations, mimeo.
[2] Goodwin, R., "The Nonlinear Accelerator and the Persistence of Business Cycles," Econometrica, vol.19, 1-17, 1951.
[3] Grasman, J., and J. Wentzel, "Coexistence of a limit cycle and an equilibrium in Kaldor’s Business Cycle Model and its Consequences," Journal of Economic Behavior and Organization, vol.24, 369-377,1994.
[4] Kuang, Y., Delay Diﬀerential Equations with Applications in Population
Dynamics, Boston, Academic Press, 1993.
[5] Leijonhufvud, A., "Eﬀective Demand Failures," Swedish Journal of Economics, 27-48, 1973.
[6] Lorenz, H-W., Nonlinear Dynamical Economics and Chaotic Motion,
Springer-Verlag, 1993.
[7] Lorenz, H-W., "Goodwin’s Nonlinear Accelerator and Chaotic Motion,"
Journal of Economics, vol. 47, 413-418, 1987.
[8] Lorenz, H-W, and H. E. Nusse, "Chaotic Attractors, Chaotic Saddles,
and Fractal Basin Boundaries: Goodwin’s nonlinear Accelerator Model
Reconsidered," Chaos, Solitions and Fractals, 13, 957-965, 2002.
[9] Matsumoto, A., "A Non-linear Macro Model of Endogenous Inventory
Oscillations," Research in Economics, vol.51, 101-129, 1997.
[10] Puu, T., Attractors, Bifurcations, and Chaos, Second Edition, SpringerVerlag, 2003.
[11] Sasakura, K., "The Business Cycle Model with a Unique Stable Limit
Cycle," Journal of Economic Dynamics and Control, vol.20, 1763-1773,
1996.
[12] Szidarovszky, F., and A. Matsumoto, "Asymptotic Stability of the linear
diﬀerential equaion and the delay diﬀerential equation with continuously
distributed lag," mimeo, 2007.
[13] Yoshida, H., and T., Asada, "Dynamic Analysis of Policy Lag in a
Keynesian-Goodwin Model: Stability, Instability, Cycles and Chaos,"Journal of Economic Behavior and Organization, vol.62, 441-469, 2007.
23

Caption of Figures

Figure 1. Limit cycle without time lag.

Figure 2. Limit cycles with time lags.

Figure 3. Limt cycle and stable trajectory in 3D system.

Figure4. Partion lines.
Figure 5. Stable and unstable regions.

Figure 6. Coexistence of multiple limit cycles.

Figure 7. Bifurcation diaram.

24

yHtL

yHt L

A

B

yH
I0

yL

I0

5

10

15

20

y° HtLt
C

Figure 1:

25

D

yHtL

q= 1

q=0.5

q=0.25q=0.125

Figure 2:26

y° HtL

zHtL
0

zHt L
0

2

1

-1

-2

2

1

wHtL 0wHtL 0

-2

-1

0

2
yHt L

-2
-1

4
6

Figure 3:

27

0

1
yHt L

2

3

n

m=1
m=2

m=3
m=4
m=5

e

n=e

q
0.2

0.4

0.6

Figure 4:

28

0.8

1

v

@MU D

e

@OUD

@OSD

@MS D
eê1 -aq

Figure 5:

29

y
A
y =0

E

z= 0

B

D
C

Figure 6:

30

z

y

6

4

2

νu

ν0

-2

Figure 7:31ν