\begin{document}

\title{f Hopf bifurcation in a dynamic IS-LM model with time delay}

\begin{tabular}{cccccccc}
\scriptsize{$^{a}$Department of Economic Informatics and Statistics,}\\
\scriptsize{Faculty of Economics, West University of Timi\c soara,}\\
\scriptsize{str. Pestalozzi, nr. 16A, 300115, Timi\c soara, Romania,}\\
\scriptsize{E-mails:mihaela.neamtu@fse.uvt.ro, cchilarescu@rectorat.uvt.ro}\\
\scriptsize{$^{b}$ Department of Applied
Mathematics, Faculty of Mathematics,}\\
\scriptsize{West University of Timi\c soara, Bd. V. Parvan, nr. 4, 300223, Timi\c soara, Romania,}\\
 \scriptsize{e-mail: opris@math.uvt.ro}\\

\end{tabular}

\begin{quote} \small{{\bf Abstract.} In this paper we investigate the impact of
delayed tax revenues on the fiscal policy out-comes.
Choosing the delay as a bifurcation parameter we study the
direction and the stability of the bifurcating periodic solutions.
With respect to the delay we show when the system is stable. Some
numerical examples are finally given for justifying the
theoretical results.}
\end{quote}

\noindent{\small{{\it Keywords:} delay differential equation,
stability, Hopf bifurcation, IS-LM model.

\noindent{\it 2000 AMS Mathematics Subject Classification: 34K20,
34C25, 37G15, 91B64} .}}

\section*{\normalsize\bf 1. Introduction}

\hspace{0.6cm}

The differential equations with time delay play an important role
for economy, engineering, biology and social sciences, because a
lot of phenomena are described with their help. In this paper, we
consider a model from economy, of the IS-LM type with time delay
and we study how the delay affects the macroeconomic stability.
The Hopf bifurcation and normal form theories are tools for
establishing the existence and the stability of the periodic
solutions. Similar ideas can be found in [5], [9], [10].

In [2] Cesare and Sportelli taking into account the papers [6],
[7], [8], study a dynamic IS-LM model of the following type:
$$\begin{array}{l}
\dot Y(t)\!=\!\alpha[I(Y(t),r(t))\!+\!g\!-\!S(Y(t)\!-\!T(Y(t), Y(t\!-\!\tau))))\!-\!T(Y(t), Y(t\!-\!\tau)))]\\
\dot r(t)=\beta[L(Y(t),r(t))-M(t)]\\
\dot M(t)=g-T(Y(t), Y(t\!-\!\tau))),\end{array}$$ with Y as
income, I as investment, g as government expenditure (constant), S
as savings, T as tax revenues, r as rate of interest, L as
liquidity, M as real money supply and $\alpha, \beta$ as positive
constants. The time delay $\tau$ appears in function T:
$$T(Y(t), Y(t\!-\!\tau)))=d(1-\varepsilon)Y(t)+d\varepsilon Y(t-\tau),\eqno(1)$$
where $d\in (0,1)$ is a common average tax rate and
$\varepsilon\in (0, 1)$ is the income tax share.

Based on the papers [2], [1], [11], we consider the following
IS-LM model:
$$\begin{array}{l}
\dot Y(t)\!\!=\!\!\alpha[I(Y(t),r(t))\!+\!g\!-\!S(Y(t)\!\!-\!\!T(Y(t), Y(t\!-\!\tau)))\!-\!T(Y(t), Y(t\!-\!\tau))]\\
\dot r(t)=\beta[L(Y(t), r(t))-M(t)]\\
\dot K(t)=I(Y(t-\tau), r(t))-\delta K(t)\\
\dot M(t)=g-T(Y(t), Y(t\!-\!\tau)),\end{array}\eqno(2)$$ with
$\deltalta>0$ and with the initial conditions: $$Y(t)=\varphi(t),
t\in[-\tau, 0], r(0)=r_1, M(0)=M_1, K(0)=K_1.$$

In the following analysis we will consider the function T given by
(1), the investment, the saving  and the liquidity of the form:
$$\begin{array}{l}I(Y(t), r(t))=aY(t)^{\alpha_1}r(t)^{-\alpha_2},\quad a>0, \alpha_1>0,
\alpha_2>0,\\
S(Y(t)-T(Y(t), Y(t-\tau)))=s(Y(t)-T(Y(t), Y(t-\tau))),\quad s\in(0,1),\\
L(r(t))=mY(t)+\displaystyle\frac{\gamma_0 }{r(t)-r_2}, \quad m>0, \gamma_0>0,
r_2>0.\end{array}\eqno(3)$$

