HOW INTEREST RATE INFLUENCES A BUSINESS CYCLE MODEL

. We study the eﬀect of interest rate on phenomenon of business cycle in a Kaldor-Kalecki model. From the information of the People’s Bank of China and the Federal Reserve System, we know the interest rate is not a constant but with remarkable periodic volatility. Therefore, we consider periodically forced interest rate in the model and study its dynamics. It is found that, both limit cycle through Hopf bifurcation in unforced system and periodic solutions generated by period doubling bifurcation or resonance in periodically forced system, can lead to cyclical economic ﬂuctuations. Our analysis reveals that the cyclical ﬂuctuation of interest rate is one of a key formation mechanism of business cycle, which agrees well with the pure monetary theory on business cycle. Moreover, this ﬂuctuation can cause chaos in a business cycle system.


1.
Introduction. Business cycle, one of the attractive economic phenomenon, plays an important role in economic decisions and market regulation [5,6,19]. Some mathematical models are proposed to investigate this phenomenon, while the earliest work can be traced back to Kaldor in 1940 [4,13]. Kaldor uses a differential equation system to study the trade cycle in steady-state economy, in which people's preferences and technological levels remain the same, while population and capital stocks may grow. After many years' development, the representative Kaldor-Kalecki business cycle model [16] was proposed Ẏ = α(I(Y, K) − S(Y, K)), where Y and K are national income and capital stock at time τ respectively, α > 0 is the adjustment coefficient in the goods market, δ ∈ (0, 1) is the depreciation rate of the capital stock. I(Y, K) and S(Y, K) are investment function and saving function respectively. The familiar forms of these two functions are where γ represents the propensity to save, β measures the sensitivity of investments respect to a variation in capital stock.
The dynamic behaviors of model (1) are widely studied such as stability, delay, bifurcations and stochastic sensitivity analysis [1,10,12,26]. In spite of these efforts, the effect of interest rate on investment has not been considered in this model. According to the classical Keynes's theory [14], interest rate has strong influence on investment and thus a financial system. To be specific, a rise (fall) in interest rates will shift the investment down (up). Based on this argument, we take the following investments function I(Y, K) = I(Y ) − βK − d.
(3) d(> 0) represents the variation induced by interest rate, which is directly proportional to interest rate. Meanwhile, by a general consensus among economists, function I(Y ) is believed to have the form of "S-shape" to variable Y [1,3,20]. I(Y ) is an increasing function of the national income Y , which cannot reach too high or too low values. In this paper we use the following form of I(Y ) I(Y ) = C 0 + nY 2 mY 2 + 1 (4) to model the "S-shape"(see Fig. 1), where C 0 , m, n are positive constants. In this case, by performing a new time transformation dτ = (mY 2 + 1)dt, the model (1) can be expressed as a planar polynomial system Ẏ = αnY 2 + α(mY 2 + 1)(C 0 − d − βK − γY ), It is well known that the operation of any economic system can not avoid being influenced by external periodic action, such as cyclical government regulation and seasonal changes [11]. Not only that, the data from the People's Bank of China (PBC) and the Federal Reserve System (FED) indicates that interest rate also shows cyclical fluctuation, see Fig. 2, which can be divided into three stages: growth, reduction and regrowth. Naturally, a question arise: How the periodically forced interest rate affects a business cycle system. Motivated by these, we consider parameter d as a function d(t) = d(1 + sinωt).
Then system (5) becomes non-autonomous, where ω and are the frequency and amplitude of the forced term, respectively. Without loss of generality, assume ω > 0 and 0 ≤ ≤ 1. We investigate bifurcations of equilibrium in unforced system ( = 0) and periodic solution in periodically forced system (0 < ≤ 1) respectively.
The aim of this paper is to study the influence of interest rate on dynamics of the business cycle system (5). We find that external periodic actions, such as the cyclical fluctuation of interest rate, can lead to business cycle and even chaos. The stable periodic solutions with different periods in periodically forced system, generated by period doubling bifurcations and resonance, can explain cyclical economic fluctuations. This can be seen as a viewpoint to support the pure monetary theory on business cycle proposed by Hawtrey [7][8][9] who suggests the interest rate shows major effect on business cycle. Moreover, the volatility of interest rate will give rise to business cycle with nearly a 9-years period, which is consistent with the Juglar cycle [2,17].

