IERCU
Institute of Economic Research, Chuo University
50th Anniversary Special Issues

Discussion Paper No.208

Dynamics in Delay IS-LM Model with Tax
Collections
Akio Matsumoto

Ferenc Szidarovszky

Chuo University

University of Pécs

September 2013I

INSTITUTE OF ECONOMIC RESEARCH
Chuo University
Tokyo, JapanDynamics in Delay IS-LM Model with Tax
Collections
Akio Matsumoto
Chuo Universityy

Ferenc Szidarovszky
University of Pécsz

Abstract
An IS-LM model is developed for the dynamics of income, interest rate
and money stock with delay in tax revenue. The main aim is to show thatthe delay matters in macro dynamics. Two di¤erent delays, …xed timedelays and continuously distributed time delays, are considered explicitly
and described by delay-di¤erential equations and integro-di¤erential equations, respectively. Conditions for the local stability of the two models arederived and compared. The destabilizing e¤ects caused by the delay are
numerically examined. Emergence of wide spectrum of dynamics ranging from simple cyclic oscillations to complex dynamics involving chaos is
described through Hopf bifurctions.

Keywords: Fixed time delay, Continuously distributed time delay,
Hopf bifurcation, Stability switch, IS-LM model

JEL number: E12, E32, E37, E62

The authors highly appreciate the …nancial suppports from the MEXT-Supported Program for the Stratgic Research Foundation at Private Universities 2013-2017, the Japan Society for the Promotion of Science (Grant-in-Aid for Scienti…c Research (C) 24530202 and
25380238) and Chuo University (Grant for Special Resarch). The usual disclaimers apply.
y Professor, Department of Economics, Chuo Univeresity, 742-1, Higahi-Nakano, Hachioji,
Tokyo, 192-0393, Japan; akiom@tamacc.chuo-u.ac.jp
z Professor, Department of Applied Mathematics, University of Pécs, Ifúság, u. 6, H-7624,
Pécs, Hungary; szidarka@gmail.com

1

1

Introduction

Since the pioneering work of Kalecki (1935) and the seminal work of Goodwin(1951), it has been recognized that economic dynamic systems usually incorporate delays in their actions and delay is one of the essentials for macroeconomic‡uctuation. Nevertheless, little attention has been given to studies on delay ineconomic variables over the past few decades. After a long "gestation" period oftime, however, the number of studies on delay gradually increases and various
attempts have been done on the impact of delays on macro dynamics.1 Amongothers, we draw attention to the papers of De Cesare and Sportelli (2005) andFanti and Manfredi (2007), (DS and FM henceforth) and re-examine the stability/instability e¤ects caused by a delay in the tax collection from a slightlydi¤erent point of view. Both papers introduce time delay into a simple IS-LM
model with a pure money …nancing de…cit, which are used to show the existenceof cyclic ‡uctuations of the macro variables (Schinasi (1981, 1982) and Sasakura(1994)). Noticing the established fact that there are delays in collecting tax,
DS concerns "economic situations where a …nite time delay cannot be ignored"
and investigates how the …xed time delay (…xed delay henceforth) in tax collection a¤ects the …scal policy outcomes. Two results are shown: the emergenceof limit cycle through Hopf bifurcation when the length of the delay becomes
longer and the co-existence of multiple stable and unstable limit cycles when
the steady (equilibrium) point is locally stable. On the other hand, FM replaces
the …xed time delay with the continuously distributed time delay (distributed
delay henceforth), emphasizing the evidence that there is "a wide variation incollection lag" and demonstrates the possibility that in the same IS-LM framework, complex dynamics involving chaos is born through an ála period-doublingbifurcation with respect to the length of the delay.
The purpose of this paper is to further explore what dynamics the delay
IS-LM model can generate.2 Specifying some of parameters, we analyticallyand numerically examine delay e¤ects caused by changing the length of delay
and tax e¤ects caused by changing the value of the tax coe¢ cient to show that"delay really matters." The followings are demonstrated:
1 In the literature of time delay, there are two di¤erent modeling ways in continuous-timescale: …xed time delay and continuously distributed time delay The former is applicablewhen an institutionally or socially determined …xed period of time delay is presented for theagents involved. The latter is appropriate for economic situations in which di¤erent lengthsof delays are distributed over heterogeneous agents or the length of the delay is random
depending on unforeseeable circumstance. The choice of the type of delay has situationdependency and results in the use of di¤erent analytical tools. In the cases of …xed timedelay, dynamics is described by a delay di¤erential equation whose characteristic equation isa mixed polynomial-exponential equation with in…nitely many eigenvalues. In the cases ofcontinuously distributed time delay, Volterra-type integro-di¤erential equations are used tomodel the dynamics and with appropriate weighting function their characteristic equations
have …nitely many eigenvalues. A comprehensive summary of the theory of …xed delays is given
in Bellman and Cooke (1956), and the mathematical methodology of dealing with distributed
delays is presented in Cushing (1977).
2 In Matsumoto and Szidarovszky (2013), we unify these two models and consider dynamics
in an IS-LM model augmented with both …xed and distributed delays.

2

