Hypercrater Bifurcations , Attractor Coexistence , and Unfolding in a 5 D Model of Economic Dynamics

Complex dynamical features are explored in a discrete interregional macrodynamic model proposed by Asada et al., using numerical methods. The model is five-dimensional with four parameters. The results demonstrate patterns of dynamical behaviour, such as bifurcation processes and coexistence of attractors, generated by high-dimensional discrete systems. In three cases of two-dimensional parameter subspaces the stability of equilibrium region is determined and its boundaries, the flip and Neimark-Hopf bifurcation curves, are identified by means of necessary coefficient criteria. In the first case closed invariant curves CICs are found to occur through 5D-cratertype bifurcations, and for certain ranges of parameter values a stable equilibrium coexists with an unstable CIC associated with the subcritical bifurcation, as well as with an outer stable CIC. A remarkable feature of the second case is the coexistence of two attracting CICs outside the stability region. In both these cases the related hysteresis effects are illustrated by numerical simulations. In the third case a remarkable feature is the apparent unfolding of an attracting CIC before it evolves to a chaotic attractor. Examples of CICs and chaotic attractors are given in subspaces of phase space.


Introduction
We consider a 5D nonlinear discrete time model of the economic transactions between two regions with fixed exchange rates.The model is an extension of the original idea of Kaldor 1 which may be relevant, for example, to the analysis of the economic interdependency between two countries under currency integration such as two countries of the Eurozone.The two regions interact economically through trade and capital movement.We note that an analogous 5D interregional model in continuous time was studied in Asada 2 and also in Malik and Zimka 3 .The present model was first studied in Asada et al. 4 , mainly by analytical methods, while the corresponding model with flexible exchange rates was presented and explored numerically in Asada et al. 5 .Other two-regional Kaldorian models of business cycles, based on trade interaction between the regions, have been studied by Lorenz 6 and Puu 7 .
In a previous paper see 8 , hereafter called paper I we explored this model numerically and determined a trade threshold for the emergence of business cycles.In the present paper we pursue further our numerical exploration to reveal different aspects of the model dynamics.One of the aims is to illustrate the feasibility and effectiveness of the numerical approach for dynamical systems of high dimensionality and several parameters.
The present exploration concerns three characteristic cases of parameter ranges.In each case we determine numerically the stability region in parameter space and identify the flip and Neimark-Hopf bifurcation curves as parts of the boundary of this region.We also compute bifurcation and Lyapunov exponent diagrams providing information on the occurring closed invariant curves CICs and on the type of bifurcation responsible for their occurrence.
The specific aim here is to illustrate certain remarkable patterns of asymptotic dynamical behaviour that occur in this model, which to the authors' knowledge have not previously been presented for high-dimensional systems.Indeed, in the first case considered CICs are found to occur through "catastrophic" 5D-crater-type bifurcations, as a result of which there exists a region of parameter values of the model for which a stable equilibrium of the economy coexists with an attracting CIC, as well as with the repelling CIC characteristic of the subcritical bifurcation.We explore the event of this bifurcation process in some detail and illustrate the associated hysteresis effect by numerical simulations.In the second case CICs first appear by supercritical bifurcations, and a remarkable feature is the coexistence and discontinuous succession of attracting CICs, occurring outside the stability region.Again, we explore this case of attractor coexistence and illustrate the associated hysteresis effect.In the third case CICs similarly appear by supercritical bifurcations, and a remarkable feature in this case is the unfolding of an attracting CIC.Examples of fluctuations corresponding to CICs and chaotic behaviour of the economic variables are presented in sections of phase space.
The paper proceeds as follows.In Section 2 we present briefly the equations of the model.In Section 3 we present the position of equilibrium and briefly discuss the stability region in sections of the parameter space and the boundary curves of that region, namely the flip and Neimark-Hopf bifurcation curves.The results of our exploration in the three cases considered are presented, in the form of stability region and bifurcation and Lyapunov exponent diagrams, and discussed in Sections 4 to 6. Section 7 summarizes and concludes.

