Describing the Dynamics and Complexity of Matsumoto-Nonaka's Duopoly Model

Abstract

We study the duopoly model proposed by Matsumoto and Nonaka (2006), in which two firms produce two complementary goods, and there are externalities of different signs. We analyze the topological complexity of the model by computing its topological entropy with prescribed accuracy, and, in addition, we show when such topological complexity is physically observable by characterizing the attractors. Finally, we exploit the fact that reaction maps have negative Schwarzian derivative to show the existence of absolutely continuous (with respect to Lebesgue measure) ergodic measures, and, as an economic application, we compute the average profit for almost all initial conditions.

1. Introduction: The Model and Our Aims

In [1] a two-market model consisting in two firms, which produce differentiated goods is introduced. The first firm produces good 𝑥 in the first market, and the second one produces good 𝑦 in the second market. It is assumed that externalities of different signs exist. An externality occurs when the actions of firms directly affect the production capabilities of other firms other than through the price mechanism of the market. In this case, the positive externality comes from the market demand; that is, the sales possibilities of one firm are positively influenced by the production of the other firm. The negative externality is due to the cost functions of each firm, since the cost depends not only on its own production but also on the other firm’s production.

In this frame, we are ready to introduce the equations of the model. Although in [1] a more general model is introduced, we are interested in the following particular case. Inverse demand functions are given by 𝑝₁𝑥,𝑦=(𝛼−1)²−𝑥2+𝛼𝑦²,𝑝₂𝑦𝑥,𝑦=1−2+𝛽𝑥², (1) where 𝑝₁ and 𝑝₂ are the market prices of goods 𝑥 and 𝑦, respectively, and 𝛼∈[1,2]⊂ℝ and 𝛽∈[0,2]⊂ℝ. Each firm decides future production depending on the other firm’s choice and production externalities. Cost functions are given by 𝐶₁()𝐶𝑥,𝑦=2𝛼−1𝛼𝑥𝑦,₂𝑥,𝑦=2𝛽𝑥𝑦, (2) and thus, the profit function of each firm is given by Π₁𝑥,𝑦=𝑝₁𝑥,𝑦𝑥−𝐶₁=()𝑥,𝑦𝛼−1²−𝑥2+𝛼²𝑦²()Π𝑥−2𝛼−1𝛼𝑥𝑦,₂𝑥,𝑦=𝑝₂𝑥,𝑦𝑦−𝐶₂=𝑦𝑥,𝑦1−2+𝛽²𝑥²𝑦−2𝛽𝑥𝑦. (3) It is assumed that each firm tends to maximize its profit. In order to get it, the firms can choose their production levels which would affect the other firm. The first firm maximizes the profit with respect to 𝑥, and the same occurs for the second firm respect to 𝑦, that is, 𝜕Π₁𝜕𝑥𝑥,𝑦=0,𝜕Π₂𝜕𝑦𝑥,𝑦=0. (4)

Hence, by solving (4), and under naive future expectations, the reaction functions of the firms are given by 𝑓_𝛼𝑦=𝛼𝑦−𝛼+1²,𝑔_𝛽(𝑥)=𝛽𝑥−1². (5) This is a static situation, but we are interested in the dynamic interactions between the two firms along the time. The function 𝑅_𝛼,𝛽𝑥,𝑦=𝑓_𝛼𝑦,𝑔_𝛽(𝑥=)𝛼𝑦−𝛼+1²,𝛽𝑥−1² (6) is the reaction function for the outcome (𝑥,𝑦). The iterations are given by 𝑥_𝑡+1=𝑓_𝛼𝑦_𝑡=𝛼𝑦_𝑡−𝛼+1²,𝑦_𝑡+1=𝑔_𝛽𝑥_𝑡=𝛽𝑥_𝑡−1²; (7) that is, 𝑥_𝑡+1,𝑦_𝑡+1=𝑅_𝛼,𝛽𝑥_𝑡,𝑦_𝑡=𝑓_𝛼𝑦_𝑡,𝑔_𝛽𝑥_𝑡, (8) where 𝛼∈[1,2]⊂ℝ and 𝛽∈[0,2]⊂ℝ. The parametric region, [1,2]×[0,2]⊂ℝ², is not chosen in an arbitrary way. Observe that the reaction map 𝑅_𝛼,𝛽 is a map from [0,1]×[0,1] into [0,1]×[0,1]. Note that 𝑥_𝑡+2,𝑦_𝑡+2=𝑅²_𝛼,𝛽𝑥_𝑡,𝑦_𝑡=𝑓_𝛼𝑔_𝛽𝑥_𝑡,𝑔_𝛽𝑓_𝛼𝑦_𝑡, (9) or, in other words, 𝑅²_𝛼,𝛽=𝑓_𝛼∘𝑔_𝛽×𝑔_𝛽∘𝑓_𝛼, (10) and therefore, one might expect that the dynamics of the whole system can be derived from the one-dimensional maps 𝐹_𝛼,𝛽∶[0,1]→[0,1] and 𝐺_𝛼,𝛽∶[0,1]→[0,1] given by 𝐹_𝛼,𝛽(𝑥)=𝑓_𝛼∘𝑔_𝛽(𝑥)=1+𝛼𝛽𝑥𝛽𝑥−2²,𝐺_𝛼,𝛽(𝑥)=𝑔_𝛽∘𝑓_𝛼𝑦=𝛽𝛼𝑦−𝛼+1²−1². (11)