2.
Bifurcations of the unforced system. In this section, we will consider bifurcations of system (5) for case = 0 (unforced). For simplicity, we write the system as following where System (7) has for Y > 0 the same orbits as the polynomial system (5), and the equilibrium remains unchanged. We get the equilibria of the unforced system by It is easy to find from equations (8) that K = γY /δ, and the coordinate Y of equilibria for the unforced system should satisfy Obviously, the unforced system does not exist trivial equilibrium (Y = 0) and the number of equilibria is determined by the number of real roots of f (Y ) = 0. Here we choose not to present the exact expression of equation (9). The unforced system (7) always exists equilibrium for f (−∞) < 0 and f (+∞) > 0. Moreover, if d − C 0 < 0, then the unforced system has at least one positive equilibrium due to f (0) < 0.
We denote any equilibrium of the unforced system as E(Y 0 , K 0 ) where the coordinate Y 0 should satisfy (9). The Jacobi matrix of system (7) at equilibrium E can be given as Let λ be the eigenvalue of J(E), then we obtain the corresponding characteristic equation by straightforward calculation Clearly, E(Y 0 , K 0 ) is a hyperbolic equilibrium if conditions β +δ −α(I (Y 0 )−γ) = 0 and βγ + δγ − δI (Y 0 ) = 0 hold. For system (5) we have the following theorem.
Proof. We will compute center manifold for the fold bifurcation. When βγ + δγ − δI (Y 0 ) = 0, one can find the eigenvalues of Jocabi matrix of system (5) at equilibrium E(Y 0 , K 0 ) are If αβγ − δ(β + δ) = 0, one can find by straightforward calculation that V 1 , V 2 are two eigenvectors of system (5) We derive the normal form on one dimensional center manifold as follows. Let E(Y 0 , K 0 ) be the equilibrium of system (5) and satisfies the relation (8).
can bring E(Y 0 , K 0 ) to the origin, and we get a new system of coordinatesŶ ,K.
Under the translation T , we can diagonalize the linear part and obtain a system of coordinatesỸ ,K. Rewriting the system and dropping the "hats", it becomes Ẏ = a 20 Y 2 + a 11 Y K + a 02 K 2 + o(|Y, K| 2 ), where Here, we do not present other coefficients of system (13) and they can be obtained by straightforward calculation.
In the following, we give the form of one dimensional center manifold. For Y ∼ 0, there exist a center manifold System (13) reduced on the one dimensional center manifold (14) is given bẏ Hence, E is a saddle-node of codimension 1 if a 20 = 0. It is a cusp if a 20 = 0, then we should compute the coefficient of term Y 3 .
We will drive the normal form of the Hopf bifurcation of system (5). If β + δ − α(I (Y 0 ) − γ) = 0 and αβγ > δ(β + δ) hold, we can get the eigenvalues of Jocabi matrix of system (5) λ ± = ±iµ, where µ = (1 + mY 2 0 ) (αβγ − δ(β + δ)), i is an imaginary unit. Let the eigenvectors corresponding to the eigenvalues be U 1 ± iU 2 , where U 1 and U 2 are real vectors. By simple calculations we obtain Bring E(Y 0 , K 0 ) to the origin by translation (11), then we obtain system (12). Define Under the translation P we get a system of coordinatesỸ ,K. Dropping the "hats", we obtain the normal form of system (5) at the Hopf bifurcation point where a ij , b ij (i, j = 1, 2, 3) are the coefficients of system (15), we choose not to present them here. In order to study the stability and direction of bifurcating periodic solutions, we should compute the first Lyapunov coefficient l 1 of the Hopf Bifurcation. By the formula in [21], we can get l 1 Generally, we have the following three cases for system (5) (1) If l 1 < 0, the system undergoes supercritical Hopf bifurcation; (2) If l 1 > 0, the system undergoes subcritical Hopf bifurcation; (3) If l 1 = 0, the system undergoes degenerate Hopf bifurcation.
Here we cannot give the concise explicit condition to determine the sign of l 1 due to its complexity. Therefore we show that l 1 can have alternate sign or vanish by simulation as follows.
Choose α = 3, β = 0.6, δ = 0.3, m = 3, n = 5.6, C 0 = 0.6, we can obtain bifurcation diagram of system (5) in d − γ plane, see Fig. 4. The Hopf bifurcation curve (green curve) passes though two degenerate Hopf bifurcation points, labeled as DH 1 and DH 2 . The first Lyapunov coefficient of these two degenerate Hopf bifurcation points is l 1 = 0, and the codimension of these two point is at least 2. The curve between the two degenerate Hopf bifurcation points (DH 1 and DH 2 ) is supercritical Hopf bifurcation curve where the corresponding first Lyapunov coefficient is l 1 < 0. The rest of the green curve is subcritical Hopf bifurcation curve where l 1 > 0. The fold bifurcation curve (red curve) of system (5) forms a closed region with three cusp, labeled as CP 1 , CP 2 and CP 3 . The three points at the cusp of the fold bifurcation curve are cusp bifurcation points at which the coefficient of quadratic term of the reduced system on the center manifold is 0. The points BT 1 and BT 2 are the Bogdanov-Takens bifurcation point of codimension 2, which are tangent points of Hopf bifurcation curve and fold bifurcations, and the corresponding eigenvalues at these points are λ 1,2 = 0. Figure 5(a) is phase portrait of system (5) near a supercritical Hopf biurcation for α = 3, β = 0.6, γ = 1.175, d = 0.195, δ = 0.3, m = 3, n = 5.6, C 0 = 0.6. The red curve in Fig. 5(a) is a stable limit cycle bifurcate from the Hopf bifurcation point (0.151, 0.582). One can see that as t → +∞, all of the trajectories (green curve and blue curve) spiral to the stable limit cycle (red curve). Figure 5(b) is phase portrait of system (5) near a subcritical Hopf biurcation for α = 3, β = 0.6, γ = 0.66, d = 0.195, δ = 0.3, m = 3, n = 5.6, C 0 = 0.6. One can find there are two limit cycles (red curve) exist in system (5), the large one is stable and the small one is unstable. As as t → +∞, the trajectories (green curve and blue curve) starting near the large cycle spiral to the stable limit cycle (large one). Since the Hopf bifurcation occurs on the 2-dimensional center manifold, so there is an unstable limit cycle on the center manifold. The trajectory (black curve) starting inside the unstable cycle spirals inward to the stable focus (0.849, 1.868). Moreover, we give the corresponding time series of the stable and unstable limit cycles generated by Hopf bifurcation, see Fig. 6. Figure 6(a) is time series of the stable limit cycle generated by the supercritical Hopf bifurcation where the initial value is selected as (0.3, 0.75). Figure 6(b) and 6(c) are time series of limit cycles when system (5) undergoes supercritical Hopf bifurcation where the initial values are selected as (2, 2.3) and (0.9136, 1.8439) respectively. Proof. For system system (7), the Jacobi matrix at equilibrium E can be written as hold, one can find the characteristic equation at E has the form λ 2 = 0.
Following we will drive the normal form of the Bogdanov-Takens bifurcation at E(Y 0 , K 0 ) of system (5). Bring E(Y 0 , K 0 ) to the origin by translation (11), then we obtain system (12). For β + δ = 0, the generalized eigenvectors of λ = 0 are , then under the relation f (Y 0 ) = 0 and the linear transformation we get a system of coordinatesỸ ,K. Dropping the "hats", system (5) becomes Here Using the near-identity transformation and rewrite u, v into Y, K, we obtain whereL 20 = L 20 andL 11 = L 11 + 2H 20 . IfL 11 = 0, we make a change of coordinates and time and preserve the orientation by time then system (5) is topologically equivalent to the normal form (16).
3. Bifurcations of the periodically forced system. In this section, we investigate the dynamical behaviors of Kaldor-Kalecki model with periodically forced interest rate. Some bifurcation curves and phase diagrams are given by us to reveal the complex dynamics of the periodically forced system. We consider system (5) with periodically forced interest rate d(t), and the system can be given as where d(t) = d(1 + sinωt). Here d, and ω are the mean value, amplitude and frequency of the oscillation of interest rate respectively. The value of the frequency ω can be selected according to reality. Without loss of generality, we assume ω = 2π which implies the period is 1 for the oscillation. In order to ensure d(t) > 0 we assume 0 ≤ ≤ 1.
Notice that the forced system is nonautonomous, which makes plenty of results on the autonomous system no longer available. In order to study the dynamics of a periodically forced system, we transform the system into a higher dimensional autonomous system by augmenting the forced system with a decoupled pair of ODEs which have an oscillatory solution of the form required for the forcing. System (19) can be done by adding a nonlinear oscillator with the period-1 forcing as one of the solution components. We take such a form of nonlinear oscillator which has the asymptotically stable solution v = sin2πt, w = cos2πt. Therefore, system (19) can be transform into a autonomous four dimensional system, with initial conditions v(0) = 0 and w(0) = 1.
In this way the equilibrium E(Y 0 , K 0 ) of the two dimensional unforced system changes as the periodic solutionĒ(Y (t), K(t), sin2πt, cos2πt) of the four dimensional system (21). We can study whether the periodic solution can survive and whether the bifurcations of the equilibrium of the unforced system can be extended to the forced system as bifurcations of periodic solutions when = 0. We use the software package AUTO to study bifurcations of the continuous four dimensional system by investigating the Poincaré map of system (21). For our case, first return map of the continuous four dimensional system can be defined as P : (Y (0), K(0), v(0), w(0)) −→ (Y (1), K(1), v(1), w(1)).
As we know, the stable (unstable) fixed points of the kth iterate of the map correspond to the stable (unstable ) periodic solutions with period k of system (21). More details about the Poincaré map and notations in obtained bifurcation diagrams can be found in [18,22,23].
It is obvious that the Poincaré map P has four multipliers for the dimension of system (21) is four. Two of the four multipliers are always constants in parameter space, due to the asymptotically stable periodic solution generated by the nonlinear oscillator (20). Denote the unchangeable multipliers as |µ 3 | = 1, |µ 4 | = c, c is a constant. We focus on the variation of multipliers µ 1 and µ 2 in parameter region, and give the following illustrations in this paper. If |µ 1 | > 1, |µ 2 | > 1, the corresponding fixed point is locally unstable and it is called a repelling point. If |µ 1 | < 1, |µ 2 | < 1, the corresponding fixed point is locally asymptotically stable and it is called a attracting point. If |µ 1 | > 1, |µ 2 | < 1 (or |µ 1 | < 1, |µ 2 | > 1), the corresponding fixed point is locally unstable and it is called a saddle. For convenience, we regard the periodic solutions of (21) as stable ( saddle or repelling) periodic solutions, if the corresponding fixed point of map P is stable (saddle or repelling). We use the notations h (k) , f (k) , t (k) to denote the Hopf (Neimark-Sacker) bifurcation curve, flip (Period doubling) bifurcation curve and fold (Tangent) bifurcation curve   respectively, which are generated by the fixed point with period k. Some other notations, A i , B i , C i , D i (i ∈ N ) are used to denote strong 1 : 1 resonance, strong 1 : 2 resonance, strong 1 : 3 resonance, strong 1 : 4 resonance, respectively. Case 1. Taking α = 3, β = 0.6, γ = 0.4, δ = 0.3, m = 3, n = 5.6, C 0 = 0.2 the unforced system (5) undergoes two subcritical Hopf bifurcation at d = d + and d = d − (d + > d − ) when d is varied. For this group of parameter values, the unforced system generates two unstable limit cycles. The asymptotic period of the cycle (evaluated numerically) is T = 2.343 and T = 9.373 when d approaches d + and d − respectively. The bifurcation diagrams of the forced system in the ( − d) planes are given in Fig. 7.
In Fig. 7(b) (or partial enlargements Fig. 7(c)), the points H 1 and H 2 on daxis correspond to two subcritical Hopf bifurcations in the unforced system (5) and they are roots of the curve h  respectively. A flip bifurcation curve f (1) passes though point B 1 , a 1:2 resonance, which is the terminal of curve h 1 1 and the staring point of curve h (2) . The curve h (2) contains other two codimension two bifurcation points, a 1:3 resonance C, a 1:4 resonance D, and terminates at a 1:2 resonance B 2 . The curve f (1) also passes though another codimension two bifurcation point P , which is the staring point of curve t (2) .   A pair of period-one fixed points will appear in region 1, one is saddle and another is unstable point, when the fold bifurcation curve t (1) 1 is crossed from above to region 1. The unstable fixed point will change as a saddle fixed point and a pair of unstable period-two fixed points will appear if the upper branch of f (1) (above point B 1 ) is crossed from region 1 to region 3. The generated saddle fixed point will become a stable fixed point and a pair of saddle period-two fixed points will appear if curve f (1) is crossed from region 3 to region 5. That is to say there are two pair of period-two fixed points in region 5. These two pair of period-two fixed points will collide on t (2) then disappear while t (2) is crossed from region 5 to the region 2. The unstable fixed point will become stable and an unstable closed invariant curve appears if curve h 1 is crossed to the below. The pair of unstable period-two fixed points become stable and an unstable closed invariant curve appears when the curve h (2) is crossed from above the below; they change as a pair of saddle fixed points and period-four fixed point will appear if the upper branch of curve f (2) is crossed from region 3 to region 4. And the pair of saddle period-two fixed points change as stable and another period-four fixed point appear while the curve f (2) is crossed from region 4 to the below. Besides, flip bifurcation curve f (4) , f (8) , ... , exist and this cascade of period doublings results in strange attractors. On the other hand, the unstable closed invariant curve of Poincaré map can be destroyed via homoclinic structure near 1:2 resonance giving rise to chaotic solutions.
In Fig. 8, the points T 1 and T 2 correspond to fold bifurcations in the unforced system, which are roots of curve t 1 passed though a 1:3 resonance C and terminates at a 1:2 resonance B, on curve f 1 , which is root of the curve h (2) . h (2) starts from B then terminates at a 1:1 resonance A on curve t (2) 3 which terminates at another codimension two bifurcation point P on flip bifurcation curve f (1) (see Fig.8(b)). At point T 3 the period of the unstable cycle is 2 and two branches of period two fold bifurcation curves, t  Fig.  7, so we do not repeat again. A pair of unstable period-two fixed points will appear when the upper branch curve f (1) is crossed from region 1 to the below. It will change as saddle fixed points and period-four fixed point will appear when the curve f (2) is crossed from the above to region 4. The pair of saddle fixed points become unstable when curve f (2) is crossed from region 4 to the below. Two pairs of period-two fixed point, the another pair of them is generated by the lower part of curve f (1) , will collide on t  is crossed to the below, and they will disappear due to a fold bifurcation curve (not given) which is very close to the lower branch of curve f (2) . Case 3. Taking α = 5, β = 0.6, γ = 0.85, δ = 0.3, m = 3.211, n = 3.368, C 0 = 1 the unforced system (5) undergoes two Supercritical Hopf bifurcation at d = d + and d = d − (d + > d − ) when d is varied. For this group of parameter values, the unforced system generates two stable limit cycles. The asymptotic period of the cycle (evaluated numerically) is T = 3.738 and T = 2.292 when d approaches d + and d − respectively. The bifurcation diagrams of the forced system in the ( − d) plane are given in Fig. 9.
In Fig. 9, the points H 1 and H 2 on the d-axis correspond to the supercritical Hopf bifurcation and are the root of curve h  passed though a 1:3 resonance C and terminates at a 1:2 resonance B 1 , on curve f 1 , which is the start point of the curve h (2) . h (2) starts from B 1 then terminates at a another 1:2 resonance B 2 where curve h (1) 2 terminates at it. f (1) and f (2) are flip bifurcation curves of period-one fixed point and period-two fixed point respectively. A stable period-one fixed point changes as unstable and a stable closed invariant curve appears in region 2 when the curve h 2 ) is crossed from region 1 (region 3) to region 2. On the other hand, the stable fixed point will become a saddle fixed point and a pair of stable period-two fixed points if the upper branch of f (1) is crossed from region 1 to region 4. The generated unstable fixed point will also become saddle fixed point if f (1) is crossed from region 2 to region 6. The saddle fixed point will change as stable and another pair of period-two fixed points appear if the lower branch of f (1) is crossed to region 3.
The pair of stable period-two fixed points will change as unstable and a stable closed invariant curve appears when h (2) is crossed from region 4 to region 6; And they become a pair of saddle fixed points and period-four fixed point appear if upper branch of curve f (2) is crossed from region 4 to region 5; The pair of saddle fixed points change as unstable and another period-four fixed point appear if f (2) is crossed from region 5 to the below. Case 4. Taking α = 5, β = 0.6, γ = 0.8, δ = 0.3, m = 3.211, n = 3.368, C 0 = 1 the unforced system (5) undergoes a supercritical Hopf bifurcation and a subcritical Hopf bifurcation at d = d + and d = d − (d + > d − ) respectively. For this group of parameter values, the unforced system generates two limit cycles, one is stable and the other is unstable. The asymptotic period of the cycle (evaluated numerically) is T = 3.915 and T = 2.735 when d approaches d + and d − respectively. The bifurcation diagrams of the forced system in the ( − d) plane are given in Fig. 10.
In Fig. 10, the points H 1 ( H 2 ) on the d-axis correspond to the supercritical (subcritical) Hopf bifurcation and are the root of curve h passed though a 1:3 resonance C 1 and terminates at a 1:2 resonance B 1 on curve f 1 , which is the start point of the curve h 2 ) passes though two codimension two bifurcation points,1:4 resonance D 2 (D 1 ) and 1:3 resonance C 3 (C 2 ), then terminates at a 1:2 resonance B 3 (B 4 ) on curve f (2) . Curves f (1) and f (2) are flip bifurcation curves of period-one fixed point and period-two fixed point respectively.
A stable period-one fixed point changes as unstable and a stable closed invariant curve appears in region 2 when the curve h (1) 1 is crossed from region 1 to region 2. And the generated unstable fixed point becomes stable and a unstable closed invariant curve appears in region 2 when the curve h (1) 2 is crossed from region 2 to region 3. On the other hand, the stable fixed point will become a saddle fixed point and a pair of stable period-two fixed points if the upper branch of f (1) is crossed from region 1 to region 4. The generated unstable fixed point will also become saddle fixed point if f (1) is crossed from region 2 to region 6. The saddle fixed point will change as stable and another pair of period-two fixed points appear if the lower branch of f (1) is crossed from region 7 to region 3.
The pair of stable period-two fixed points will change as unstable and a stable closed invariant curve appears when h (2) 1 is crossed from region 4 to region 6. The generated pair of unstable period-two fixed points will become stable and a unstable closed invariant curve appears when h In the periodically forced system (19), the equilibrium of system (5) becomes periodic solution of period T = 1 due to the adding nonlinear oscillator with frequency ω = 2π. Periodic forcing can induce different bifurcations of periodic solution, fold bifurcations, period-doubling bifurcations and torus bifurcations, which can lead to the appearance of various solutions.
In the following we will present periodic solutions with different periods where the parameter values, α, β, δ, m, n, C 0 , are the same as them in Case 3. Take α = 5, β = 0.6, δ = 0.3, γ = 0.85, m = 3.211, n = 3.368, C 0 = 1, and select parameters γ, d, as free parameters. Fig. 11(a) is a stable period-two orbit for d = 0.3, = 0.8, we can find from the time series that the period of this solution is 2. When period doubling (flip) bifurcation of period-two orbit occurs, period-two orbit changes the stability and period-four orbit appears. Fig. 11(d) is a stable period-four orbit for d = 0.18, = 0.8, where the period of the solution is 4. Period-eight orbits also exist in some regions of bifurcation diagrams as a result of bifurcation of period-four orbit. Figure 11(g) shows a stable torus for d = 0.2, = 0.1, which arises from a torus bifurcation. The period of the torus cannot be judged from the corresponding time series for it has an infinite period. Thus, we give the Poincaré section of the torus, see Fig. 11(h).  FIG. 13(b). The Largest Lyapunov exponents is calculated based on the algorithm in Ref [24]. Obviously, the spectrums show that the time series in Fig. 12(c) and Fig. 12  attractors due to the corresponding largest Lyapunov exponent of the time series greater than zero. Moreover, the spectrums also indicate there exist some other chaotic areas in parameter spaces of system (19). On the other hand, some ribbonlike structures with self-similarity in the diagram of Poincaré section also indicate that the corresponding attractor shows chaotic behavior, and it is easily found in Fig. 12(c) and Fig. 12(f). Figure 14 is the bifurcation diagram in (Y − d) plane where = 0.5 in Fig. 14(a) and = 0.2 in Fig. 14(b). Period doubling (flip) bifurcations can be seen intuitively in these two subgraphs, which leads to the appearance of periodic solutions with different periods in system (19) (see Fig. 12). The dense areas of points in Fig.  14 indicate there exist some attractors in this area, and the attractors should be torus or chaos. It is easy to find three lines exist in some areas in Fig. 14(b)   implies some stable period-three solutions (see Fig. 15(a)), generated by the strong 1:3 resonance C (see Fig. 9), exist in this region. In addition, the period-five solutions also exist in some parameter regions, which may be generated by other resonance. These periodic solutions can also undergo different bifurcations, and cascade of period doublings will lead to strange attractors.  attractors, through torus destruction and cascade of period doublings, can be detected in the periodically forced system. Some scholars ignored the effect of external periodic actions and attributed the phenomenon of cyclical economic fluctuations to limit cycle generated by Hopf bifurcation in the business cycle model. Our analysis reveals that the cyclical fluctuation of interest rate can also lead to the appearance of business cycle in a economic system. Moreover, the period of the generated business cycle is multiples of external period due to period doubling bifurcation and resonance. That is, if the period of external periodic actions is T , then the period of generated business cycle will be T, 2T, 3T, ... In 1862, Clément Juglar [2,17] suggested that there exists the business cycle of 9 to 11 years in economic system, and it is called the Juglar cycle. Generally, a complete business cycle can be divided into three stages by its trend: growth, reduction and regrowth. From Fig. 2 one can know that the period of interest rate of FED and PBC is near 9 years (from 1986.9 to 1995.4, from 1995.4 to 2004.1 in Fig. 2(a), from 1998.7 to 2008.12, from 2002.2 to 2011.7 in Fig. 2(b)). Thus we know from the bifurcation results in section 3 that there are period solutions of 9 years in the periodically forced business cycle system, which agrees with the Juglar cycle. In 1935, Kondratieff [15] analyzed a large number of economic statistics from different countries and found there was business cycle of 54 years in economic system. In later years, Schumpeter [25] developed the theory of Kondratieff and suggested there exists some long-term cycles of 48-60 years which can be divided into 6 mid-term cycle of 9 to 11 years. Actually, the period-6T solutions generated by period doubling bifurcation in the periodically forced system (19) can explain the long-term cycle in economic system. The period of a period-6T solutions is near 54 years because the period of the interest rate is near T = 9 (years).