The paper is organized as follows. In section 2 we investigate the
local stability of the equilibrium point associated to system (2).
Choosing the delay as a bifurcation parameter some sufficient
conditions for the existence of Hopf bifurcation are found. In
section 3 there is the main aim of the paper, namely the
direction, the stability and the period of a limit cycle solution.
Section 4 gives some numerical simulations which show the
existence and the nature of the periodic solutions. Finally, some
conclusions are given.

\section*{\normalsize\bf 2. The qualitative analysis of system (2).}

\hspace{0.6cm}

Using functions (1) and (3), system (2) becomes:
$$\begin{array}{l}
\dot Y(t)\!=\!\alpha[((s\!-\!1)(1\!-\!\varepsilon)d\!-\!s)Y(t)\!+\!d\varepsilon(s\!-\!1)Y(t\!-\!\tau)\!+\!aY(t)^{\alpha_1}r(t)^{-\alpha_2}\!+\!g]\\
\dot r(t)=\beta[mY(t)+\displaystyle\frac{\gamma_0}{r(t)-r_2}-M(t)]\\
\dot K(t)=aY(t-\tau)^{\alpha_1}r(t)^{-\alpha_2}-\delta K(t)\\
\dot M(t)=g-d(1-\varepsilon)Y(t)-d\varepsilon
Y(t-\tau),\end{array}\eqno(4)$$ $\alphapha, \beta>0$, $\alphapha_1,
\alphapha_2>0$, $m>0$, $\gamma_0>0$, $r_2>0$, $\deltalta>0$, $s\in
(0,1)$, $d\in (0,1)$, $\varepsilon\in (0,1)$.

System (4) is a system of equations with time delay. The
qualitative analysis is done using the methods from [3].

The equilibrium point of system (4) has the coordinates $Y_0, r_0,
K_0, M_0$, where:
$$Y_0\!=\!\displaystyle\frac{g}{d}, r_0\!=\![\displaystyle\frac{s(1\!-\!d)}{a}Y_0^{1\!-\!\alpha_1}]^{-\frac{1}{\alpha_2}},
K_0\!=\!\displaystyle\frac{s(1\!-\!d)Y_0}{\deltalta},
M_0\!=mY_0\!+\!\displaystyle\frac{\gamma_0}{r_0\!-\!r_2}.\eqno(5)$$

Using the translation
$$Y=x_1+Y_0, r=x_2+r_0, K=x_3+K_0, M=x_4+M_0$$
in (4) and considering the Taylor expansion of the right members
from (4) until the third order, we have:
$$\dot x(t)=Ax(t)+Bx(t-\tau)+F(x(t), x(t-\tau))\eqno(6)$$
where
$$A\!\!=\!\!\left(\!\!\begin{array}{cccc}

\alpha a_1 & a\alphapha\rho_{01} & 0 & 0\\

\beta m & \beta\gamma_0\gamma_1 & 0 & -\beta\\

 0 & a\rho_{01} & -\deltalta &0\\
a_4 & 0 & 0 & 0
\end{array}\!\!\right)\!\!, B\!\!=\!\!\left(\!\!\begin{array}{cccc}

\alpha b_1 & 0 & 0 & 0\\

0 & 0 & 0 & 0\\

a\rho_{10} & 0 & 0 & 0\\

 b_4 & 0 & 0 & 0\end{array}\!\!\right)\eqno(7)$$
where $x(t)=(x_1(t), x_2(t), x_3(t), x_4(t))^T$, and
$$\begin{array}{l}
a_1=(s-1)d(1-\varepsilon)-s+a\rho_{10}, b_1=(s-1)d\varepsilon,
a_4=-d(1-\varepsilon), b_4=-d\varepsilon\\
\rho_{01}=-\alpha_2Y_0^{\alpha_1}r_0^{-(\alpha_2+1)},
\rho_{10}=\alpha_1Y_0^{\alpha_1-1}r_0^{-\alpha_2}\end{array}$$ and

$$F(x(t), x(t-\tau))=(F^1(x(t), x(t-\tau)), F^2(x(t), x(t-\tau)),
F^3(x(t), x(t-\tau)),\eqno(8)$$ $F^4(x(t), x(t-\tau)))^T$,