Our aims are, on one hand, to give an analytical proof of the existence of chaos in the model. This proof will be given by computations of topological entropy of the system with prescribed accuracy. As positive topological entropy will imply the existence of chaos in the sense of Li and Yorke, we will be able to give a proof of existence of chaos for a wide range of parameter values. On the other hand, we will be able to describe the nature of attractors of the system, and we will go beyond to the seminal work [1] proving the existence of absolutely continuous ergodic measures. The existence of such measures will allow us to obtain some consequences from both dynamic and economic points of view. Finally, we must emphasize that our approach is somewhat different from that made in [2]. Moreover, our scheme can be adapted to analyze the dynamics of models given by maps with the form of 𝑅_𝛼,𝛽.

The paper is organized as follows. In Section 2 we study the complexity of the model. For that, we compute the topological entropy using different algorithms depending on the number of pieces of monotonicity of the functions involved in the model. As an approach to the attractors, their number is studied in Section 3. On one hand, the fact that the Schwarzian derivative of the function is negative allows us to know the possibilities of the metric attractors that can appear. On the other hand, a real computation of the attractors allows us to realize that topological chaos need not be physically observed, with the Lyapunov exponents being the key for analyzing these phenomena. The Lyapunov exponent gives us the key to study that situation. Invariant measures and average profits are considered in Section 4. Finally, Section 5 is devoted to conclusions.

2. Computing Topological Entropy

The goal of this section is to compute the topological entropy of the system. Topological entropy was introduced in [3] by Adler, Konheim, and McAndrew. When 𝑓 is a continuous piecewise monotone function, the Misiurewicz-Szlenk's theorem, [4], gives us the following characterization of the topological entropy  ℎ(𝑓).

Theorem 1. Let 𝑓:[0,1]→[0,1] be a piecewise continuous monotone functions and 𝑐(𝑓) the number of pieces of monotonicity. Then, ℎ𝑓=lim_𝑛→∞1𝑛log𝑐𝑓^𝑛. (12)

Although our system is two-dimensional, the computation of topological entropy can be reduced to a one-dimensional problem. For that, we use the power formula [3], ℎ𝑅_𝛼,𝛽=12ℎ𝑅²_𝛼,𝛽, (13) and the commutativity formula [5], ℎ𝑓_𝛼∘𝑔_𝛽=ℎ𝑔_𝛽∘𝑓_𝛼, (14) to obtain ℎ𝑅_𝛼,𝛽=ℎ𝐹_𝛼,𝛽=ℎ𝐺_𝛼,𝛽. (15)

Next, we must recognize that Misiurewicz-Szlenk formula is not useful for computing topological entropy of our model, so we are forced to use numerical algorithms which are based on the number of monotone pieces of 𝐹_𝛼,𝛽 and 𝐺_𝛼,𝛽. Therefore, we will make a precise description of monotonicity regions for these maps.

2.1. Monotonicity Regions of 𝐹_𝛼,𝛽 and 𝐺_𝛼,𝛽

We study the shape of the composition map, where shape means in this case the number of pieces of monotonicity that each map has depending on the values of the parameters 𝛼 and 𝛽. This will be used in the next section to compute the topological entropy.

The map 𝐹_𝛼,𝛽∶[0,1]→[0,1], 𝐹_𝛼,𝛽(𝑥)=(𝑓_𝛼∘𝑔_𝛽)(𝑥), in terms of monotonicity, has the following structure in the parametric space (see Figure 1).(i)𝐹_𝛼,𝛽 is a monotone decreasing map if (𝛼,𝛽)∈𝐹₀, where 𝐹₀∶=𝛼,𝛽∈[1,2]×[0,2]:0≤𝛽≤1−11−𝛼; (16) see example for 𝛼=1.5 and 𝛽=0.3 in Figure 2(a). (ii)𝐹_𝛼,𝛽 is a unimodal map with a minimum for 𝑥=(1/𝛽)(1−1−1/𝛼) if (𝛼,𝛽)∈𝐹₁, where 𝐹₁∶=𝛼,𝛽∈[1,2]×[0,2]:1−11−𝛼<𝛽≤1; (17) see example for 𝛼=1.5 and 𝛽=0.8 in Figure 2(b). (iii)𝐹_𝛼,𝛽 is a bimodal map with a minimum for 𝑥=(1/𝛽)(1−1−1/𝛼) and a maximum for 𝑥=1/𝛽 if (𝛼,𝛽)∈𝐹₂, where 𝐹₂∶=𝛼,𝛽∈[1,2]×[0,2]:1<𝛽≤1+11−𝛼; (18) see example for 𝛼=1.5 and 𝛽=1.2 in Figure 2(c). (iv)𝐹_𝛼,𝛽 is a trimodal map with two minima (for 𝑥=(1/𝛽)(1−1−1/𝛼) and 𝑥=(1/𝛽)(1+1−1/𝛼)) and a maximum for 𝑥=1/𝛽 if (𝛼,𝛽)∈𝐹₃, where 𝐹₃∶=𝛼,𝛽∈[1,2]×[0,2]:1+11−𝛼<𝛽; (19) see example for 𝛼=1.5 and 𝛽=1.8 in Figure 2(d).

(a)

(b)

Figure 1 

Parametric regions for the composition maps 𝐹𝛼,𝛽 (a) and 𝐺𝛼,𝛽 (b) according to the number of pieces of monotonicity.

(a)

(b)

(c)

(d)

Figure 2 

Graphics of the map 𝐹𝛼,𝛽 for 𝛼=1.5 and 𝛽=0.3, 0.8, 1.2, and 1.8 for each graphic.