(I) A …xed delay IS-LM model can generate not only cyclic ‡uctuations which
are qualitatively di¤erent from those in the non-delay model but also complex dynamics;(II) Comparing the distributed model with the …xed model reveals that the
distributed model is more stable than the …xed model in the sense thatthe former has a larger stable region than the latter;(III) An increase of the shape parameter of the weighting function of the distributed model has a destabilizing e¤ect by decreasing the stability region;(IV) If the shape parameter converges to in…nity, then the stability region
converges to that of the …xed delay model.
This paper develops as follows. In Section 2, the …xed delay in tax collection
is introduced and the delay and tax e¤ects are examined. In Section 3, the…xed delay is replaced with the distributed delay. After conducting stability
analysis, we numerically show that stability is sensitive to the shape of thedelay kernel and analytically show that the distributed model converges to the…xed model when the delay kernel tends to the Direc delta function. Section 4o¤ers conclusions and further research directions.

2

Fixed Delay IS-LM Model

As a benchmark for macro stability analysis, we start with a non-delay IS-LM
model with a pure money …nancing de…cit:
8
Y_ (t) = [I(Y (t); R(t))
>
>
>
>
<
_
(MN ) :R(t)= [L(Y (t); R(t))
>
>
>
>
: _
M (t) = g T (t);s (Y (t)T (t)) + g

T (t)] ;M (t)] ;where the three state variables, Y; R and M; respectively represent income,interest rate and real money supply, the parameters, a; ; g and s are positive
adjustment coe¢ cients in the markets of income and money, constant government expenditure and the constant marginal propensity to save and I( ) andL( ) denote the investment and liquidity preference functions. Tax revenue isdenoted by T and is collected as a lump sum with a constant rate, 0 < < 1;
T (t) = Y (t):Following DS, we specify the investment and money demand functions as
I(Y; R) = A

Yand L(Y; R) = Y +
R
R

3

(1)

_
with positive parameters A; and : The conditions, Y_ (t) = R(t)= M_ (t) = 0
determine the unique steady state (Y ; R ; M ) such that
Y =

g

; R =

A
s(1

)

and M = Y +

R

:

(2)