$ F^1(x(t), x(t-\tau))= \displaystyle\frac{\alphapha a}{2}\rho_{20}x_1^2(t)+\alphapha
a\rho_{11}x_1(t)x_2(t)+\displaystyle\frac{\alphapha a}{2}\rho_{02}x_2^2(t)+$

           $\displaystyle\frac{\alphapha a}{6}(\rho_{30}x_1^3(t)+3\rho_{21}x_1(t)^2x_2(t)+3\rho_{12}x_1(t)x_2(t)^2+\rho_{03}x_2(t)^3)$

$$F^2(x(t), x(t-\tau))= \displaystyle\frac{\gamma_0}{2}\gamma_2x_2(t)^2
+\displaystyle\frac{\gamma_0}{6}\gamma_3x_2(t)^3$$

 $F^3(x(t), x(t-\tau))=
\displaystyle\frac{a}{2}\rho_{20}x_1(t-\tau)^2+a\rho_{11}x_1(t-\tau)x_2(t)+\displaystyle\frac{a}{2}\rho_{02}x_2(t)^2+$

$+\displaystyle\frac{a}{6}(\rho_{30}x_1^3(t-\tau)+3\rho_{21}x_1(t-\tau)^2x_2(t)+3\rho_{12}x_1(t-\tau)x_2(t)^2+\rho_{03}x_2(t)^3)$

$F^4(x(t), x(t-\tau))=0$

and

$$ \rho_{20}=\alpha_1(\alpha_1-1)Y_0^{\alpha_1-2}r_0^{-\alpha_2}, \quad
\rho_{11}=-\alpha_1\alpha_2Y_0^{\alpha_1-1}r_0^{-(\alpha_2+1)}$$
$$ \rho_{02}=\alpha_2(\alpha_2+1)Y_0^{\alpha_1}r_0^{-(\alpha_2+2)}, \quad
\rho_{30}=\alpha_1(\alpha_1-1)(\alpha_1-2)Y_0^{\alpha_1-3}r_0^{-\alpha_2}$$
$$ \rho_{21}=-\alpha_1\alpha_2(\alpha_1-1)Y_0^{\alpha_1-2}r_0^{-(\alpha_2+1)},
\quad \rho_{12}=\alpha_1\alpha_2(\alpha_2+1)Y_0^{\alpha_1-1}r_0^{-(\alpha_2+2)}$$
$$ \rho_{03}=-\alpha_2(\alpha_2+1)(\alpha_2+2)Y_0^{\alpha_1}r_0^{-(\alpha_2+3)},$$
$$\gamma_1=-\displaystyle\frac{1}{(r_0-r_2)^2},
\gamma_2=\displaystyle\frac{2}{(r_0-r_2)^3}, \gamma_3=-\displaystyle\frac{6}{(r_0-r_2)^4}.$$

    The characteristic equation of linear part from (6) is:
$$det(\lambdabda I-A-Be^{-\lambdabda\tau})=(\lambdabda+\deltalta)\Delta(\lambdabda,
\tau)=0\eqno(9)$$ where
$$\Delta(\lambdabda,
\tau)=P(\lambdabda)+e^{-\lambda\tau}Q(\lambda)\eqno(10)$$

$$P(\lambda)=\lambdabda^3+p_2\lambdabda^2+p_1\lambda+p_0, \quad Q(\lambda)=q_2\lambdabda^2+q_1\lambda+q_0$$ and
$$\begin{array}{l}
p_2=-(\alpha a_1+\beta\gamma_0\gamma_1),\quad
p_1=\alpha\beta(\gamma_0a_1\gamma_1+ma\rho_{01}), \quad p_0=\alpha a\beta a_4\rho_{01}\\
q_2=-\alpha b_1, \quad q_1=\alpha\beta\gamma_0 b_1\gamma_1, \quad
q_0=\alpha\beta b_4\rho_{01}.\end{array}$$

To investigate the local stability of the equilibrium point, we
begin by considering, as usual, the case without delay ($\tau=0$
). In this case the characteristic polynomial is:
$$(\lambda+\deltalta)(P(\lambdabda)+Q(\lambda))=0$$
hence, according to the Hurwitz criterion, the equilibrium point
is stable if and only if:
$$p_2+q_2>0, \quad (p_1+q_1)(p_2+q_2)>p_0+q_0.$$

When  $\tau>0$, standard results on stability of systems of delay
differential equations postulate that a equilibrium point is
asymptotically stable if an only if all roots of equation (10)
have a negative real part. It is well known that equation (10) is
a transcendental equation which has an infinite number of complex
roots and the some possible roots with positive real part are
finite in number.