The map 𝐺_𝛼,𝛽:[0,1]→[0,1],𝐺_𝛼,𝛽(𝑥)=(𝑔_𝛽∘𝑓_𝛼)(𝑥), has the following structure in the parametric space (see Figure 1).(i)𝐺_𝛼,𝛽 is a unimodal map with a maximum (for 𝑥=1−1/𝛼) if (𝛼,𝛽)∈𝐺₁, where 𝐺₁∶=𝛼,𝛽∈[1,2]×[0,2]∶𝛽≤1; (20) see example for 𝛼=1.5 and 𝛽=0.6 in Figure 3(a). (ii)𝐺_𝛼,𝛽 is a bimodal map with a maximum for 𝑥=1−1/𝛼 and a minimum for 𝑥=1−1/𝛼+1/𝛼𝛽 if (𝛼,𝛽)∈𝐺₂, where 𝐺₂∈[]×[]1∶=𝛼,𝛽1,21,2:1≤𝛼≤1+𝛽; (21) see example for 𝛼=1.3 and 𝛽=1.5 in Figure 3(b). (iii)𝐺_𝛼,𝛽 is a trimodal map with a maximum for 𝑥=1−1/𝛼 and two minima in the points 𝑥=1−1/𝛼−1/𝛼𝛽 and 𝑥=1−1/𝛼+1/𝛼𝛽 if (𝛼,𝛽)∈𝐺₃, where 𝐺₃∈[]×[]1∶=𝛼,𝛽1,21,2:1+𝛽<𝛼; (22) see example for 𝛼=1.9 and 𝛽=1.9 in Figure 3(c).

(a)

(b)

(c)

Figure 3 

Graphics of the map 𝐺𝛼,𝛽 for (a) 𝛼=1.5 and =0.6, (b) 𝛼=1.3 and 𝛽=1.5, and (c) 𝛼=1.9 and 𝛽=1.9.

2.2. Practical Computation of Topological Entropy

Algorithms for computing the topological entropy for the special case of unimodal maps are given in [6, 7]. An algorithm for bimodal maps (with three pieces of monotonicity) is described in [8], and an algorithm for a particular class of piecewise continuous maps with four pieces of monotonicity is described in [9]. These three types of algorithms cover all the cases involved in the dynamical system that we are considering for the composition maps 𝐹_𝛼,𝛽 and 𝐺_𝛼,𝛽. The results obtained for the model can be observed in Figure 4.

(a)

(b)

Figure 4 

Topological entropy with accuracy 10−3 for each (𝛼,𝛽)∈[1,2]×[0,2], 3D view (a) and projection (b).

Figures 5(a), 5(b), 5(c), and 5(d) show the topological entropy with accuracy 10⁻⁴, where 𝛼 takes the values 1.25, 1.30, 1.50, and 1.85.

(a) Topological entropy for 𝛼=1.25 and accuracy 10−4

(b) Topological entropy for 𝛼=1.30 and accuracy 10−4

(c) Topological entropy for 𝛼=1.50 and accuracy 10−4

(d) Topological entropy for 𝛼=1.85 and accuracy 10−4

Figure 5 

Topological entropy for different fixed 𝛼 while 𝛽∈[0,2].

For 𝛽=1.1, and 𝛽=1.15 the topological entropy with accuracy 10⁻⁴ is given by Figures 6(a) and 6(b).

(a) Topological entropy for 𝛽=1.1 and accuracy 10−4

(b) Topological entropy for 𝛽=1.15 and accuracy 10−4

Figure 6 

Topological entropy for fixed 𝛽, 1.1 (a) and 1.15 (b) while 𝛼∈[1,2].

2.3. Existence of Chaotic Maps

Before stating the existence of chaos, we must recall what we understand for chaos. A continuous interval map 𝑓∶𝐼→𝐼 is chaotic in the sense of Li and Yorke (LY-chaotic) if there is an uncountable set 𝑆⊂𝑋 (called scrambled set of 𝑓) such that for any 𝑥,𝑦∈𝑆, 𝑥≠𝑦, it is held that liminf_𝑛→∞𝑑𝑓^𝑛(𝑥),𝑓^𝑛𝑦=0,limsup_𝑛→∞𝑑𝑓^𝑛(𝑥),𝑓^𝑛𝑦>0. (23)

Positive topological entropy implies chaos in the sense of Blanchard et al. [10], and therefore, the above computations of topological entropy are an analytical proof of the existence of chaos in the model. Even more, we can prove that positive topological entropy is equivalent to chaos in the sense of Li-Yorke. For that, we just need to notice that when the topological entropy is equal to zero, the map is not chaotic because the relative extreme points of the composition maps 𝑓_𝛼∘𝑔_𝛽 and 𝑔_𝛽∘𝑓_𝛼 are nonflat since the second derivative in these points is nonzero, (remember that a critical point is nonflat if some higher derivative is nonzero). Then, by [11, Theorem A], 𝐶^∞ maps with nonflat turning points have no wandering intervals (i.e., for a continuous interval map 𝑓, an interval 𝐽 is called wandering if 𝐽, 𝑓(𝐽), 𝑓(𝐽)²,… are disjoint and no point 𝑥∈𝐽 is asymptotically periodic), and by [12, Lemma 2.7], a map with zero topological entropy is chaotic if and only if it has wandering intervals. Hence, our model is not chaotic when the topological entropy is zero.

Therefore, we have two regions the nonchaotic one where every trajectory has the property that for any 𝜖>0 there is a periodic trajectory which is 𝜖 close to it [13]. In practice, one might check that a computer simulation shows the convergence to a periodic orbit when parameters are in the nonchaotic region. For the second region one could expect a similar situation, but as we will show in the next section, the topological chaos we have shown to exist, could not be observed. For instance, for 𝛼=1.8 and 𝛽=1.85, the topological entropy of 𝑔_𝛽∘𝑓_𝛼 is positive, but the time series of several orbits reveal a periodic motion (see Figure 7).

Figure 7 