Model Equations
We consider the 5D model of interregional Kaldorian dynamics under fixed exchange rates proposed by Asada et al. 4 .The equations of the model in the form explored in paper I are the following five nonlinear difference relations: where we have abbreviated

2.2
Total nominal money supply of the two regions is fixed at M 1 t M 2 t M 600, and the function f is a particular case of the Kaldorian S-shaped dependence of the investment function on income given by The subscript i 1, 2 is the index number of a region and t denotes the time period.The price levels are considered fixed, and the meanings of the system variables are as follows.Y i is the real regional income, K i the real physical capital stock, and M i the nominal money stock.The parameters of the model are α i > 0 the adjustment speed of the goods market of each region, β > 0 the degree of capital mobility, and δ the degree of interregional trade 0 ≤ δ ≤ 1 .For a fuller description of the model and its economic foundations we refer to paper I.

Equilibrium Stability and Bifurcations
It is easily found that the system has a unique equilibrium at

3.1
All equilibrium values except K * 2 are positive, for β > 0 and 0 ≤ δ ≤ 1.K * 2 is also positive and we have a positive equilibrium if Stability of the equilibrium is determined by the roots of the characteristic polynomial of the Jacobian of the mapping 2.1 : and for stability all roots of P 5 λ must be inside the unit circle in the complex plane.Our basic tool for the numerical determination of the region of stability is the grid-search technique.This involves computing the characteristic polynomial and its roots at the node points of a dense grid in a two-dimensional space of parameters and plotting the points at which the equilibrium is stable.
The flip bifurcation curve is represented by the coefficient relation as a necessary condition, while a coefficient relation that must be satisfied as a necessary condition by the Neimark-Hopf bifurcation curve can be obtained by requiring that two of the roots of the characteristic polynomial have unit product see paper I points of a dense grid : Each one of 3.4 and 3.5 represents a relation among the model parameters α 1 , α 2 , β, δ.
For fixed values of two of the parameters it represents a relation between the remaining two parameters providing one or more curves in the plane of the remaining parameters.Segments of the curves arising from 3.4 form the parts of the boundary of the stability region that are flip bifurcation curves, and similarly segments of the curves arising from 3.5 form the parts of the boundary that are Neimark-Hopf bifurcation curves.In such a parameter plane CICs, the discrete equivalent of cyclical behaviour, are generated locally by bifurcation on a Neimark-Hopf bifurcation curve.
The bifurcation is established theoretically by the Neimark-Hopf bifurcation theorem for discrete systems see 9, 10 .It is employed successfully for the prediction of bifurcations in the present high-dimensional system under the same assumptions regarding the critical pair of eigenvalues of the Jacobian at the critical value of the bifurcation parameter, with the additional plausible requirement that all other eigenvalues satisfy the stability criterion |λ| < 1, the case of interest being the occurrence of CICs when stability of equilibrium is lost.
In the case of a supercritical bifurcation an attracting CIC is born out of the equilibrium when the latter loses stability.In the case of a subcritical bifurcation, on the other hand, equilibrium loses stability because a repelling CIC collapses on it.Both these bifurcations are local.The subcritical bifurcation, however, is often preceded by a global bifurcation in which an attracting CIC of finite amplitude is also generated inside the stability region.In this case we have coexistence of two attractors of different type, the equilibrium point itself and a stable CIC, with the unstable CIC also present.This process begins with the generation of the two CICs inside the stability region and ends when equilibrium loses stability on the bifurcation curve.Outside the stability region only the stable outer CIC remains.This twophase bifurcation process is sometimes called a crater bifurcation see 11 .It can be thought of as a catastrophic event occurring at the loss of equilibrium stability, due to the lack of continuity of the evolution of the steady state of the economy involved in the process.Note that as stated in Benhabib and Miyao 12 both supercritical and subcritical orbits seem plausible from an economic point of view.In economic dynamics, however, it is usually the supercritical stable CICs that are seen as stylized business or growth cycles.
In the next sections we proceed to discuss some interesting dynamics in three specific cases of parameter ranges of the model.