We want to obtain the values $\tau_0$  such that the equilibrium
point (5) changes from local asymptotic stability to instability
or vice versa. We need the imaginary solutions of equation
$\Delta(\lambda, \tau)=0$. Let $\lambda=\partialm i\omega$ be these solutions and
without loss of generality we assume $\omega>0$. We suppose that
$P(i\omega)+Q(i\omega)\neq0$, for all $\omega\in{\rm{I \! R}}$. The previous
conditions are equivalent to $(p_1+q_1)(p_2+q_2)\neq p_0+q_0$.

A necessary condition to have $\omega$ as a solution of
$\Delta(i\omega,\tau)=0$ is that $\omega$ must be a root of the
following equation:
$$f(\omega)=\omega^6+a_F\omega^4+b_F\omega^2+c_F=0\eqno(11)$$
where $a_F=p_2^2-q_2^2-2p_1$, $b_F=2q_0q_2-2p_0p_2-q_1^2+p_1^2$,
$c_F=p_0^2-q_0^2$.

Let $k=-\displaystyle\frac{a_F}{3}$ and
$f_D=\displaystyle\frac{1}{4}[f(k)]^4+\displaystyle\frac{1}{9}[f^{'}(k)]^3$ .

Using the results from [2],  it results:

{\bf Proposition 1.} {\it 1. Let $\varepsilon>\displaystyle\frac{1}{2}$. Then
the following cases can be discerned:

(i) If $a_F\geq0$ or $b_F\leq 0$, then equation (11) has only one
real positive root;

(ii) If $a_F<0$, $b_F>0$  and $f_D\leq0$, then equation (11) has
only one real positive root.

 2. If $\varepsilon<\displaystyle\frac{1}{2}$,
$a_F<0$, $b_F<0$ and $f_D\leq0$ then equation (11) has two real
positive roots which are distinct if $f_D\neq0$.

3. If $\varepsilon=\displaystyle\frac{1}{2}$ and:

(i) $a_F<0$, $b_F=0$ then equation (11) has only one real positive
root;

(ii)$a_F<0$, $b_F>0,$ then equation (11) has two real positive
roots;

(iii)$b_F<0,$ then equation (11) has only one real positive root.}

Also, we have:

{\bf Theorem 1.} {\it  If we suppose that the equilibrium point
$(Y_0, r_0, K_0, M_0)$ is locally asymptotically stable without
time delay, then in conditions of Proposition 1 there exists only
one stability switch.}

{\bf Theorem 2.} {\it If $\tau_0$ is a stability switch and
$f_D\neq0$, then a Hopf bifurcation occurs at $\tau_0$, where
$$\tau_0=\displaystyle\frac{1}{\omega_0}arctg\left (\displaystyle\frac{\omega_0(\omega_0^4q_2-
\omega^2(q_0-q_1p_2+q_2p_1)+q_0p_1-p_0q_1)}{\omega_0^4(q_1-q_2p_2)+\omega_0^2(q_0p_2-q_1p_1+p_0q_2)
-p_0q_0}\right)$$ and $\omega_0$ is a root of (11).}

\section*{\normalsize\bf 3. The normal form for system (4). Cyclical behavior.}

In this section we describe the direction, stability and the
period of the bifurcating periodic solutions of system (4). The
method we use is based on the normal form theory and the center
manifold theorem introduced by Hassard  [4]. Taking into account
the previous section, if $\tau=\tau_0$ then all roots of equation
(9) other than $\partialm i\omega_0$  have negative real parts, and any
root of equation (9) of the form $\lambda(\tau)=\alpha(\tau)\partialm
i\omega(\tau)$ satisfies $\alpha(\tau_0)=0$, $\omega(\tau_0)=\omega_0$ and
$\displaystyle\frac{d\alpha(\tau_0)}{d\tau}\neq0$. For notational convenience let
$\tau=\tau_0+\mu, \mu\in{\rm{I \! R}}$. Then $\mu=0$ is the Hopf bifurcation
value for equations (4).

Define the space of continuous real-valued functions as
$C=C([-\tau_0,0],{\rm{I \! R}}^4).$

In $\tau=\tau_0+\mu, \mu\in{\rm{I \! R}}$, we regard $\mu$ as the
bifurcation parameter. For $\Phi\in C$ we define a linear
operator:
$$L(\mu)\Phi=A\Phi(0)+B\Phi(-\tau)$$ where A and B are given by
(7) and a nonlinear operator $F(\mu, \Phi)=F(\Phi(0)$,
$\Phi(-\tau))$, where $F(\Phi(0), \Phi(-\tau))$ is given by (8).
By the Riesz representation theorem, there exists a matrix whose
components are bounded variation functions, $\eta(\thetaeta, \mu)$
with $\thetaeta\in[-\tau_0, 0]$ such that:
$$L(\mu)\Phi=\int\limits_{-\tau_0}^0d\eta(\thetaeta,\mu)\partialhi(\thetaeta),
\quad \thetaeta\in[-\tau_0,0].$$

For $\Phi\in C^1([-\tau_0, 0], {\rm{I \! R}}^{4})$ we define:
$${\cal A}(\mu)\Phi(\theta)=\left\{\begin{array}{ll} 

\displaystyle\frac{d\Phi(\theta)}{d\theta}, & \theta\in[-\tau_0,0)\\
\int\limits_{-\tau_0}^0d\eta(t,\mu)\partialhi(t), &
\theta=0,\end{array}\right.$$
$$R(\mu)\Phi(\theta)=\left\{\begin{array}{ll} 

0, & \theta\in[-\tau_0,0)\\
F(\mu, \Phi), & \theta=0.\end{array}\right.$$

We can rewrite (6) in the following vector form:
$$\dot u_t={\cal A}(\mu)u_t+R(\mu)u_t\eqno(12)$$
where $u=(u_1, u_2, u_3, u_4)^T$, $u_t=u(t+\thetaeta)$ for
$\thetaeta\in[-\tau_0, 0]$.

For $\Psi\in C^1([0,\tau_0], {\rm{I \! R}}^{*4})$, we define the adjunct
operator ${\cal A}^*$ of ${\cal A}$ by:
$${\cal A}^*\Psi(s)=\left\{\begin{array}{ll} 

-\displaystyle\frac{d\Psi(s)}{ds}, & s\in(0, \tau_0]\\
\int\limits_{-\tau_0}^0d\eta^T(t,0)\partialsi(-t), &
s=0.\end{array}\right.$$ We define the following bilinear form:
$$<\Psi(\theta),
\Phi(\theta)>=\bar\Psi^T(0)\Phi(0)-\int_{-\tau_0}^0\int_{\xi=0}^0\bar\Psi^T(\xi-\thetaeta)d\eta(\thetaeta)\partialhi(\xi)d\xi,$$
where $\eta(\thetaeta)=\eta(\thetaeta,0)$.

We assume that $\partialm i\omega_0$ are eigenvalues of ${\cal A}(0)$. Thus,
they are also eigenvalues of ${\cal A}^*$. We can easily obtain:
$$\Phi(\theta)=ve^{\lambda_1\theta},\quad \theta\in[-\tau_0, 0]\eqno(13)$$
where $v=(v_1, v_2, v_3, v_4)^T$,
$$v_1=1, v_2=-\displaystyle\frac{\beta(a_4+b_4e^{-\lambda_1\tau_0}-m\lambda_1)}{\lambda_1(\lambda_1-\beta\gamma_0\gamma_1)}, v_4=\displaystyle\frac{a_4+b_4e^{-\lambda_1\tau_0}}
{\lambda_1},$$
$$v_3=\displaystyle\frac{a}{\lambda_1+\deltalta}[\rho_{10}e^{-\lambda_1\tau_0}+\rho_{01}\displaystyle\frac{\beta(a_4+b_4e^{-\lambda_1\tau_0}-m\lambda_1)}{\lambda_1(\beta\gamma_0\gamma_1-\lambda_1)}]$$
is the eigenvector of ${\cal A}(0)$ corresponding to $\lambda_1=i\omega_0$
and
$$\Psi(s)=we^{\lambda_2s},\quad s\in[0,\infty)$$ where
$w=(w_1, w_2, w_3, w_4)$,
$$w_1=\displaystyle\frac{\lambda_2-\beta\gamma_0\gamma_1}{\alphapha a\rho_{01}}\displaystyle\frac{1}{\bar\eta}, w_2=\displaystyle\frac{1}{\bar\eta},
w_3=0, w_4=-\displaystyle\frac{\beta}{\lambda_2}\displaystyle\frac{1}{\bar\eta}$$
$$\eta=\displaystyle\frac{\lambda_1-\beta\gamma_0\gamma_1}{a\alphapha\rho_{01}}(1+\alpha b_1\displaystyle\frac{\lambda_1\tau_0e^{-\lambda_1\tau_0}\!-\!e^{-\lambda_1\tau_0}+1}{\lambda_1^2})+v_2-
\displaystyle\frac{\beta}{\lambda_1}
(v_4+b_4\displaystyle\frac{\lambda_1\tau_0e^{-\lambda_1\tau_0}\!-\!e^{-\lambda_1\tau_0}+1}{\lambda_1^2})$$
is the eigenvector of ${\cal A}(0)$ corresponding to $\lambda_2=-i\omega_0.$

We can verify that: $<\Psi(s), \Phi(s)>=1$, $<\Psi(s),
\bar\Phi(s)>=<\bar\Psi(s), \Phi(s)>=0$, $<\bar\Psi(s),
\bar\Phi(s)>=1.$

Using the approach of Hassard [4], we next compute the coordinates
to describe the center manifold $\Omegaega_0$ at $\mu=0$. Let
$u_t=u_t(t+\theta), \theta\in[-\tau_0,0)$ be the solution of equation
(12) when $\mu=0$ and
$$z(t)=<\Psi, u_t>,
\quad w(t,\theta)=u_t(\theta)-2Re\{z(t)\Phi(\theta)\}.$$

On the center manifold $\Omegaega_0$, we have:
$$w(t,\theta)=w(z(t), \bar z(t), \theta)$$ where
$$w(z,\bar z, \theta)=w_{20}(\theta)\displaystyle\frac{z^2}{2}+w_{11}(\theta)z\bar
z+w_{02}(\theta)\displaystyle\frac{\bar z^2}{2}+w_{30}(\theta)\displaystyle\frac{z^3}{6}+\dots$$
in which $z$ and $\bar z$ are local coordinates for the center
manifold $\Omegaega_0$ in the direction of $\Psi$ and $\bar\Psi$ and
$w_{02}(\theta)=\bar w_{20}(\theta)$. Note that $w$ and $u_t$ are real.

For solution $u_t\in \Omegaega_0$ of equation (12), since $\mu=0$, we
have:
$$\dot z(t)=\lambda_1z(t)+g(z, \bar z)\eqno(14)$$ where
$$\begin{array}{ll}
g(z, \bar z)& =\bar\Psi(0)F(w(z(t),\bar z(t), 0)+2Re(z(t)\Phi(0)))=\\
& =g_{20}\displaystyle\frac{z(t)^2}{2}+g_{11}z(t)\bar z(t)+g_{02}\displaystyle\frac{\bar
z(t)^2}{2}+g_{21}\displaystyle\frac{z(t)^2\bar z(t)}{2}+\dots\end{array}$$
where
$$g_{20}=F^1_{20}\bar w_1+F^2_{20}\bar w_2,
g_{11}=F^1_{11}\bar w_1+F^2_{11}\bar w_2, g_{02}=F^1_{02}\bar
w_1+F^2_{02}\bar w_2,\eqno(15)$$ with
$$\begin{array}{l}F_{20}^1=\alphapha a(\rho_{20}+2\rho_{11}v_2+\rho_{02}v_2^2), F_{11}^1=\alphapha a(\rho_{20}+\rho_{11}(\bar v_2+v_2)+\rho_{02}v_2\bar v_2), \\
F_{20}^2=\gamma_0\gamma_2v_2^2, F_{11}^2=\gamma_0\gamma_2v_2\bar
v_2, F_{02}^1= \bar F_{20}^1, F_{02}^2=\bar F_{20}^2,
\end{array}$$ and
$$g_{21}=F^1_{21}\bar w_1+F^2_{21}\bar w_2\eqno(16)$$ where
$$\begin{array}{l}
F_{21}^1= \alphapha a\rho_{20}(2w_{11}^1(0)+w_{20}^1(0))+\alphapha
a\rho_{11}(2w_{11}^2(0)+w_{20}^2(0)+2w_{11}^1(0)\bar
v_2 +\\
w_{20}^1(0)v_2)+ \alphapha a\rho_{02}(2w_{11}^2(0)v_2+w_{20}^2(0)\bar
v_2)+\alphapha
a(\rho_{30}+2\rho_{21}v_2+2\rho_{12}v_2^2+\\
+\rho_{03}v_2^2\bar
v_2 +\rho_{12}v_2^2+\rho_{21}\bar v_2)\end{array}$$
$$F_{21}^2=\gamma_0\gamma_2(2w_{11}^2(0)v_2+w_{20}^2(0)\bar
v_2)+\gamma_0\gamma_3v_2^2\bar v_2.$$

The vectors $w_{20}(\thetaeta)$, $w_{11}(\thetaeta)$ with
$\thetaeta\in[-\tau,0]$ are given by:
$$\begin{array}{l}
w_{20}(\thetaeta)=-\displaystyle\frac{g_{20}}{\lambda_1}ve^{\lambda_1\thetaeta}-\displaystyle\frac{\bar
g_{02}}{3\lambda_1}\bar v_e^{\lambda_2\thetaeta}+E_1e^{2\lambda_1\thetaeta}\\
w_{11}(\thetaeta)=\displaystyle\frac{g_{11}}{\lambda_1}ve^{\lambda_1\thetaeta}-\displaystyle\frac{\bar
g_{11}}{\lambda_1}\bar ve^{\lambda_2\thetaeta}+E_2\end{array}\eqno(17)$$
where
$$
E_1=-(A+e^{-\lambda_1\tau_0}B-2\lambda_1I)^{-1}F_{20},\quad
E_2=-(A+B)^{-1}F_{11},$$ where $F_{20}=(F_{20}^1, F_{20}^2,
F_{20}^3, 0)^T$, $F_{11}=(F_{11}^1, F_{11}^2, F_{11}^3, 0)^T$.

Based on the above analysis and calculation, we can see that each
$g_{ij}$ in (15), (16) are determined by the parameters and delay
from system (4). Thus, we can explicitly compute the following
quantities:
$$\begin{array}{l}
C_1(0)=\displaystyle\frac{i}{2\omega_0}(g_{20}g_{11}-2|g_{11}|^2-\displaystyle\frac{1}{3}|g_{02}|^2)+\displaystyle\frac{g_{21}}{2}\\
\mu_2=-\displaystyle\frac{Re(C_1(0))}{Re\lambda'(\tau_0)},
T_2=-\displaystyle\frac{Im(C_1(0))+\mu_2Im\lambda'(\tau_0)}{\omega_0},
\beta_2=2Re(C_1(0)).\end{array}\eqno(18)$$ and
$$\lambda'=\displaystyle\frac{\lambda(q_2\lambda^2+q_1\lambda+q_0)e^{-\lambda \tau}}
{3\lambda^2+2p_2\lambda+p_1+e^{-\lambda\tau}(2q_2\lambda+q_1-\tau(q_2\lambda^2+q_1\lambda+q_0))}$$
In summary, this leads to the following result:

 {\bf Theorem 3.} {\it In formulas (18), $\mu_2$
determines the direction of the Hopf bifurcation: if $\mu_2>0
(<0)$, then the Hopf bifurcation is supercritical (subcritical)
and the bifurcating periodic solutions exit for $\tau>\tau_0
(<\tau_0)$; $\beta_2$ determines the stability of the bifurcating
periodic solutions: the solutions are orbitally stable (unstable)
if $\beta_2<0 (>0)$; and $T_2$ determines the period of the
bifurcating periodic solutions: the period increases (decreases)
if $T_2>0 (<0)$.}

\section*{\normalsize\bf 4. Numerical example.}

\hspace{0.6cm} In this section we find the waveform plots through
the formula:
$$X(t+\theta)\!=\! z(t)\Phi(\theta)\!+\!\bar
z(t)\bar\Phi(\theta)\!+\!\displaystyle\frac{1}{2}w_{20}(\theta)z^2(t)+w_{11}(\theta)z(t)\bar
z(t)\!+\!\displaystyle\frac{1}{2}w_{02}(\theta)\bar z(t)^2+X_0,$$ where $z(t)$ is
the solution of (14), $\Phi(\theta)$ is given by (13), $w_{20}(\theta),
w_{11}(\theta), w_{02}(\theta)=\bar w_{20}(\theta)$ are given by (17) and
$X_0=(Y_0, r_0, K_0, M_0)^T$ is the equilibrium state.

For the numerical simulations we use Maple 9.5. We consider system
(4) with $a=0.38$, $\alphapha=0.96$, $\beta=1$, $\alphapha_1=0.5$,
$\alphapha_2=0.83$, $\gamma_0=1$, $d=0.1$, $s=0.3$, $r_2=0.003$,
$\deltalta=0.2$, $m=0.005$, $g=50$. The equilibrium point is:
$Y_0\!=\! 500$, $r_0\!=\! 0.03572181612$, $K_0\!=\! 675$,
$M_0\!=\! 33.06065092$. In what follows we consider two different
shares $\epsilon$ of delay tax revenues: $\varepsilon=0.3$ and
$\varepsilon=0.8$.

For $\varepsilon=0.3$ we obtain: $\mu_2\!=\!
1.654628706\cdot10^{-8}$, $\beta_2\!=\! 2.224294680\cdot10^{-9}$,
$T_2\!=\! 2.092652051\cdot10^{-9}$, $\omegaega_0\!=\! 0.6685954740$,
$\tau_0\!=\! 4.965007916$. Then the Hopf bifurcation is
supercri\-ti\-cal, the solutions are orbitally unstable and the
period of the solution is increasing. The wave plots are given in
the following figures:

\begin{center}
{\small \begin{tabular}{c|c|c} \hline
Fig.1.Waveplot $(t,Y(t))$&Fig.2.Waveplot $(t,r(t))$&Fig.3.Waveplot $(t,K(t))$\\&&\\
\cline{1-3} \epsfxsize=4cm

\epsfysize=4cm

\epsffile{fig1.eps} &

\epsfxsize=4cm

\epsfysize=4cm

\epsffile{fig2.eps} &

\epsfxsize=4cm

\epsfysize=4cm

\epsffile{fig3.eps}
\\
 \hline
\end{tabular}}
\end{center}

\begin{center}
{\small \begin{tabular}{c|c|c} \hline
Fig.4.Waveplot $(t,M(t))$&Fig.5.Waveplot $(t,I(t))$&Fig.6.Waveplot $(t,L(t))$\\&&\\
\cline{1-3} \epsfxsize=4cm

\epsfysize=4cm

\epsffile{fig4.eps} &

\epsfxsize=4cm

\epsfysize=4cm

\epsffile{fig5.eps} &

\epsfxsize=4cm

\epsfysize=4cm

\epsffile{fig6.eps}
\\
 \hline
\end{tabular}}
\end{center}

For $\varepsilon=0.8$ we obtain: $\mu_2\!=\!
-9.160756314\cdot10^{-8}$, $\beta_2\!=\!
-1.968119398\cdot10^{-8}$, $T_2\!=\! -1.608068638\cdot10^{-8}$,
$\omegaega_0\!=\! 0.8553440397$, $\tau_0\!=\! 3.918246696$. Then the
Hopf bifurcation is subcritical, the solutions are orbitally
stable and the period of the solution is decreasing. The wave
plots are given in the following figures:

\begin{center}
{\small \begin{tabular}{c|c|c} \hline
Fig.7.Waveplot $(t,Y(t))$&Fig.8.Waveplot $(t,r(t))$&Fig.9.Waveplot $(t,K(t))$\\&&\\
\cline{1-3} \epsfxsize=4cm

\epsfysize=4cm

\epsffile{fig7.eps} &

\epsfxsize=4cm

\epsfysize=4cm

\epsffile{fig8.eps} &

\epsfxsize=4cm

\epsfysize=4cm

\epsffile{fig9.eps}
\\
 \hline
\end{tabular}}
\end{center}

\begin{center}
{\small \begin{tabular}{c|c|c} \hline
Fig.10.Waveplot $(t,M(t))$&Fig.11.Waveplot $(t,I(t))$&Fig.12.Waveplot $(t,L(t))$\\&&\\
\cline{1-3} \epsfxsize=4cm

\epsfysize=4cm

\epsffile{fig10.eps} &

\epsfxsize=4cm

\epsfysize=4cm

\epsffile{fig11.eps} &

\epsfxsize=4cm

\epsfysize=4cm

\epsffile{fig12.eps}
\\
 \hline
\end{tabular}}
\end{center}

\section*{\normalsize\bf 5. Conclusions.}

From the analysis of the model with continuous time, it results
that the model accepts a limit cycle. The nature of the limit
cycle is given by the coefficients (18) which include the
parameters of the model. We establish the nature of the limit
cycle. Because the expressions of the parameter and the
coefficients (18) are difficult to analyze directly, the use of
Maple 9.5 was essential. The paper's results confirm that a series
of economical processes, where a variable with time delay
intervenes, have a limit cycle, thus, allowing a prediction
concerning the evolution of the model.

\end{document}