We draw a sample of 100 points of three orbits with initial values 0.25, 0.5, and 1 when the parameter values are 𝛼=1.8 and 𝛽=1.85. The picture shows the convergence to a two-periodic orbit which maps 0.033582 to 0.000140067 and vice versa.

3. Chaos, Attractors, and Schwarzian Derivative

We start this section with some useful definitions. Let 𝑓∶𝐼→𝐼 be a map. A probabilistic measure 𝜇 is said to be invariant for 𝑓 if 𝜇(𝐴)=𝜇(𝑓⁻¹(𝐴)) for any Borel set 𝐴⊂𝐼. In addition, the measure is ergodic if the equality 𝑓(𝐴)=𝐴 implies that 𝜇(𝐴)=0 or 1. Denote by ℳ(𝑋,𝑓) and ℰ(𝑋,𝑓) the set of invariant and ergodic measures of 𝑓, respectively.

Given 𝑥∈𝐼, define its 𝜔-limit set 𝜔(𝑥,𝑓) as the set of limits points of its orbit. Recall that a metric attractor is a subset 𝐴⊂[0,1] such that 𝑓(𝐴)⊆𝐴, 𝑂(𝐴)={𝑥∶𝜔(𝑥,𝑓)⊂𝐴} has positive Lebesgue measure, and there is no proper subset 𝐴⫅^′_/𝑞𝐴 with the same properties. 𝑂(𝐴) is called the basin of the attractor.

By [14, Theorem 4.1], for a multimodal map 𝑓∶𝐼→𝐼 without wandering intervals, there are three possibilities for its metric attractors.(A1)A periodic orbit (recall that 𝑥 is periodic if 𝑓^𝑛(𝑥)=𝑥 for some 𝑛∈ℕ).(A2)A solenoidal attractor, which is basically a Cantor set in which the dynamics is quasiperiodic. More precisely, the dynamics on the attractor is conjugated to a minimal translation, in which each orbit is dense on the attractor. The dynamics of 𝑓 restricted to the attractor is simple; neither positive topological entropy nor Li-Yorke chaos can be obtained. Its dynamics is often known as quasiperiodic.(A3)A union of periodic intervals, 𝐽₁,…,𝐽_𝑘, such that 𝑓^𝑘(𝐽_𝑖)=𝐽_𝑖 and 𝑓^𝑘(𝐽_𝑖)≠𝐽_𝑗, 1≤𝑖<𝑗≤𝑘, and such that 𝑓^𝑘 is topologically mixing. Topologically mixing property implies the existence of dense orbits on each periodic interval (under the iteration of 𝑓^𝑘).

Moreover, if 𝑓 has an attractor of types (A2) and (A3), then it must contain the orbit of a turning point, and therefore its number is bounded by the turning points. In the case of our model, we can strengthen the above result by noticing that 𝐹_𝛼,𝛽 and 𝐺_𝛼,𝛽 have negative Schwarzian derivative. Namely, the Schwarzian derivative [15–17] at a point 𝑥 is given by 𝑆𝑓(𝑥)=𝑓^′′′(𝑥)𝑓^′(𝑥)−32𝑓^′′(𝑥)𝑓^′(𝑥)². (24)

Since 𝑆(𝑓∘𝑔)(𝑥)=𝑆𝑓(𝑔(𝑥))⋅(𝑔^′(𝑥))²+𝑆𝑔(𝑥) and the maps 𝑓_𝛼 and 𝑔_𝛽 have negative Schwarzian derivative the same occurs to the composition maps 𝐹_𝛼,𝛽 and 𝐺_𝛼,𝛽. Since attractors of type (A1) must also attract the orbit of a turning point when the map has negative Schwarzian derivative, we have that the number of attractors of 𝐹_𝛼,𝛽 and 𝐺_𝛼,𝛽 is bounded to at most 3. Even more, when 𝐹_𝛼,𝛽 and 𝐺_𝛼,𝛽 have three monotone pieces, then their value in two turning points agree, which implies that in fact, the number of attractors of 𝐹_𝛼,𝛽 and 𝐺_𝛼,𝛽 can be at most 2. The following result shows that the attractors of 𝐹_𝛼,𝛽 and 𝐺_𝛼,𝛽 are deeply connected.

Proposition 2. Let 𝑓,𝑔∶𝐼→𝐼 be 𝐶¹ maps such that 𝜔(𝑥,𝑓∘𝑔) is an attractor of 𝑓∘𝑔. Then 𝜔(𝑔(𝑥),𝑔∘𝑓)=𝑔(𝜔(𝑥,𝑓∘𝑔)) is an attractor of 𝑔∘𝑓.

Proof. As usual, denote by 𝜆 the one-dimensional Lebesgue measure. Let 𝐴⊂𝐼 with 𝜆(𝐴)>0 such that for any 𝑦∈𝐴 we have 𝜔(𝑦,𝑓∘𝑔)=𝜔(𝑥,𝑓∘𝑔). Since 𝑓 is 𝐶¹, we have that 𝜆(𝑓⁻¹(𝐴))>0. Now, let 𝑧∈𝑓⁻¹(𝐴), and fix 𝑦∈𝐴 such that 𝑦=𝑓(𝑧). For any 𝑛∈ℕ we have 𝑔∘𝑓^𝑛(𝑧)=𝑔𝑓∘𝑔^𝑛−1𝑓(𝑧)=𝑔𝑓∘𝑔^𝑛−1𝑦, (25) from which we conclude that 𝜔(𝑔(𝑥),𝑔∘𝑓)=𝑔(𝜔(𝑥,𝑓∘𝑔)) is an attractor of 𝑔∘𝑓.