Hypercrater Bifurcation, Coexistence of Point and CIC Attractors
In this section we assume equal speeds of adjustment of the goods markets of the two regions α 1 α 2 α and adopt the value α 1 for the common speed of adjustment.We are then left with the two parameters β and δ for the model.This is the case exhibiting most clearly the economic interaction of the regions, as it leaves the regional interaction parameters free for variation enabling one to explore the influence of regional interaction on the dynamics of the economy.The stability region in this case as found by our grid search technique is shown in Figure 1.All branches of the curves produced by relations 3.4 and 3.5 are included but it is easy to see the relevant segments forming the boundary curves of the stability region.The curve representing condition 3.2 for positive equilibrium is also drawn in as a continuous line at the bottom of the diagram .The region's upper right-hand "corner" is at δ 0.5633, β 1.0863.For the economic interpretation of stability diagrams in the δ, β parameter plane we refer to paper I.The bifurcation diagrams of the regional incomes Y 1 and Y 2 are shown in Figure 2 for β 0.8, as computed with starting values very near to the equilibrium.The sudden appearance of finite amplitude at the value of the bifurcation parameter δ ∼ 0.568, at which equilibrium loses stability, suggests that the stability loss is caused by a subcritical Neimark-Hopf bifurcation and that the finite size attractor appearing exists before the loss of equilibrium stability.
In Figure 3 we present the results of some numerical experiments aimed at obtaining numerical-visual evidence of the existence of the unstable CIC, characteristic of the subcritical bifurcation, which causes equilibrium to lose stability by collapsing on it at the critical value of the bifurcation parameter.Each bifurcation diagram, and corresponding Lyapunov exponent diagram, of Figure 3 has been computed with the same starting values, equal to the equilibrium values for all variables, except Y 2 the starting value of which is marked in each case by a dot.
It is seen that the larger the distance of the starting Y 2 value is from its equilibrium value, the earlier the CIC appears in the bifurcation and Lyapunov exponent diagram.The stable CIC appears first earliest at δ ∼ 0.535 for a starting distance of 10 or more units in Y 2 , in the first diagram.In the second diagram for a starting distance of 8 units in Y 2 the stable CIC appears later at δ ∼ 0.546 , and in the third diagram for a starting distance of 6 units in Y 2 it appears even later at δ ∼ 0.559 .Finally, in the last diagram for near-zero starting distance in Y 2 the stable CIC appears when equilibrium loses stability at δ ∼ 0.568.
Repeating the experiment for several starting distances in Y 2 , above and below the equilibrium value heavy line , we obtain the curves shown in all bifurcation diagrams of Figure 3.For starting distances in Y 2 between these curves orbits converge to the equilibrium, and for starting distances in Y 2 above the upper or below the lower curve orbits converge to the outer stable CIC.These curves provide the sought visualization of the Y 2 amplitude of the unstable CIC which causes equilibrium to lose stability by collapsing on it at the critical value of the bifurcation parameter.It is established by these results that in the interval of the bifurcation parameter values 0.535 < δ < 0.568 the stable equilibrium coexists with an inner unstable CIC and with an outer stable CIC, as in the crater bifurcation case of two-dimensional systems.For a high-dimensional system as the present one, we call this a "hypercrater" bifurcation process.The outer stable CIC is shown for two values of δ in Figure 4.The subcritical-crater bifurcation process has interesting implications from the mathematical as well as the economic point of view.These have been studied extensively in the case of two dimensions see, e.g., 11, 13-15 .The mathematical interest relates to the study of the first phase of the bifurcation process in which the equilibrium is still stable and the attracting CIC coexists with the repelling CIC.From the economic point of view the interest relates to the interpretation of the dynamics involved in the case of attractor coexistence.When the stable equilibrium coexists with the attracting CIC, the coexisting repelling CIC bounds the basin of attraction of the equilibrium and ensures that small shocks do not affect the dynamical behaviour of the system.However, large shocks may lead to permanent fluctuations.In economic policy terms the hysteresis effect, which may be present in cases of attractor coexistence, creates uncertainty since adjustment of the parameter by policy cannot guarantee a return to stable equilibrium.
This effect has been illustrated for an important two-dimensional economic model in Agliari 16 .Here we illustrate the same effect in the case of the present 5D model for α 1, β 0.8. Figure 5 shows the trajectory under two shocks applied to the system by means of sudden changes of the parameter δ parametric shocks .Initially the value of the bifurcation parameter is δ 0.55 for which the equilibrium is stable and coexists with a stable CIC see Figure 2 .We start the iteration near the equilibrium value, and after a transient phase the value of Y 2 quickly settles down to the equilibrium value Y * 2 ∼ 178.03.At t 200 a shock on the system is applied, by changing the parameter value to δ 0.6 for which the equilibrium is unstable.The trajectory begins to fluctuate and after a transient phase attains regular large fluctuations.It has been attracted by the stable CIC.
A second shock is applied, at some later time, by restoring the parameter to its previous "stable" value.There are now two possibilities as follows.
If at the time of application of the second shock the variables have values sufficiently near the stable equilibrium, then the trajectory is attracted by the equilibrium again.In our example this happens with the second shock applied at t 515 as shown in Figure 5 a .If, however, at the time of application of the second shock the variables have values sufficiently far from the stable equilibrium, then the trajectory is attracted by the coexisting stable CIC instead.In our example this happens with the second shock applied at t 500 as shown in Figure 5 b .The trajectory does not return to the equilibrium as could be expected but continues to fluctuate soft line as before the second shock.The trajectory continues to be attracted by the coexisting CIC because at this different time of application of the second shock for the new value of the parameter the variables have values corresponding to a point in phase space which still belongs to the basin of attraction of the coexisting CIC.This is the hysteresis effect.