Observing the facts that any economic system has delays in collecting tax
and the length of such delay can be thought as the same for all agents in thesame tax system, DS makes a simpli…ed assumption that tax revenue at time
t consists of two complementary components; one based on the current income
and the other on past income,T (t) = (1") Y (t) + " Y (t)(3)

where is a time delay and " < 1 is a positive constant. Substituting (3) into
(MN ) transforms the non-delay model into a …xed-delay model,
8
Y (t)
>
>
Y_ (t) =A[s + (1 s)(1 ") ]Y (t) + g (1 s)" Y (t
) ;
>
>
R(t)
>
>
>
>
<
(MF ) :_R(t)=Y (t) +
M (t) ;
>
>
R(t)
>
>
>
>
>
>
: _
M (t) = g (1 ") Y (t) " Y (t):
DS has already shown two main results, the existence of limit cycle throughHopf bifurcations and the co-existence of multiple cycles. In this section, we
will show that the …xed delay model (MF ) gives rise a much wider variety ofdynamic behavior including chaos.
The stationary state of the modi…ed model (MF ) is the same as in the nondelay model, (Y ; R ; M ). The linearized version of model (MF ) is
1 0
Y_ (t)a11@ R_ (t) A = @ a21a31
M_ (t)
0

a12
a22
0

10
1 0
0
Y (t)
a23 A @ R (t) A + @M (t)0b1 Y (t0b2 Y (t

)
)

1
A

(4)

where the elements of the matrices are evaluated at the stationary point andgiven bya11 =(1 + (1s)"); a12 =

a22 =

s2 (1
)2
; a23 =
A2

s2 (1

)2
A

; a31 =

g; a21 =

;
(5)

(1

")

and
b1 = (1s)" and b2 = " :4(6)

We also look for the exponential solutions, Y (t) = e t u; R (t) = e t v andM (t) = e t w; and substitute these into equation (4) to obtain the characteristic
equation,P ( ) + Q( )e
=0
(7)
where
P( ) =

3(a11 + a22 ) 2 + (a11 a22

a12 a21 )a12 a23 a31

(8)

and
Q( ) = b1 2

a22 + a12 a23 b2 :(9)

If we introduce notation
p2 =(a11 + a22 ); p1 = a11 a22

a12 a21 ; p0 =a12 a23 a31

andq 2 = b 1 ; q1 =a22 ; q0 = a12 a23 b2 ;then equation (7) can be rewritten as3+ p2 2 + p1 + p0 + q 2 2 + q 1 + q 0 e

= 0:

(10)

Suppose that = i!; ! > 0; is a root of this equation. Substituting it into
the equation and then separating the real and imaginary parts, we havep2 ! 2 + p0 + (q0! 2 q2 ) cos + !q1 sin ! = 0(11)! 3 + !p1 + !q1 cos !

(q0! 2 q2 ) sin ! = 0:

Thus
(q0

! 2 q2 )2 + (!q1 )2 = (p2 ! 2p0 )2 + (! 3

!p1 )2 :Hence ! must be a positive root of the following equation! 6 + (p222p1q22 )! 4 + (p21

q12 + 2q0 q2 )! 2 + (p20

2p0 p2q02 ) = 0:

Let z = ! 2 : Then the sextic equation can be reduced to a cubic equation,(z) = z 3 + az 2 + bz + c = 0
with

a = p22

q22

2p1 ;

b = p21 + 2q0 q2
c = p20
The discriminant of

2p0 p2

q12 ;

q02 :

(z) isD = a2 b2 + 18abc

4b3
5

4a3 c

27c2 :(12)

(z) has one real root and a conjugate pair of imaginary roots if D < 0 andthree distinct real roots if D > 0.3 Notice that
p0

s2 (1

q0 =

)2
A

(1

2")

(13)

which is negative if " > 1=2 and positive if " < 1=2. The rest of this section isdivided into two subparts: Dynamics in the case of " > 1=2 is …rst considered
and then in the case of " < 1=2:

2.1

Delay tax is dominant: " > 1=2

We take " = 3=5 in this section.4 Returning to equation (3), we see that morethan half of the tax revenue is delayed when " > 1=2: We specify the parametervalues and formulate this speci…cation as an assumption since we repeatedly use
this speci…cation of the parameters in numerical studies to be employed below.
Assumption== A = 1;

= 4=5;= 3; s = 1=5 and g = 10:

Under Assumption, clearly b > 0 and c < 0 while it depends on a value of
whether a is positive or negative. In Figure 1(A), the upward-sloping black
curve is the a = 0 curve indicating that a < 0 to its left and a > 0 to its
right. The sign of the discriminant is also -dependent. D > 0 in the light blue
regions, D < 0 in the ocher region and D = 0 on their boundaries. Along thehorizontal line at g = g, three critical values of are determined so that thediscriminant is zero. We denote them as i for i = 1; 2; 3 in ascending order.
The horizontal line also crosses the a = 0 curve at a : Having the signs of a andD with respect to and then applying Descartes’ rule, we can determine thenumber of positive roots along the horizontal line5 :(1) one positive root (and a pair of imaginary roots) if 0 < < 1 ;(2) three distinctive positive roots if 1 < < 2 ;(3) one positive root (and a pair of imaginary roots) if 2 < < a ;(4) one positive root (and a pair of imaginary roots) if a < < 3 ;(5) one positive root (and two negative roots) if 3 < < 1:
Figurep 1(B) illustrates how the number of the positive roots and the valueof ! = z depend on : A di¤erent color curve corresponds to a di¤erent root.
3 Introducing a new variable, x = z + a=3; reduces(z) to the simpli…ed form,

3x + px + q
in which the discriminant is

q 2
p 2
+
2
3
and p and q are de…ned accordingly. We can arrive at the same results to be obtained even
if this form of the discriminant is used.
4 The results to be obtained in this section do not depend on the selection of the particularvalue of " as far as " > 1=2:
5 In the same way, we can determine this number for any other value of g:6

Solving (z) = 0 yields at most three positive solutions, zi for i = 1; 2; 3:6 Thepp
red curve is associated to ! 1 = z1 , the blue curve to ! 2 = z2 and the greenpcurve to ! 3 = z3 . It can be seen that three di¤erent-color curves are illustrated
for 2 ( 1 ; 2 ) and two di¤erent color curves are coincident for = 1 and= 2 when equation (12) has equal positive roots. Only one curve is illustrated
for any other values of . We can summarize the results on the number of thepositive roots and the associated purely imaginary roots of equation (12):Proposition 1 Given Assumption and " = 3=5, there are three purely imaginary solutions, j = i! j for j = 1; 2; 3 with ! 3 > ! 2 > ! 1 > 0 if 1 < < 2 ;two imaginary solutions with ! 3 = ! 2 > ! 1 or ! 3 > ! 2 = ! 1 if = 1 or
= 2 and one imaginary solution = i! 1 if < 1 or = i! 3 if > 2 .

(A) Region division(B) Real roots

Figure 1. Positive roots of the characteristic equation (10)From (11), we obtain two solutions! 4 q2! 2 (q0 + p1 q2 p2 q1 ) + p1 q0(q0 ! 2 q2 )2 + (!q1 )2! 4 (q1

p2 q2 ) + ! 2 (p0 q2 + p2 q0 p1 q1 )(q0 ! 2 q2 )2 + (!q1 )2

sin ! = !

p0 q 1;(14)

cos ! =

p0 q 0:

There is a unique value of ! , 0 < !
2 satisfying both equations. Taking
an inverse of the second equation7 and taking account of the fact that the cosine
function has the same value for every 2 into account, the critical value of delay
is determined as= '( ; n)
6 We do not give an explicit form of each z

i as it is long and complicated.
7 Solving the …rst equation yields the same value of :7

where '( ; n) is de…ned'( ; n) =

1
cos 1
!

! 4 (q1

p2 q2 ) + ! 2 (p0 q2 + p2 q0 p1 q1 )(q0 ! 2 q2 )2 + (!q1 )2

p0 q 0+ 2n(15)for n = 1; 2; ::. Substituting a positive !; which is functional for , into thisequation determines the partition curve that partitions the ( ; ) plane by separating the stable region of it from the unstable region. In Figure 2, the stationary
state is locally stable in the yellow region and unstable in the white region. Thus
the boundary between these regions is the boundary curve which consists of thesegments of = '( ; n) for n = 0; 1. The color of the partition curve is the sameas the color of the positive root curve as depicted in Figure 1(B). The distorted
L-shaped continuous curve with three di¤erent color segments is the partition
curve with n = 0: Its lowest part colored in green crosses the horizontal axis
twice at = A and = B and becomes negative between these intersections.
So this captures the fact that if there is no delay (that is, = 0), then the steadystate is locally unstable for between these two values.8 The two slender anddeformed U -shaped curves having blue and green segments are the = '( ; n)curves for n = 1 and n = 2; respectively. The = '( ; n) curve shifts upward
when the value of n increases.9 Given the length of delay, the steady state loses
stability for the tax rate on the boundary of the yellow region. Apparently such
critical values of the tax rate are -dependent and di¤erent from A , for which
the steady state loses stability in the non-delay model.
Having determined the partition curves, we now examine the delay e¤ ect
and the tax e¤ ect on stability. The former is caused by a change in the lengthof delay and is examined by changing along a vertical line at some . On theother hand, the latter is caused by a change in the tax rate with a positive delay
8 More explanations are given later.
9 The red segments with n

1 have much larger

8

and are located outside Figure 4.

and obtained by changing

along a horizontal line at some .

Figure 2. Partition curves with " > 1=2

2.1.1

Tax E¤ect

We study the tax e¤ect. Again as a bench mark, we …rst simulate the model(MF ) with = 0 or the model (MN ) twice. Under Assumption, A ' 0:29
and B ' 0:8: The initial values of Y and R are the same in both simulations,
Y (0) = Y and R(0) = R but the initial values of M are di¤erent, M (0) =M + 1 in the …rst simulation and M (0) = M + 5 in the second simulation,
to verify the initial point dependency. The resultant bifurcation diagrams arepresented in Figures 3(A) and 3(B), in each of which the downward sloping black
curve depicts the equilibrium value of output, Y e = g= . In the simulations, thebifurcation parameter is increased from 0:2 to 1 with an increment of 1=1000;
the iterations are repeated 5000 times and the local maximum and minimumof y(t) for the last 100 iterations are plotted against each value of . In Figure
3(A), the bifurcation diagram of the …rst simulation is depicted. It is observedthat the stationary state loses stability when arrives at A and bifurcates toa periodic cycle having one maximum and one minimum for > A : It is also
observed that the oscillation disappears at = B and stability is regainedfor > B : In Figure 3(B), the bifurcation diagram in the second simulation
is illustrated. It is seen that stability is lost at = A as in Figure 3(A) but
regained at some value larger than B : Further simulations with di¤erent initialpoints have been conducted and then lead to the fact that stability is regained
not necessarily at = B but at some larger value although stability is always
lost at = A . These numerical results are summarized as follows:9Proposition 2 Given Assumption, the non-delay model (MN ) generates periodic oscillations when the steady state is destabilized and has initial point
dependency when it regains stability.

(A) m(0) = M + 1(B) m(0) = M + 5

Figure 3. Bifurcation diagrams with di¤erent initial values
We now look into the case of positive delay to examine whether changing thevalue of the positive delay essentially a¤ects the results obtained without delay.
To this end, we perform simulations with four di¤erent values, 1 = 0:5; 2 =4, 3 = 9; 4 = 14 and two di¤erent constant initial functions '1 (t) = M + 1
and '2 (t) = M + 5 for t
0: We choose
as the bifurcation parameterand increase its value from zero to unity along the horizontal dotted lines at= i for i = 1; 2; 3; 4 in Figure 2. Given i ; let s denote the critical valueof the tax rate for which the horizontal line at i crosses the partition curvefrom left for the …rst time. The resultant bifurcation diagrams are given in
Figures 4(A), 4(B), 4(C) and 4(D) where stability is lost for = s and theisoelastic black curve describes the equilibrium level of output. In each …gures,the diagram with '1 (t) is colored in red, the one with '2 (t) in blue and thered diagram is superimposed on the blue diagram. The following results are
obtained. There is no big di¤erence between Figure 3 and Figure 4(A), thesteady state loses stability at = s ; bifurcates to a limit cycle having one
maximum and one minimum and regains stability at a value depending on theselection of the initial function. Figure 4(B) shows the birth of oscillations
having several maxima and minima, indicating that the delay matters. Notice
two issues: one is that the blue parts can be seen for larger values of being
close to unity, implying that the selection of the initial function a¤ects dynamics
there; the other is that its e¤ect is economically harmless because some partsof a time trajectory become negative (i.e., economically meaningless) when it is10

e¤ective. Figure 4(C) yields di¤erent dynamics. The steady state loses stability
at = s1 ; regains it at = s2 and then loses stability again at = s3 : As aresult of a larger delay, oscillations for > s3 exhibit more ups and downs thanthose in Figure 4(B). In the …nal example, much more complicated oscillations
are observed and dynamics becomes economically infeasible when becomes large
as in Figure 4(D).
From these numerical simulations we obtain the following two results. First,
di¤erent initial functions generate only limited e¤ects on economically feasible
dynamics having positive trajectories. In other words, the tax e¤ect has only
minor initial function dependency. Second, the shape di¤erences among thebifurcation diagrams reveal that the tax e¤ect with a large value of generates
substantial e¤ects on output dynamics. Concerning the stability switch, we …ndthat there are two di¤erent ways. Stability loss and regain in the delay model
occur for very small values of in the same way as in the non-delay model.
The stability regain takes place only once for a relatively large value of : We
notice in Figure 2 that the lower end of the U -shaped curve with n = 1 breaks
into the stable yellow region and makes that region unstable. If a value of isappropriately selected such that a horizontal line at this crosses the lower partof the U -shaped curve, then increasing values of along this line induce stability
switches three times for relatively smaller values of and …nally destabilize thesteady state. The results of numerical simulations are summarized as follows:
Proposition 3 The tax e¤ ects substantially depend on the delay at least in three
ways: (i) the delay model does not have a strong initial function dependency
regardless of the length of delay; (ii) unstable trajectories exhibit more erratic
ups and downs as the length of delay becomes longer; (3) stability switch canoccur for some value of delay.

(A)

1 = 0:5

(B)

11

2 =4

(C)

3 =9

(D)

4 =4

Figure 4. Bifurcation diagrams with respect to
2.1.2

Delay E¤ect

To verify the delay e¤ects caused by changing the length of delay on stability, we
select delay as the bifurcation parameter and examine the sign of the derivativeof Re[ ( )] at the points where ( ) is purely imaginary. Di¤erentiating equation(12) with respect to yields
d
d

1

=

(

3

3 2 + 2p2 + p1
2q2 + q1+
+ p2 2 + p1 + p0 )(q2 2 + q1 + q0 )where we study, only for convenience, the inverse of d =d . Substituting =i!; rationalizing the denominator and noticing that the terms with and 3
are imaginary while the constant and the terms with 2 and 4 are real yield
the following form of the real part of the derivative of :
"
#
1
0
d
(! 2 )
Re
= 2 2
d
! (!
p1 )2 + (p2 ! 2 p0 )2=i!

where the denominator is de…nitely positive. It is seen that the sign of thereal part is positive for ! i with 0 (! 2 ) > 0 and negative for ! i with 0 (! 2 ) <0: Indeed, it can be veri…ed that if equation (! 2 ) = 0 has three positiveroots, ! 1 ; ! 2 ; ! 3 ; then the derivative of the root has the following directions ofinequality:
"
#
"
#
"
#
d
d
d
> 0; Re< 0; and Re
>0
(16)
Re
d =i!1
d =i!2d =i!312

and otherwise(! 2 ) = 0 has only one positive root with
Re

d
d

> 0:

(17)

=i!

The positive sign implies that all roots cross the imaginary axis from left toright as increases, that is, stability is lost. On the other hand, the negative
sign means crossing the axis in the reverse direction, that is, stability is gained.
We now perform two simulations with di¤erent values of . In the …rst
simulation, we increase from zero to 20 along the vertical dashed line at == '( ; n) (n = 0; 1; 2)
A (' 0:287) in Figure 2. The vertical line crosses thecurves six times. Three intersections with the lower, middle and upper blue
curves are at 1 ' 1:59, 4 ' 13:7 and 6 ' 19:1: Similarly two intersections
with the middle and upper green curves are at 2 ' 8:23 and 5 ' 16:4: Finally
3 ' 10:3 is the intersection with the red curve. The subscripts of i are labeledin the ascending order. Numerical result is illustrated in Figure 5(A). The
starting point of = A and = 0 is on the lowest green segment so that thereal part of the eigenvalue is zero and its derivative with respect to is positive.
The steady state loses stability when becomes positive and bifurcates to aperiodic cycle. As the length of delay becomes larger, the cycle expands, shrinks
and then merges with the stationary state at = 1 to gain stability. It loses
stability again when increasing crosses the green partition curve from below at= 2 : Repeating the same process, the steady state regains stability at = 3 .
The second switch from instability to stability occurs. Passing through a small
stability interval of ; the steady state is destabilized when increasing arrives
at the downward-sloping red partition curve at = 4 . For > 4 ; the steadystate remains unstable and a periodic cycle with seven extremal values emerges.
In the second simulation, is increased to 1 = 0:6: As seen in Figure 2,
the dotted vertical line at = 1 does intersect with the partition curve andthe steady state is locally unstable even at = 0. Comparing Figure 5(A)with Figure 5(B) implies that ‡uctuations becomes complicated as increases.
Period-doubling and period-halving bifurcations can be observed. Althoughcomplex ‡uctuation emerges for a larger value of ; its time trajectory in thegray region of Figure 5(B) takes negative value. Notice that the di¤erence
between these two simulations is only the value of : We summarize the results:
Proposition 4 Delay a¤ ects output dynamics and generates a shape di¤ erent
bifurcation cascade depending on the value of the tax rate.: (1) when a smaller
value of the tax rate is given, increasing …nally destabilizes the steady state
after repeating several stability switches; (2) when a larger value of the tax rate
is given, no stability switch occurs and trajectories exhibit more ups and downs
as increases

13

(A)

=

A (' 0:287)(B)=

1 (= 0:6)

Figure 5. Bifurcation diagrams with respect to

2.2

Current tax is dominant: " < 1=2

In this section we take " = 2=5 and consider the case where non-delay tax is
dominant over delay tax.10 For " = 2=5, c > 0 always, however, the signs ofa; b and D are ambiguous. In Figure 6(A), the loci of a = 0 and b = 0 are
illustrated in red and blue and both are upward-sloping. It can be checked thata > 0 and b < 0 to the right of the corresponding locus and the inequality isreversed to the left. The D = 0 locus is depicted in black and takes a L-shaped
curve. D > 0 in the light-blue region to the right of the locus and D < 0 in theyellow region to the left. Along the g = g horizontal line; three critical valuesof , a ' 0:63; b ' 0:49 and D ' 0:27; are obtained so that a = 0; b = 0
and D = 0, respectively. Applying Descartes’ rule of signs, we can determine
the number of positive roots along this horizontal line according to the signs ofthese parameters and the discriminant:(1) no positive roots if < D ;(2) two positive roots (and one negative root) if(3) two positive roots (and one negative root) if(4) two positive roots (and one negative root) if

D <

b;

b <

a;

a <

<1

where the two positive roots are equal for = D : Dependency of the two rootsp! i = zi (i = 1; 2) on the value of the tax rate is described by the blue andgreen downward sloping curves shown in Figure 6(B). We summarize the results
as follows:
1 0 As in the previous case, this particular value is selected only for convenience. We could
have the qualitatively same results with di¤erent values of " as far as " < 1=2:

14

Proposition 5 Given Assumption, there are two purely imaginary solutions,> D (or D > 0) and no imaginaryi = i! i for i = 1; 2 with ! 2 > ! 1 > 0 if
solutions if < D (or D < 0).

(A) Region division(B) Real roots with g = 10

Figure 6. Positive roots of characteristic equation (12) with " = 2=5
Similarly to the case of " = 3=5; we can obtain the partition curves in the( ; ) plane by substituting ! i (i = 1; 2) into equation (15). Five distorted U shaped curves for n = 0; 1; 2; 3; 4 are depicted in Figure 7 where the steady stateis asymptotically stable in the yellow region. Since the curve shifts upward asn increases, the lowest curve is with n = 0 and the highest curve with n = 4.
Each curve consists of blue and green segments associated with two positiveroots with the same color in Figure 6(B). For < D ; the steady state is locallyasymptotically stable irrespective of the length of delay. Since delay does nota¤ect dynamics, such delay is called harmless. For > D , a signi…cant change
in the stability region takes place. The boundary of the stability region is madeup of in…nitely many arcs of the partition curves. The most remarkable e¤ect of
the change is that as the delay is increased, many stability switches may occur
for a …xed value of ; especially for closer to D . This is in sharp contrast withthe case of " > 1=2 where, after several stability switches, further increased invalue of could not induce re-stabilization of the steady state when it crosses
the red segment of the partition curve shown Figure 2.

15

Figure 7. Partition curves with " = 2=5
In order to study the global stability properties of the steady state, we …rst
…x the parameter values, = 20 and = 0:55 and then perform two numericalsimulations.11 In the …rst simulation, we examine the tax e¤ect with positive
delay by increasing from zero to unity along the horizontal dotted line at= . The bifurcation diagram with respect to is illustrated in Figures 8(A)where stability is regained twice indicating that the tax rate a¤ects the stability
as well as the instability. It is observed that unstable trajectories converge
to periodic cycles whose periodicity increases as the tax rate becomes larger.
Period-doubling and -halving process are also observed for 2 ( 3 ; 4 ): In thesecond simulation, we detect the delay e¤ect by increasing from zero to 40
along the vertical dotted line at = . The bifurcation diagrams with respect
to is illustrated in Figure 8(B) where the delay e¤ect seems to be similar tothe delay tax e¤ect. Indeed, stability switches take place twice and unstabletrajectories exhibit simple dynamics for smaller values of ; period-doublingand period-halving occur for medium values (i.e., 2 ( 2 ; 3 )) and complex
dynamics emerges for larger values (i.e., > 3 ).
Proposition 6 Given Assumption and ! = 2=5; delay model has multiple stability switches and can generate simple dynamics to complex dynamics as thebifurcation parameter ( or ) becomes larger.
1 1 Since it is veri…ed that the initial function dependency is weak, we take the same initial
function, '(t) = M + 1; in both simulations.

16

(A) Tax e¤ect ( = 20)(B) Delay e¤ect ( = 0:55)

Figure 8. Bifurcation diagrams

3

Distributed Delay IS-LM Model

As mentioned in the Introduction, FM concerns the situation where there isa wide variation in collecting tax due to heterogenous economic agents. FM
replaces the …xed delay in tax collection with distributed delay and show occurrence of chaotic oscillations. It adopts a very special form of tax collection,
T (t) = Y e (t) where Y e (t) is the expected income at time t and is a weighted
average of past incomes. We reconsider the FM version of the IS-LM model
augmented with the continuously distributed time delay. Our version is di¤erent in the following two points. First, the general form of tax collection given inequation (3), T (t) = (1 ") Y (t) + " Y e (t); is assumed, which includes the special form when " = 1: Second, we use a bell-shaped delay kernel instead of the
exponentially declining delay kernel. Fixed delay model (MF ) is now modi…ed
as a distributed delay model,
8
Y (t)
>
>
Y_ (t) =A+ g (s + (1 s)(1 ") Y (t) (1 s)" Y e (t)
>
>
R(t)
>
>
>
>
<
(MD ) :_R(t)=Y (t) +
M (t)
>
>
R(t)
>
>
>
>
>
>
: _
M (t) = g (1 ") Y (t) " Y e (t)withY e (t) =

Z t

W (t

0

17

; ; m)Y ( )dwhere the weighting function or the delaying kernel is de…ned by
8 1 1
(t s)
>
if m = 0;
>
< e
W (t
; ; m) =>>m: 1 m m+1 (t)m e (t s) if m 1:
m!

(18)

Here m is a nonnegative integer and is a positive real parameter, which corresponds to the average length of delay. The expected output formed at time
t is the weighted average of the actual demand in the past. According to definition (18), the shape of the weighting function is determined by the value ofthe shape parameter m: For m = 0; weights are exponentially declining withthe most weight given to the most current data (i.e., a fading memory). FM
con…nes attention to this case and show that the distributed delay IS-LM model
can display stability switch from stability to instability through a Hopf bifurcation and gives rise to chaotic motions when the average length of the delays
is relatively large. For m
1; zero weight is given to the most current data,rising to maximum as = t
and declining exponentially thereafter. Thus
the weights take a bell-shaped form which becomes taller and thinner as mincreases. We draw attention mainly to the case with m 1:
To examine the local stability of this dynamical system in a neighborhood of
the steady state, we linearize the model to relate it to the relevant parameters.
We start with the case of m = 1 in which the expected output is
Z t
2
1
1
e
(t
)e (t s) Y ( )dY (t) =
0and its rate of change is1
Y_ e (t) = (X(t)

Y e (t)) :(19)

Here the new variable X(t) is de…ned as
Z t
1 1 (t s)
X(t) =e
Y ( )d0

and its time derivative is1_
X(t)= (Y (t)X(t)) :

(20)

Adding di¤erential equations (19) and (20) to the integro-di¤erential system(MD ) yields a 5D dynamical system of ordinary di¤erential equations. Its Jacobi
matrix is
10
a11 a12 0b1
0
B a21 a22 a23
0
0 C
B
C
B
0
b2
0 C
J m=1 = B a31 0
C
1 A
1
@ 0
0
0
1
1
0
0
0
18

where aij and bi are the same as those given in (5) and (6), respectively. Expanding det (J m=1
I) = 0 yields the characteristic equation of the Jacobi
matrix J m=1 ;
1+1

P( ) 1 ++ Q( ) = 0;1

where P ( ) and Q( ) are already given in (8) and (9). In the same way, thecharacteristic equations with m = 2 and m = 3 are obtained as
2+1

P( ) 1 +2

3+1+ Q( ) = 0 and P ( ) 1 +3+ Q( ) = 0:

By induction, the characteristic equation has the general form
m+1

P( ) 1 ++ Q( ) = 0;mwhich is an (m + 4)-th order polynomial equation
a0 m+4 + a1 m+3 + ::: + am+3 + am+4 = 0
where the coe¢ cients ai are de…ned accordingly.

3.1

Stability Analysis

Since it is di¢ cult to obtain a general solution of the (m + 4)-th degree characteristic equation, we con…ne attention to special cases with m = 1; 2; 3 andm ! 1: We will analytically consider local stability of the steady state and thennumerically examine global dynamics. Before proceeding to stability analysis,
we construct the Routh-Hurwitz determinant:10
a1 a0 0 0
0
B a3 a2 a1 a0
0 C
CB
B a5 a4 a3 a2
0 C
C:
Jm = det B
B a7 a6 a5 a4
0 C
C
B
@
0 A
0 0 0 0am+4k
Let Jmdenote the k-th order leading principal minor of Jm :In the case of m = 1; the characteristic function is a 5th degree polynomial
and its coe¢ cients are numerically con…rmed to be all positive for 0"119and 0
a1

=

a2

=

a3

=

a4

=

a5

=

1 under Assumption,
50 + (3 + (19 20") + 3 2 )> 0;
25
125 + 10 (3 + (19 20 ) + 3 2 ) + 2 (1)2 (40 + 3(5 4") 2 )
> 0;
125 2
2
50(1 ") 2 (1)2 + 5 (3 + 19 + 3 2 ) + 2 (1)(40 + 3(5 4") 2 )
> 0;
2
125
(1
)2 (8 + 20(1 ") + 3 2 )
> 0;
25 2
2(1
)2
> 0;
5 2

where a0 = 1. Furthermore it is also numerically veri…ed that " = 3=5 and" = 2=50
1
a1 a0 0aa1
0J12 = det
> 0 and J13 = det @ a3 a2 a1 A > 0:
a3 a2
a5 a4 a3

Notice that J14 > 0 implies J15 > 0 since a5 > 0: So, according to the RouthHurwitz stability criterion, the following determinant needs to be positive forpreserving local stability,01
a1 a0 0 0
B a3 a2 a1 a0 C
C
J14 = det B@ a5 a4 a3 a2 A :0 0 a5 a4

Under the set of parameters given in Assumption, the J14 = 0 locus is thepartition curve. In Figures 9(A) and 9(B), the partition curves with " = 3=5
and " = 2=5 correspond to the right most curves arrowed by m = 1, respectively.
The steady state is locally stable to the left of the curve and unstable to the right(i.e., the white region). Two other curves are the loci of J25 = 0 and J36 = 0,
that is, the partition curves with m = 2 and m = 3: The yellow regions are thestability regions under the …xed time delay.12 From these numerical examples,
we obtain the following three results. The …rst is that the distributed time delay
models are more stable than the …xed delay model in the sense that the stability
region of the former is larger than the one of the latter. The partition curve
shifts leftward as m becomes larger, shrinking the stable region. So the second
result is that increasing the value of m has a destabilizing e¤ect. The third isthat di¤erent values of the weight " between the current and past taxes do not

1 2 In Figure 9(A), only the locus ofnience.

= '( ; n) for n = 0 is depicted for graphical conve-

20

generate qualitatively di¤erent dynamics as is seen in Figures 9(A) and 9(B).

(A) " = 3=5(B) " = 2=5Figure 9. Region division and partition curves
The stable region with distributed delay becomes smaller as the value of mincreases and converges to the region de…ned by the …xed delay when m tends
to in…nity. This result is natural if we notice the properties of the delay kernel.
The kernel with m 1 is bell-shaped and becomes more peaked around t
as m increases. Furthermore it tends to the Dirac delta function if m ! 1. In
consequence, for su¢ ciently large values of m; the delay kernel may be regarded
as very close to the Dirac delta function and the dynamic behavior under the
distributed delay is very similar to that under the …xed delay. We may explain
this result mathematically by noticing that the characteristic equation with
continuous delay can be written as(m+1)

P ( ) + Q( ) 1 +

= 0:

m

If m tends to in…nity, then the left hand side converges to
P ( ) + Q( )e=0which is identical with the characteristic equation (7) of the delay di¤erential
equation with a single …xed delay.13 In short, under distributed delay, although
we comprehensively use all delayed or past output data, the stability domain issensitive to the shape of the delay kernel.
Proposition 7 Increasing the value of the shape parameter m has a destabilizing e¤ ect in the sense that it decreases the stability region of the distributed
1 3 Notice that the same notation
models.

is used to describe the delay in the …xed and distributed

21

delay model, and m ! 1 makes the distributed delay model converge to the…xed delay model.
Taking " = 3=5; we conduct the …nal simulations of the delay and tax e¤ects
in the distributed model. In the …rst example, we increases the value of ; thelength of delay, from 0 to 10 along the dotted vertical line at = 1 that crosses
the partition curve with m = 1 for = 1 in Figure 9(A). Figure 10(A) depicts
a bifurcation diagram with respect to ; in which cyclic oscillations around theunstable steady state emerge for < 1 ; the steady state regains stability at= 1 and is locally asymptotically stable for > 1 : In the second example,
we increase the value of ; to verify the tax e¤ect, along the dotted horizontal
line at = 1 that crosses the partition curve with m = 1 for = 1 : Asseen in Figure 10(B), the steady state loses stability on the partition curve andbifurcates to cyclic oscillations for > 1 :(A) Delay e¤ect ( = 0:55)(B) Tax e¤ect ( = 20)

Figure 10. Bifurcation diagrams with " = 3=5

4

Concluding Remark

Dynamic IS-LM model was considered with …xed and continuously distributed
delays in past income. The non-delay model generates periodic oscillations withrespect to the tax rate when the steady state loses stability and has initial point
dependency when stability is regained. If the delay is …xed and its length isthe bifurcation parameter, then the dynamic properties depend on whether thecurrent tax is dominant or not. If it is, then the model has multiple stability, switches and can generate simple dynamics to complex dynamics as thebifurcation parameter increases. If not, then multiple stability switches occur
and …nally further increase leads to instability. If the delay is continuously distributed, then the stability region decreases as the shape parameter m of the22

weighting function increases, and as m ! 1; it converges to the stability regionof the …xed delay model.
The existence of …xed and continuously distributed delays or multiple …xed orcontinuously distributed delays can be examined in a similar manner, however
both the analytical study and the numerical simulations could become more
complicated.

23

References[1] Bellman, R. and Cooke, K. L. Di¤ erential-di¤ erence Equations. Academic
Press, New York, 1956.
[2] Cushing, J., Integro-di¤ erential Equations and Delay Models in Population
Dynamics, Springer-Verlag, Berlin/Heidelberg/New York, 1977.
[3] De Cesare, L. and Sportelli, M., A Dynamic IS-LM Model with Delay
Taxation Revenues, Chaos, Solitions and Fractals, 25, 233-244, 2005.
[4] Fanti, L. and Manfredi, P., Chaotic Business Cycles and Fiscal Policy: An
IS-LM Model with Distributed Tax Collection Lags, Chaos, Solitions
and Fractals, 32, 736-744, 2007.
[5] Gandolfo, G., Economic Dynamics,
Berlin/Heidelberg/ New York, 2010.

4th Edition,

Springer-Verlag,[6] Goodwin, R. The Nonlinear Accelerator and the Persistence of Business
Cycles, Econometrica 19, 1-17, 1951.
[7] Kalecki, M., A Macrodynamic Theory of Business Cycles, Econometrica,
3, 327-344, 1935.
[8] Matsumoto, A. and Szidarovszky, F., Dynamics in an IS-LM model augmented with Fixed and Distributed Time Delays, mimeo, 2013.
[9] Sasakura, K., On the Dynamic Behavior of Schinasi’s Business Cycle
Model, Journal of Macroeconomics, 16, 423-444, 1994.
[10] Schinasi, G., A Nonlinear Dynamic Model of Short Run Fluctuations, Review of Economic Studies, 48, 649-659, 1981.
[11] Schinasi, G., Fluctuations in a Dynamic, Intermediate-run IS-LM model:
Applications of the Poincare-Bendixon Theorem, Journal of Economic
Theory, 28, 369-375, 1982.

24