The attractors of the map 𝑅_𝛼,𝛽 can be obtained by combining the attractors of 𝑓_𝛼∘𝑔_𝛽 and 𝑔_𝛽∘𝑓_𝛼. Basically, we must use the inclusion [18, Theorem 2] 𝜔𝑥,𝑦,𝑅_𝛼,𝛽⊆𝜔𝑥,𝑓_𝛼∘𝑔_𝛽∪𝜔𝑦,𝑔_𝛽∘𝑓_𝛼×𝜔𝑔_𝛽𝑦,𝑓_𝛼∘𝑔_𝛽∪𝜔𝑓_𝛼(𝑥),𝑔_𝛽∘𝑓_𝛼, (26) and, in view of Proposition 2, check whether the attractors of 𝑓_𝛼∘𝑔_𝛽 and 𝑔_𝛽∘𝑓_𝛼 we are considering are linked or not. In any case, it should be noticed that even when 𝑓_𝛼∘𝑔_𝛽 and 𝑔_𝛽∘𝑓_𝛼 have just one attractor, the map 𝑅_𝛼,𝛽 may have several of them depending on their distribution in the 𝑋𝑌 plane (see, e.g., [18, 19] and the examples of attractors in Figure 12).

Obviously, if 𝑓_𝛼∘𝑔_𝛽 and 𝑔_𝛽∘𝑓_𝛼 have 2 attractors, the number of different attractors of 𝑅_𝛼,𝛽 increases. So, still we have a practical work to do to decide whether 𝑓_𝛼∘𝑔_𝛽 and 𝑔_𝛽∘𝑓_𝛼 have one or two attractors. We fix the map 𝑓_𝛼∘𝑔_𝛽 and realize that if it has three or four monotone pieces, then the image of its turning points is in {0,1}. We take 1𝑛^𝑛−1_𝑖=0𝑓_𝛼∘𝑔_𝛽^𝑖(1)−1𝑛^𝑛−1_𝑖=0𝑓_𝛼∘𝑔_𝛽^𝑖(0) (27) and converge it to 𝑓_𝛼∘𝑔_𝛽𝑑𝑢−𝑓_𝛼∘𝑔_𝛽𝑑𝜈, (28) where 𝑢 and 𝜈 are mutually independent ergodic measures of 𝑓_𝛼∘𝑔_𝛽 (see, e.g., [20]). Since different ergodic measures must be supported in different attractors, when the expression (27) tends to zero when 𝑛 tends to infinity we have an evidence that just one attractor exists. Hence, we make simulations to compute (27) showing the results in Figure 8. Note that we have to be cautious when 0 or 1 is eventually periodic, because in such case its trajectory ranges all the attractors if the periodic orbit is an attractor, which need not happen (e.g., we can address that the well-known example of the map 𝑓(𝑥)=4𝑥(1−𝑥) has a dense orbit, and hence, its attractor is the whole interval, and the image of its turning point is 1, which is eventually the fixed point 0). Fortunately, the probability of finding such pathological examples is 0, but this is the case for 𝛽=2 and 𝛼>1.

(a) 𝜖=0.1

(b) 𝜖=0.01

Figure 8 

We compute (27) for samples of length 𝑛=10000, 1<𝛼, 𝛽≤2, and constant step size ℎ=0.01 for both 𝛼 and 𝛽. The parameter 𝛼 is represented in 𝑋-axis. We draw in red the values of (𝛼,𝛽) such that (27) is greater than a fixed value 𝜖. In (a) 𝜖=0.1, and for (b) 𝜖=0.01.

Taking into account the results obtained for generating Figure 8, we see that for 𝛼=1.3 there are values of 𝛽 around 1.5 for which it seems that two different attractors of 𝐺_𝛼,𝛽 may coexist. Figure 9 shows a refinement of Figure 8 for 𝛼=1.3 and 1.35. The bifurcation diagrams (see Figure 10) also suggest the existence of two different attractors for 𝛼=1.3 but not for 1.35.

(a)

(b)

Figure 9 

(a) We compute (27) for samples of length 𝑛=100000, 𝛼=1.3, 1<𝛽≤2, and constant step size ℎ=0.01 for 𝛽. We see that (27) is clearly positive when 𝛽 is in a neighborhood of 1.5. (b) We do the same when 𝛼=1.35, but now we cannot conclude clearly the existence of different attractors.

(a)

(b)

(c)

(d)

Figure 10 

For 𝛼=1.3, (a) and (b) show the bifurcation diagram for initial conditions 0 and 1, respectively. The bifurcation diagrams are constructed with generating orbits of length 100200 and representing just the last 200 points. (c) and (d) show the same when 𝛼=1.35.

As an example, we show different types of attractors and limit sets (but nonmetric attractors) obtained for the parameter values 𝛼=1.3 and 𝛽=1.56 (see Figure 11 and the explanations therein).

(a) 𝑥0=1/𝛽,𝑦0=0

(b) 𝑥0=1/𝛽,𝑦0=0.001

(c) 𝑥0=1/𝛽,𝑦0=1

(d) 𝑥0=1/𝛽,𝑦0=𝐺(1)

(e) 𝑥0=1/𝛽(1−1−1/𝛼),𝑦0=1

(f) 𝑥0=1/𝛽(1−1−1/𝛼),𝑦0=0.999

Figure 11 

For 𝛼=1.3 and 𝛽=1.56, we make 200000 iterations and draw the last 100000 in the 𝑋𝑌-axis. Except for (b) and (e), the initial conditions are taken to be the extremum of 𝐹𝛼,𝛽 and 𝐺𝛼,𝛽. (a) and (e) are not attractors because the initial conditions are synchronized; in (a) 𝑔𝛽(𝑥0)=𝑦0, and in (e) (𝑔𝛽∘𝑓𝛼∘𝑔𝛽(𝑥0))=𝑦0. These are 𝜔-limit sets from initial conditions in a set of zero two-dimensional Lebesgue measures. If we change the initial conditions slightly, we obtain the attractors (b) and (f). Finally, (c) and (d) are two different attractors with the same projection in the 𝑋- and 𝑌-axis, while for (b) and (f) the projections are different. Attractors (b) and (f) are possible because 𝐹𝛼,𝛽 and 𝐺𝛼,𝛽 have two different attractors consisting of two-periodic transitive intervals (of type (A3)).