Coexistence of Attracting Closed Invariant Curves
Let us now specify a fixed speed of adjustment of the goods market of region 2, α 2 0.1, and the extreme value of the degree of interaction of the two economies through trade, δ 1.In this case we are left with the two free parameters α 1 and β.For the specification δ 1 the equilibrium expressions 3.1 are simplified to

5.1
The stability region of the equilibrium in this case is shown in Figure 6.The part of the stability region in which the roots of the characteristic equation are all real is shown dark-shaded.The part in which some of the roots are complex conjugate is shown light-shaded.The upper right-hand "corner" is at β 1.0099, α 1 0.03663.We see that for the specifications adopted in the present case of the model, increase of the parameter α 1 , the speed of adjustment of the goods markets of region 1, destabilizes the equilibrium quickly for small values of α 1 .Equilibrium stability is possible only when α 1 is roughly less than half of α 2 .On the other hand, increase of the parameter β, the degree of capital mobility between the regions, destabilizes the equilibrium less quickly and does not appreciably affect the destabilization process due to α 1 .The Neimark-Hopf bifurcation curve can be identified as part of the boundary of the stability region and is shown as a continuous line.
As in the previous case, we have employed numerical simulations to obtain appropriate bifurcation and Lyapunov exponent diagrams.In Figure 7 we present the bifurcation diagrams of the variables Y 1 , Y 2 , and K 1 together with the corresponding Lyapunov exponent diagram, for β 0.1 with α 1 as the bifurcation parameter.These bifurcation diagrams provide visual estimates of the amplitude size of the occurring fluctuations.
In the present case of the model the Neimark-Hopf bifurcation curve gives rise to supercritical bifurcations.The CICs generated initially increase in amplitude in a continuous manner.But, for small values of the parameter β, such as the present value β 0.1, after a certain value of α 1 we meet with a discontinuous succession of attractors.This appears as a sudden jump observed clearly in the bifurcation diagrams where we have a sudden decrease of amplitude size.As seen in Figure 7 this sudden jump occurs at α 1 ∼ 0.40.
The sudden change of amplitude implies a discontinuous succession from one attractor existing before this succession value of α 1 to another attractor existing after this value.We found, however, that the succession value of the bifurcation parameter at which the discontinuity appears depends on the starting values used for the computation of the bifurcation diagram.While the discontinuity appears at α 1 ∼ 0.40 if the starting values are near the equilibrium, it appears at α 1 ∼ 0.45 if the starting values are far from the equilibrium.This is clear in Figure 8 showing in detail the results of the two calculations.It is the same phenomenon as in Section 4 Figures 2 and 3    This phenomenon has similar characteristics to the hypercrater bifurcation process illustrated in the previous section.But, instead of a "point-to-CIC" interplay inside the stability region of the equilibrium, we now have a "CIC-to-CIC" interplay outside the stability region.The second attracting CIC appears suddenly when the first one is still attracting in this example at α 1 ∼ 0.40 , just as in the crater bifurcation the attracting CIC appears when the equilibrium point is still attracting.The two CICs then coexist until the first CIC becomes unstable at α 1 ∼ 0.45 , and for larger values of the parameter the second CIC dominates the system remaining as the only attractor.
Implications of this coexistence of attracting CICs are similar to those of the crater bifurcation discussed in Section 4. Economic interpretation is that small shocks do not generally affect the fluctuations but larger shocks may lead to quite different fluctuations.Among the characteristics of both types of attractor coexistence is the associated hysteresis effect.
In Figure 11 we now illustrate this effect for the CIC-to-CIC interplay of the coexisting attracting CICs of the present example β 0.1 .Initially the value of the parameter is α 1 0.43.The first CIC is still stable but the second CIC already exists see Figure 8 .With starting values far from the equilibrium the trajectory settles on the first CIC, after a transient phase.At t 300 we apply a shock on the parameter, changing its value to α 1 0.45 for which the first CIC has become unstable.The trajectory begins to deviate from the first CIC and settles on the second CIC, the latter being the only attractor for the new value of the parameter.At t 600 a second shock on the parameter is applied, restoring it to its initial value.The hysteresis  and α 1 0.43 1st CIC stable , with two parametric shocks, shown as black dots, at times t 300 when the value of α 1 is changed to 0.45 1st CIC unstable and t 600 when α 1 is restored to its previous value.
effect is now revealed as the trajectory does not return to the first CIC but remains on the second CIC characterized by a different frequency.

Unfolding of an Attracting Closed Invariant Curve
We consider now parameter values as in Section 5, except for the value of the degree of interaction of the two economies through trade which is now fixed at δ 0.  The stability region of the equilibrium in this case is shown in Figure 12.As in the previous cases, the part of the stability region in which the roots of the characteristic equation are all real is shown dark-shaded, and the part in which some of the roots are complex conjugate is shown light-shaded.We see that for the specifications adopted in this case increase of the parameter α 1 , the speed of adjustment of the goods market of region 1, does not destabilize the equilibrium quickly, contrary to the case δ 1.
The Neimark-Hopf bifurcation curve can be identified as part of the boundary of the stability region as in the previous cases and is shown as a continuous line in Figure 12.The remaining part of the boundary is a flip bifurcation curve.The stability region is now much larger than in the previous case.Its upper right-hand "corner" is at β 1.1059, α 1 3.6587.
Again, numerical simulations are employed, and in Figure 13  with α 1 as the bifurcation parameter.The CICs occur through supercritical bifurcation and increase in amplitude size in a continuous manner.But, after a certain value of the bifurcation parameter α 1 we meet with the phenomenon of attractor unfolding as an intermediate phase between the initial phase of continuously increasing regular CICs and the subsequent phase of chaotic behaviour.The situation is depicted in Figure 14.Initially, the CIC develops for increasing values of α 1 after its generation at the supercritical bifurcation.The unfolding then occurs at α 1 ∼ 2.868.The evolution of the attracting CIC is shown in Y 1 , Y 2 and Y 1 , M 1 plane "projections."After unfolding, as α 1 is further increased the unfolded CIC becomes chaotic.