(a)

(b)

Figure 12 

For 𝛼=1.3 and 𝛽=1.2 we make 200000 iterations and draw the last 100000 in the 𝑋𝑌 axis. We take initial conditions 𝑥0=1/𝛽 in both cases and 𝑦0=0.001 for (a) and 𝑦0=(𝑔𝛽∘𝑓𝛼)(0.001) for (b).

When 𝛼=1.35 the map 𝑔_𝛽∘𝑓_𝛼 seems to have only one attractor, and, for 𝛽=1.2, we see in Figure 12 that several of them may exist (see [18] for an analytical explanation of this fact). Figure 13 shows different types of attractors that can appear when we fix 𝛼=1.86.

(a) 𝛼=1.86, 𝛽=1.25, 𝑥0=1/𝛽, and 𝑦0=1−1/𝛼

(b) 𝛼=1.86, 𝛽=1.54, 𝑥0=1/𝛽, and 𝑦0=1−1/𝛼

(c) 𝛼=1.86, 𝛽=1.56, 𝑥0=1/𝛽, and 𝑦0=1−1/𝛼

(d) 𝛼=1.86, 𝛽=1.74, 𝑥0=1/𝛽, and 𝑦0=1−1/𝛼

Figure 13 

For 𝛼=1.86 we have drawn several different 𝜔-limits. We make 200000 iterations and draw the last 100000 in the 𝑋𝑌 axis. We can see that for these parameter values the function is chaotic, but it could happen that the chaos is not physically observable.

Once we have described the attractors of 𝑅_𝛼,𝛽, we go back to the idea of explaining when topological chaos cannot be physically observed. To that end, it is worth to notice that for a 𝐶¹ interval map 𝑓, a periodic point 𝑥∈𝐼 of period 𝑝 is an attractor when |∏^𝑝_𝑖=1𝑓^′(𝑓^𝑖(𝑥))|≤1, and therefore, its Lyapunov exponent (see, e.g., [14]) is 𝜆𝑓(𝑥)=limsup_𝑛→∞1𝑛^𝑛_𝑖=1log𝑓^′𝑓^𝑖(𝑥)≤0. (29) Since turning points characterize the attractors of 𝑓_𝛼∘𝑔_𝛽 and 𝑔_𝛽∘𝑓_𝛼, we conclude that attractors of type (A3) may exist if and only if 𝜆(𝑔_𝛽∘𝑓_𝛼)(𝑐)>0 for some turning point 𝑐 of 𝑔_𝛽∘𝑓_𝛼. The following graphic shows the maximum Lyapunov exponent for the parametric region. Figure 14 shows the maximum of {𝜆(𝑔_𝛽∘𝑓_𝛼)(𝑐_𝑖),𝑖=1,2}, where 𝑐₁ and 𝑐₂ are the turning points that have different forward images by 𝑔_𝛽∘𝑓_𝛼.

Figure 14 

We take orbits of length 𝑛=10000 of 𝑔𝛽∘𝑓𝛼 with initial conditions on the image of the turning points and estimate the Lyapunov exponent for (𝛼,𝛽)∈[1,2]×[0,2] with step size of 0.01. In red we draw those values of the parameter for which the maximum of Lyapunov exponents for the turning points is smaller than or equal to 0.

It is observed that there is a white band in the upper right side of the projection; this corresponds with the region where the fixed point is a local attractor. The topological entropy in these points is positive, but the chaos is not physically observable. Figure 14 explains the results obtained in Figure 7 because the parameter values 𝛼=1.8 and 𝛽=1.85 are clearly in a region with negative Lyapunov exponents.

4. On Invariant Measures and Average Profits

One of the main results in [1] states that for 𝛼=2 and 𝛽=2 the map 𝑅_𝛼,𝛽 has an invariant measure 𝜇 which is absolutely continuous with respect to the two-dimensional Lebesgue measure. Such invariant measure is supported on the whole square [0,1]², and therefore, for almost every initial condition, (𝑥,𝑦) (up to a set of zero Lebesgue measure) holds that lim_𝑛→∞1𝑛^𝑛−1_𝑖=0Π𝑅^𝑖_𝛼,𝛽𝑥,𝑦=_[0,1]2Π𝑑𝜇, (30) where Π=(Π₁,Π₂) is the profit function. Therefore, even when its dynamics is chaotic, with probability one, the average profit can be estimated, and it decided whether the chaotic regime is more profitable than profits at equilibrium points or not. The same estimations can be done also for other economic quantities like production or costs.

The aim of this section is to prove that we can state the same result for a wide range of parameter values of 𝛼 and 𝛽. Note that now we have to take into account several facts. (i)The map 𝑅_𝛼,𝛽 may have several attractors. Therefore, we can have several average values. (ii)When we are in the basis of a periodic attractor 𝐴={(𝑥₁,𝑦₁),…,(𝑥_𝑘,𝑦_𝑘)}, the equality lim_𝑛→∞1𝑛^𝑛−1_𝑖=0Π𝑅^𝑖_𝛼,𝛽=1𝑥,𝑦𝑘^𝑘_𝑖=1Π𝑥_𝑖,𝑦_𝑖 (31) trivially holds. Therefore, we must concentrate on attractors constructed from types (A2) and (A3). (iii)Since, computationally, attractors of type (A2) are not easily detected and their probability is small, we must concentrate our efforts on deciding when attractors of type (A3) support an invariant ergodic measure which is absolutely continuous with respect to the Lebesgue measure. Note that, by, for example, [21], the existence of dense orbits is not enough to guarantee the existence of absolutely continuous ergodic measures.

We start with the next result, which states that absolutely continuous ergodic measures of 𝑓_𝛼∘𝑔_𝛽 and 𝑔_𝛽∘𝑓_𝛼 are deeply connected.

Proposition 3. Let 𝑓,𝑔∶𝐼→𝐼 be 𝐶¹ maps, and let 𝜇∈ℰ(𝐼,𝑓∘𝑔) be absolutely continuous with respect to the Lebesgue measure. Then 𝑔𝜇 defined by 𝑔𝜇(𝐴)=𝜇(𝑔⁻¹(𝐴)) for every Borel set 𝐴 is an absolutely continuous ergodic measure of 𝑔∘𝑓.

Proof. The fact that 𝑔𝜇∈ℰ(𝐼,𝑔∘𝑓) can be seen in [22]. Now, let ℎ be a density function for 𝜇; that is, 𝜇(𝐴)=∫_𝐴ℎ𝑑𝜆 for any Borel set 𝐴. Since 𝑓𝜇(𝐴)=_{𝑓−1(𝐴)}ℎ𝑑𝜆=_𝐴ℎ∘𝑔⁻¹𝑑𝜆, (32) we conclude that ℎ∘𝑓 is the density function of 𝑓𝜇, and the proof ends.

Proposition 3 is a great help in detecting the absolutely continuous ergodic measures of 𝑅_𝛼,𝛽; we just consider one map 𝐹_𝛼,𝛽 or 𝐺_𝛼,𝛽 and decide whether they have absolutely continuous ergodic measures. The absolutely continuous ergodic measures of 𝑅_𝛼,𝛽 will be constructed by applying Proposition 3 and taking the product measure, which will be absolutely continuous but probably not ergodic. If we denote the product measure by 𝜇₁×𝜇₂, we take the ergodic decomposition 𝜇₁×𝜇₂=∑^𝑘_𝑖=1𝑝_𝑖𝜈_𝑖, where ∑^𝑘_𝑖=1𝑝_𝑖=1, and 𝑘 is the number of different attractors constructed by the attractors of 𝐹_𝛼,𝛽 and 𝐺_𝛼,𝛽 supporting the measures 𝜇₁ and 𝜇₂. The absolutely continuous ergodic measures are therefore 𝜈₁,…,𝜈_𝑘. In the case studied in [1] we find that 𝑘=1.

So, we fix a turning point of 𝐹_𝛼,𝛽 for instance and compute its Lyapunov exponent 𝜆(𝐹_𝛼,𝛽)(𝑐). If it is negative or zero, there is nothing to do. So, we assume that it is positive and let 𝐽₁,…,𝐽_𝑘 be periodic intervals such that ∪^𝑘_𝑖=1𝐽_𝑖 is an attractor of type (A3). Any ergodic measure 𝜇 supported on ∪^𝑘_𝑖=1𝐽_𝑖 must satisfy 𝜇(𝐽_𝑖)=1/𝑘, 𝑖=1,…,𝑘, and then we can consider 𝐹^𝑘_𝛼,𝛽, the invariant subinterval 𝐽_𝑖 for some 1≤𝑖≤𝑘, and the ergodic measure (for 𝐹^𝑘_𝛼,𝛽) 𝜈 given by the formula 𝜈(𝐴)=𝜇(𝐴∩𝐽_𝑖)⋅𝑘. Since 𝐹^𝑘_𝛼,𝛽 has negative Schwarzian derivative, by [23], we conclude that 𝜈, and therefore, 𝜇, can be chosen to be absolutely continuous with respect to the one-dimensional Lebesgue measure.

For instance, we consider the values 𝛼=1.3 and 𝛽=1.56. The possible attractors were studied in detail in Figure 11, and we obtained the following average profits for the initial values as we showed in Table 1.

Table 1 

We compute average profits for the initial conditions considered in Figure 11 and following the iteration processes we make there we check that both firms should produce starting from the initial conditions on the basis of the attractor (𝑓). Note that the profits in the fixed point (Cournot equilibrium) are 0.0140236 for firm 1 and 0.148848 for firm 2, and therefore, the only firm interested in stabilizing the Cournot point will be firm 2.

We compute the average profit for each firm. This result could be different if we change initial conditions when the number of attractors is different from one, as we studied in the previous section. Thus, for each firm we are going to present two different graphics in the square (𝛼,𝛽)∈[1,2]² with two different initial conditions, and we observe the differences between them. In the square (𝛼,𝛽)∈[1,2]×[0,1] there is only at most one turning point. Then we study the number of fixed Cournot points depending on (𝛼,𝛽)∈[1,2]×[0,2]. When there are more than a fixed point we determine in which fixed point the profit is maximum for each firm. Finally we compare the average profit with the profit in the fixed point (in which it is maximum) and compute the residual profit.

Figure 15 shows the average profit of both firms. For that, we use orbits of length 𝑛=10000 and, as initial conditions, (𝑥,𝑦), where 𝑥 is the minimum of 𝑓_𝛼∘𝑔_𝛽 when it exists and 0.4, otherwise, and 𝑦 is the maximum of 𝑔_𝛽∘𝑓_𝛼. We will denote this election of initial conditions as min-max conditions. When 𝛽 is greater than one there are more turning points which eventually may produce other average profits as it is shown in Figure 16. Here, we show the profit of firm 1, taking as initial conditions (𝑥,𝑦), where 𝑥 is the maximum of 𝑓_𝛼∘𝑔_𝛽 and 𝑦 is the maximum of 𝑔_𝛽∘𝑓_𝛼, which will be denoted as max-max conditions, and compare the situation with the initial conditions min-max previously computed. Figure 17 shows the same scenario for the second firm.