Summary-Conclusions
We have explored numerically a nonlinear 5D discrete interregional Kaldorian macro-dynamic model with fixed exchange rates and presented certain remarkable patterns of dynamical behaviour that occur.To the authors' knowledge these patterns of dynamical behaviour have not previously been presented for high-dimensional discrete systems.
In the first case we assumed equal speeds of adjustment of the goods markets of the regions α 1 α 2 α 1 and established that attracting and repelling CICs are generated inside the region of equilibrium stability.This is the first phase of a crater-type bifurcation process, the second phase of which is the subcritical Neimark-Hopf bifurcation occurring on the bifurcation curve.As a result there exists a region of parameter values of the model for which a stable equilibrium coexists with a stable CIC.The repelling CIC signposts the boundary of the basin of attraction of the equilibrium.Inside this parameter region of attractor coexistence, as well as after equilibrium loses stability by means of the subcritical bifurcation, the attracting CIC is far from the equilibrium.Economic interpretation is that small shocks do not affect the asymptotic behaviour of the system, but large shocks may lead to permanent fluctuations.As observed by Kind 11 , this is a more accurate description of the idea of corridor stability 17 , in that large shocks do not lead to totally unstable dynamics.The associated hysteresis effect was illustrated by numerical simulations.
In the second case we assumed a fixed speed of adjustment of the goods market of region 2, α 2 0.1, and an extreme fixed value of the degree of interaction of the two economies through trade, δ 1.A remarkable feature found in this case is the coexistence and dynamical interplay between two single attracting CICs, occurring outside the stability region.Again, the associated hysteresis effect was illustrated by numerical simulation.
In the third case we assumed a fixed speed of adjustment of the goods market of region 2, α 2 0.1, as before, but a smaller degree of interaction of the two economies through trade, δ 0.5.For the parameter constellation considered the CICs occur through supercritical bifurcation.A remarkable feature in this case is the eventual unfolding of the occurring stable CIC.This was illustrated graphically by means of successive 2D projection plots of the 5D attracting CIC for increasing values of the bifurcation parameter.
Since high-dimensional systems are largely unexplored, the present results are useful in demonstrating the feasibility of the numerical approach for detecting and describing patterns of complex dynamical behaviour, such as bifurcation processes and coexistence of attractors, generated by high-dimensional discrete systems of economic dynamics.

Figure 1 :
Figure 1: Region of stability of the equilibrium for α 1 α 2 α 1 in the δ, β parameter plane.Flip bifurcation curve is shown dashed.

1 hFigure 3 :
Figure 3: Bifurcation diagram of Y 2 and Lyapunov exponent for α 1, β 0.8, and various starting values of Y 2 , resulting from numerical experiments designed to detect the interval of attractor coexistence and the "invisible" repelling CIC of the subcritical bifurcation.

Figure 5 :
Figure 5: Trajectory of Y 2 for α 1 and δ 0.55 initially attracted by the stable equilibrium.After the first shock when the value of δ is changed to 0.6 the trajectory is attracted by the stable CIC. a After the second shock applied at t 515 the trajectory is attracted by the equilibrium.b After the second shock applied at t 500 the trajectory continues being attracted by the CIC.The two parametric shocks are marked by dots.
indicating coexistence of attractors, but now the coexisting attractors are both closed invariant curves.The two attractors coexisting in this example are shown in Figures 9 and 10 .

1 −0. 1 dFigure 8 : 1 Y 2 Figure 9 :Figure 10 :
Figure 8: a Details of the bifurcation diagram of K 2 for δ 1, β 0.1.a As computed with starting values far from the equilibrium.b As computed with starting values near equilibrium.

Figure 11 :
Figure 11:The hysteresis effect associated with coexistence of attracting CICs.Time path of Y 2 for β 0.1 and α 1 0.43 1st CIC stable , with two parametric shocks, shown as black dots, at times t 300 when the value of α 1 is changed to 0.45 1st CIC unstable and t 600 when α 1 is restored to its previous value.

Figure 14 :
Figure 14: Top row: development of the attracting closed invariant curve after the supercritical bifurcation.Second and third row: attractor unfolding.Fourth row: Eventually the closed invariant curve becomes a chaotic attractor.