(a)

(b)

Figure 15 

Average profits of firm 1 (a) and firm 2 (b) for 𝑛=10000 and step size for parameters 𝛼 and 𝛽 eqal to 0.01. The initial conditions are 𝑥0 the maximum of 𝑓𝛼∘𝑔𝛽 and 𝑦0 the maximum of 𝑔𝛽∘𝑓𝛼 (type max-max).

(a) Average profit of firm 1 with min-max initial conditions and 𝛽≥1

(b) Average profit of firm 1 with max-max initial conditions and 𝛽≥1

(c) Projection of the average profit of firm 1 with min-max initial conditions

(d) Projection of the average profit of firm 1 with max-max initial conditions

Figure 16 

Average profit for firm 1 with different initial conditions. Note that in projections, the darker the color is, the lower the average profit is.

(a) Average profit of firm 2 with min-max initial conditions and 𝛽≥1

(b) Average profit of firm 2 with max-max initial conditions and 𝛽≥1

(c) Projection of the average profit of firm 2 with min-max initial conditions

(d) Projection of the average profit of firm 2 with max-max initial conditions

Figure 17 

Average profit for firm 2 with different initial conditions. Note that in projections, the darker the color is, the lower the average profit is.

We are interested in analyzing whether the destabilization of Cournot points, even in the chaotic regime, is a good business for both firms. To this end, we compare the average profits of generic orbits with the expected profit at Cournot points. Firstly, we have to compute the number of fixed points of the composition 𝑓_𝛼∘𝑔_𝛽, which can be seen in Figure 18. The parameter space can be divided into two regions, 𝐴₁ and 𝐴₂. In 𝐴₁ there is one fixed point while in 𝐴₂ the number of fixed points is three, which will be ordered 𝑥₁(𝛼,𝛽)≤𝑥₂(𝛼,𝛽)≤𝑥₃(𝛼,𝛽). As an exception, for (𝛼,𝛽)=(2,2), the number of fixed points are four. In [1, Theorems 1 and 2] the stability and bifurcation diagrams of stationary points is studied in detail.

Figure 18 

Number of fixed points of the composition 𝑓𝛼∘𝑔𝛽. In 𝐴1 there is one fixed point; in 𝐴2 there are three fixed points. When (𝛼,𝛽)=(2,2) there are four fixed points.

In region 𝐴₂ we compute the profit of the first firm for each fixed point, and we observe that it is maximum for the fixed point 𝑥₃(𝛼,𝛽) and minimum for 𝑥₁(𝛼,𝛽), obtaining the converse situation for the second firm. Figure 19 (resp., Figure 20) shows the profit of the first firm (resp., second firms) against the profit obtained by the second firm (resp., second firm) when the fixed point that maximizes the profit of firm 1 (resp., firm 2) is chosen.

(a)

(b)

Figure 19 

Average profits of firms 1 (a) and 2 (b) when we choose the fixed point which gives us the maximum profit of the first firm 1 among all the fixed points. Note that in projections, the darker the color is, the lower the profit is.

(a)

(b)

Figure 20 

Average profits of firms 1 (a) and 2 (b) when we choose the fixed point which gives us the maximum profit of the first firm 2 among all the fixed points. Note that in projections, the darker the color is, the lower the profit is.

The case 𝛼=1.3 was studied in Figure 11 and Table 1. We compute the profit of both firms in this case to compare the results see Figures 21 and 22.

(a)

(b)

Figure 21 

For 𝛼=1.3 and 𝛽∈(1,2], we show the average profits of firms 1 (blue) and 2 (red) with initial conditions 𝑥0=(1/𝛽)(1−1−1/𝛼), 𝑦0=0.999. We also show the profit at the Cournot points for firms 1 (dashed and blue) and 2 (dashed and red). The first figure (a) shows in the 𝑋 axis the parameter 𝛽∈[0,2], and the second one (b) with 𝛽∈[1,2] is a zoom of the first one.

(a)

(b)

Figure 22 

In this case 𝛼=1.3, but the initial conditions are 𝑥0=1/𝛽 and 𝑦0=0.001. The code is the same as in Figure 21. With these initial conditions there are points in which the profit of firm 1 (close to 𝛽=1.56) is greater than the profit reached in the Cournot point.

5. Conclusions

A duopoly model depending on two parameters is studied in detail. By computing its topological entropy we characterize the parameter values which admit the existence of chaotic behavior of trajectories. We also analyze when the above mentioned chaotic trajectories either can be observed or remain on a set of zero Lebesgue measures and therefore are not physically observed in computer simulations. For that, we study the metric attractors of the model concluding that several of them can coexist. The existence of absolutely continuous ergodic measures is proved, and we showed that, even when chaotic maps are considered, the average of economic functions, like profit and production and along the orbits can attain a finite number of values. The existence of such values allows us to make decisions on whether the equilibrium points deserve to be stabilized by choosing strategies that are distinct from naive expectations.

Acknowledgments

This paper has been partially supported by Grants MTM2011-23221 from Ministerio de Ciencia e Innovación (Spain) and 08667/PI/08 from Programa de Generación de Conocimiento Científico de Excelencia de la Fundación Séneca, Agencia de Ciencia y Tecnología de la Comunidad Autónoma de la Región de Murcia (II PCTRM 2007-10), and COST Action IS1104 “The EU in the new economic complex geography: models, tools and policy evaluation.”