Learning Dynamics in the Cobweb Model with Heterogeneous Producers

In this study we investigate the nonlinear dynamic s of the traditional cobweb model with two types of heterogeneous producers who are risk avers e and seek to learn the distribution of asset price s, in terms of the sample mean and variance of historical prices, using the Arithmetic Learning Processes (ALP) over different window lengths. We show that h eterogeneity has a double edged effect on the dynamics in the sense that heterogeneous learning c an stabilize an otherwise unstable dynamics in some cases and destablize an otherwise stable dynamics i n other cases as well. When the steady state become s unstable, the model displays complicated dynamics t hrough a variety of types of bifurcations.


INTRODUCTION
Consider the well-known cobweb model: e t t p aq b = + for the supply and p t = αq t +β for the demand, where q t and p t are quantities and prices, respectively, at period t, e t p is the price expected at time t based on the information at t-1 and a,b,β(>0) and α<0 are constants. It is well known that, under the naive expectation scheme e t t 1 p p − = , the price either converges to the optimal market equilibrium (when |α/a|<1) or explodes (when |α/a|>1). When the producers are homogeneous, it has been shown that non-linearities in the supply or demand curves may lead the cobweb model to exhibit both stable periodic and chaotic behavior (see for example, Artsein [1] , Jensen and Urban [19] , Chiarella [7] , Holmes and Manning [17] , Hommes [18] , Puu [21] and Day [13] ). These authors consider a variety of backward looking mechanisms for the formation of the expectations e with 0≤w≤1. With risk averse producers, the traditional linear cobweb model becomes nonlinear (e.g., Boussard and Gerard [3] , Burton [6] and Boussard [4] ). By assuming that the actual price p t is uncertain so that e t p has mean t p and variance t v , Boussard [4] shows that, under the simplest learning scheme tp p = and Among various learning schemes, the properties of recursive learning processes under homogeneous expectations have been studied extensively (e.g., Evans and Ramey [15] , Balasko and Royer [2] , Evans and Honkapohja [14] ). Under the assumption of bounded rationality, producers are somewhat uncertain about the dynamics of the economic system in which they are to play out their roles, so they need to engage in some learning scheme. Apart from Boussard [4] , a great deal of the literature on learning has been devoted to schemes for the mean, but very little to schemes for the variance. Chiarella and He [7,9] extend Boussard's framework in a way that takes account of the risk aversion of producers and allows them to learn about the distribution of the unknown price over the next period. By assuming that producers' subjective estimates of mean and variance follow the Arithmetic Learning Process (ALP):  [7,9] show that the resulting cobweb dynamics form a complicated nonlinear expectations feedback structure whose dimensionality depends upon the length of the window of past prices (the lag length) used to estimate the moments of the price distributions. It is found that an increase of the window length L can enlarge the parameter region (in terms of |α/a) of the local stability of the steady state. In addition, through a detailed normal form and bifurcation analysis for window length of 2, it is revealed that, at the crossover from local stability to instability, the dynamics exhibit resonance behavior which is indicative of quite complicated dynamical behavior and even chaos (for the model with constant elasticity supply and demand functions). Given that it is more realistic to assume that the producers are heterogeneous, instead of homogeneous, in forming their expectations, this study extends the study in Boussard [4] and Chiarella and He [8,9] and considers a simplest case of heterogeneous producers by assuming that there are two types of heterogeneous producers undertaking bounded rationality learning. The aim of this study is twofold. To incorporate risk aversion of the two types of heterogeneous producers into the traditional linear cobweb model and to analyze the effect of the heterogeneous arithmetic learning process (with different window length) on the dynamics of the model.
The study is organized as follows. A general cobweb model with heterogeneous producers is established in Section 2. The heterogeneous arithmetic learning processes is introduced and the existence of the steadystate is then discussed in Section 2. In Sections 3 and 4, the dynamical behavior of the heterogeneous model, including stability and bifurcation analysis, is analyzed. In a companion paper Chiarella et. al. [12] consider the alternative heterogeneous geometric decay learning scheme, which is also widely studied in the literature.

The cobweb model with heterogeneous producers:
This section establishes a cobweb model when producers are heterogeneous in their risk and expectation formulation of both the mean and variance. In the case of linear supply and demand functions, the model may be written as Eq. 1 and 2: Supply p a q b (i 1 2) , , where q t corresponds to the aggregate supply, q i,t and e i t p , are, respectively, the quantity and price expected of producer i at time t based on the information set at t-1 and p t is the price and a i , b i , β(>0) and α<0 are constants. Our approach to the formation of expectations will be somewhat different in that we assume that the actual price p t is uncertain so that the heterogeneous producers treat e i t p , as a random variable drawn from a normal distribution whose mean and variance they are seeking to learn. It would of course be preferable (and more in keeping with models of asset price dynamics in continuous time finance) to treat e i t p , as log-normally distributed.
However this would then move us out of the meanvariance framework so we leave an analysis of this approach to future research.
The market clearing price and the heterogeneous model: Let i t p , and i t v , be, respectively, subjective mean and variance of price e i t p , of producer i formed at time t based on the information set at t-1 and q t be quantity at time t. With constant absolute risk aversion A i , the marginal revenue certainty equivalent of producer i is With constant absolute risk aversion A i we assume the certainty equivalent of the receipt r = pq is Then maximisation of this function with respect to q t leads to the marginal revenue certainty We recall that this objective function is consistent with producers having the utility of receipts function and hence the supply for producer i is given by Eq. 5: Denote by n i the proportion of type i producers in general, the proportion n i is a function of time t that is, n i,t , which could be measured by a certain fitness function and discrete choice probability, as in Brock and Hommes [5] . Because of the complexity of the dynamics, we consider only the case with fixed propitiation and leave the changing proportions problem to future work. then it follows from (1), (5) and t i i t q n q , = ∑ that the market clearing price is determined by Eq. 6: The homogeneous cobweb model: As a special case of the heterogeneous model (6), assume that producers are homogeneous, that is, a i = a, b i = b, i t which has been considered in Chiarella and He [8] .
A cobweb model with two types of heterogeneous producers following ALP: In the following discussion, the simplest heterogeneous model when there are two types of producers is considered. Let n 1 = (1+w)/2 and n 2 = (1-w)/2 with -1≤w≤1. Then (6) has the form Eq. 8: The model (8) is incomplete unless producers' expectations are specified. In this study, the Arithmetic Learning Processes (ALP) is assumed. More precisely, assume that producers form their subjective estimates of the mean and variance from the sample mean and variance, by considering past market clearing prices over a window of length L i , that is Eq. 9: where, i = 1,2 and L 1 , L 2 are integers.
Existence and uniqueness of the steady state: Denote by p* the state steady. Then, under the ALP, it is found from (6) and the relation In particular, for the two-type-producer model (8), the steady-state p* is given by Eq. 10: The following sections study the dynamics of the two-type-producer model (8) when producers follow the arithmetic learning process by using different window lengths L i . Assume L 1 ≤L 2 and denote L = max {L 1 , L 2 } = L 2 . Because of the dependence of the subjective mean t p and variance t v on price lagged L periods, equation (8) is a difference equation of order L (see system (A.1) in Appendix A.1). The local stability of the unique steady state p t = p* is determined by the eigenvalues of the corresponding characteristic equation (equation (A.2) in Appendix A.1), which is difficult to analyze in general. We therefore focus first on the case when L 1 = L 2 = L in Section 3 and then some special cases when L 1 ≠ L 2 , such as, L 1 , L 2 = 1,2,3,4 in Section 4. Set: and Eq. 11: Then γ 1 = β 1 /L 1 >0 and γ 2 = β 2 /L 2 >0. From the following discussions, we can see that the local stability of the steady state depends on the parameters, including those from supply and demand functions a 1 , a 2 , a, the proportional difference of the two types of producers w and the window lengths L 1 and L 2 used by the heterogeneous producers. Our discussion here focuses on two different aspects. On the one hand, for a fixed window length combination of (L 1 , L 2 ), we consider how the demand parameter α and the proportional difference w of producers affect the local stability of the steady state and bifurcation. On the other hand, for a set of fixed parameters, we examine how these results on the local stability and bifurcation are affected by different combination of the window lengths. Regarding the first aspect, it is found that both the local stability region and bifurcation boundary are easy to construct geometrically by using the parameters β 1 and β 2 , instead of w and α. The one-one relation (11) between (w,α) and (β 1 ,β 2 ) makes it possible to transform the results between the different sets of parameters and in addition, to preserve the geometric relation of the local stability regions. Note that the determinant of the Jacobian of the transformation (11) does not change the sign, implying the preservation of the transformation. In the following discussion, because of the geometric advantage just discussed, the results are formulated in terms of (β 1 ,β 2 ), although some of the stability regions are plotted using (w,α) as well.

Dynamics of the homogeneous ALP:
In this section, we consider the case when both types of producer use the same window length, that is L 1 = L 2 = L, then a relatively complete result on the local stability region of the steady state, the types of bifurcation and stability boundary is obtained in Proposition 1 for general lag length L. Proposition 1. For the nonlinear system (8), assume producers follow ALP and L 1 = L 2 = L. Then: • the characteristic equation of the unique steady state p* is given by Eq. 12: • p* is locally asymptotically stable (LAS) if Eq. 13: • the stability boundary β 1 +β 2 = L defines a 1: (L+1) resonance bifurcation. When L = 2 and 3, for fixed a 1 = 0.8<a 2 = 1, the stability regions and resonance bifurcation boundaries are plotted in terms of parameters (α,w) in Fig 2. In order to see the bifurcation features, numerical simulations are conducted to analyse the dynamics of the nonlinear system (8) for the cases L = 2,3,4 and 10

Proof. See Appendix A.3:
The local stability region and the resonance bifurcation boundary are plotted in Fig. 1 in the (β 1 , β 2 ) parameter space for general lag length L. In particular, there are two special cases of interest. When a 1 = a 2 , the local stability condition and the bifurcation boundary are independent of w, as expected. When a 1 ≠ a 2 , the local stability region for α becomes (i) α∈(-La 1 , 0] for w = 1; and (ii) α∈(-La 2 , 0] for w = -1. In other word, the local stability depends more on the ratio a 1 /a 2 and less on the population distribution w. Assume a 1 ≠ a 2 for L = 2the bifurcation boundary β 1 +β 2 = 2 becomes (1+w)/a 1 +(1-w)/a 2 = -4/α, which defines the 1:3 resonance bifurcation boundary w as a function of α, namely Eq. 14: Note that α<0 and, for a 1 /a 2 <(>)1,W(α) is an increasing (decreasing) function of α and hence, as w increases, the local stability region for α becomes smaller (larger). For general lag length L, the analysis in Appendix A.3 leads to the following Corollary. Corollary 2. For the nonlinear system (8), assume producers follow ALP and L 1 = L 2 = L. Then, in terms of the parameters α and w, an increase in L stabilizes the otherwise unstable steady state.
The above theoretical analysis of the local stability and bifurcation is verified by numerical analysis on the nonlinear system (8). Consider a special case when a 1 = a 2 = 1. In this case, the steady state is LAS for α∈(-2, 0] for any w∈[-1,1] and α = -2 leads to a 1:3 resonance bifurcation. Bifurcation diagrams for the nonlinear system (8) are plotted in Fig. 3 in terms of the parameter α<0 with different initial values. It is found that, when the initial values are close to the steady state (within 1% interval of the steady state), the bifurcation value α o is close to the theoretical bifurcation value α o = -2, as indicated in the upper panel in Fig. 3. When the initial values are not close to the steady state (within 400% interval of the steady state), the bifurcation value α o ≈-1.9 moves away from the theoretical bifurcation value α o = -2, as indicated in the lower panel in Fig. 3. In both cases, the nonlinear system (8) displays a simple type of bifurcation, which is a 3cycle as indicated by the phase plot in Fig. 4 and the time series plot in Fig. 5, over a wide range of the parameter α.
In order to understand the nature of the resonance bifurcation, let L = 2, in which case the instability of the steady state leads to a 1:3 resonance bifurcation. Consider the case as indicated in the lower panel in Fig. 3 and let β = 11, a 1 = a 2 = 1,A = 0.005, w = 0,b 1 = b 2 = 0. For α = -1.9 near the bifurcation value α o , a phase plot (in the space of (x t-1 ,x t ) for different initial values is plotted in Fig. 4. In this case, a strong 1:3 periodic resonance bifurcation leads to two sets of period three cycles P(p 1 ,p 2 ,p 3 ) and S(s 1 ,s 2 ,s 3 ), having the following behavior. When the initial values are close to the steady state, the solutions converge to the steady state p* and both P and S are unstable. When the initial values are not close to the steady state p*, it becomes unstable and solutions converge to one of the two sets of the period three cycles, either P or S, depending on the initial values. The dynamics of the nonlinear system (8) are very similar to those found in Chiarella and He [9,10] .
For L 1 = L 2 = 3,4, instability of the steady state leads to 1:4 and 1:5 periodic resonance bifurcations, respectively and similar dynamics along the bifurcation boundary are also found. To illustrate the periodicity of different resonance bifurcation, time series for L = 2,5 and 10 are plotted in Fig. 5. Similar outcomes (not reported here) are also found when a 1 ≠ a 2 .

Dynamics of the Heterogeneous ALP:
We now consider the case where each producer uses a different window length and let L 1 <L 2 = L. Unlike the case when L 1 <L 2 = L, the local stability regions of the steady state and bifurcation boundaries for different combination of lag lengths have less clear-cut features and become very complicated and difficult to analyse in general. To be able to see how the window lengths of the heterogeneous producers affect the stability of the steady state and bifurcation, a combination of analytical analysis and numerical simulation is used in the following discussion. Analytical results for L = max {L 1 ,L 2 }≤4 are summarized in Proposition 3, followed by a comparison of the local stability regions for various lag lengths. Some numerical simulations of various types of bifurcation are employed to demonstrate the complex of the heterogeneous ALP. • For (L 1 , L 2 )= (1,2), the steady state x* is LAS for: In addition, a flip bifurcation occurs along the boundary β 1 = 1 and a Neimark-Hopf bifurcation occurs along the boundary β 2 = 2 with two eigenvalues λ 1,2 =e ±(2πθ)I , ρ ≡ 2 cos(2πθ)∈[-2, -1] • For (L 1 , L 2 ) = (1,3), the steady state x* is LAS for": In addition, a flip bifurcation occurs along the boundary β 1 +β 2 /3 = 1 • For (L 1 , L 2 ) = (1,4), the steady state x* is LAS for: β ,β ∈ ≡ β ,β ; ≤ β ,β ,β < , In addition, a flip bifurcation occurs along the boundary β 1 = 1 and a Neimark-Hopf bifurcation occur along the boundary ∆ = 0 • For (L 1 , L 2 ) = (2,3), the steady state x* is LAS for: In addition, along β 1 = 2, a Neimark-Hopf and flip bifurcation occurs with λ 1 = -1 and λ 2,3 = e ±2π/3i along β 2 = 3 a Neimark-Hopf bifurcation occurs with λ 1 ∈[-1,1] and λ 2,3 = e ±2πθ with ρ ≡ 2 cos (2πθ)∈[-1,0] • For (L 1 , L 2 ) = (2,4), the steady state x* is LAS for: In addition, a flip bifurcation occurs along the boundary β 2 = 4 • For (L 1 , L 2 ) = (3,4), the steady state x* is LAS for: In addition, along β 2 = 4, a Neimark-Hopf and flip bifurcation occurs  Fig. 6. In all these cases, there is no saddlenode bifurcation and the nature of the Neimark-Hopf bifurcation is characterized by the value of θ and therefore of ρ, as indicated by subsequent discussion.

Comparison of the local stability regions:
To see the effect of the various learning processes, the local stability regions for different (L 1 , L 2 ) are combined together in Fig. 7. Comparing the stability regions of the steady state for different combination of (L 1 , L 2 ) leads to the following observations. Let D L = D LL = {(β 1 , β 2 ): 0<β 1 , β 2 , β 1 +β 1 <L} then D L ⊂D L' for L<L'. Implying that an increase of the lag length enlarges the parameter region of the local stability of the steady state.
For L = 1,2,3, D LL ⊂D L,L+1 . In addition, the local stability region D L,L+1 is significantly enlarged compared with D LL .
Based on these observations, regarding to the local stability region, one may draw the following conclusions: • When L 2 ≠L 1 +1, an increase of window length (either L 1 or L 2 ) enlarges the parameter region of the local stability of the steady state in general (e.g., D 11 ⊂D 13 ⊂D 14 ⊂ 24 , D 22 ⊂ 24 ) • When L 1 = L 2 +1, an increase of L 1 to L 2 does not necessarily stabilise an otherwise unstable steady state for certain regions of the parameters (e.g., D 12 D 22 , D 23 D 33 , D 34 D 44 ). In other words, homogeneity of the lag length (L 1 and L 2 ) may not have a stabilising effect • When the difference of the different lag lengths is small, in particular, when L 2 -L 1 = 1 (e.g., (L 1 ,L 2 ) = (1,2), (2,3) and (3,4)), the stability regions are significantly enlarged, compared with the homogeneous case of L 1 = L 2 • As indicated by Fig. 8, an increase in both the lag length and the population proportion for type 2 producers enlarges the stability region of the parameter α in general. However, when L 2 = L 1 ± 1 (e.g., L 1 = 2, L 1 = 1,3 in Fig. 8(a) and L 1 = 7, L 1 = 6,8 in Fig. 8(b)), an increase of the population proportion of the type 2 producers does not necessarily enlarge the stability region for α. In those cases, there is an optimal value of w leading to the largest stability region in α Complex dynamics under the heterogeneous ALP: Proposition 3 indicates that heterogeneous ALP can lead to various types of bifurcation. The variety of types of bifurcation and complicity of the dynamics are demonstrated through the case (L 1 ,L 2 ) = (1,2) in the following discussion. For (L 1 ,L 2 ) = (1,2) the characteristic equation of the steady state is given by Γ(λ) ≡ λ 2 +(γ 1 +γ 2 )λ+ γ 2 = 0 where γ 1 = β 1 and γ 2 = β 2 /2. Based on the analysis in Appendix A.4 (i), along the boundary β 1 = 1, β 2 ∈[0,2], there is an eigenvalue λ = -1, implying that a flip bifurcation occurs along this boundary. Along the other boundary β 2 = 2, β∈[0,1], there occur two eigenvalues λ 1,2 = e ±2πθI , satisfying ρ ≡ λ 1 +λ 2 = 2cos(2πθ) = -(β 1 +β 2 /2), λ 1 λ 2 = β 2 /2 = 1 and hence, the Neimark-Hopf bifurcation boundary is defined by β 1 = -1-ρ and β 2 = 2. It follows from β 1 ∈[0,1] that ρ∈[-2,-1]. The types of Neimark-Hopf bifurcation are determined by the value of θ and hence of ρ. If θ = p/q is a rational fraction, then so-called p:q-periodic resonance occurs. If θ is an irrational number, then one obtains quasiperiodic orbits. Therefore, the types of Neimark-Hopf bifurcation along the boundary are determined by the values of ρ∈[-2,-1]. The corresponding values of ρ required for p:q resonances to occur can be found from the table in Sonis [22] . The local stability region D 12 is transformed from the parameter space (β 1 , β 2 ) in Fig. 6 (a) to the parameter space (α,w) in Fig. 9(a) with the corresponding flip and Neimark-Hopf boundaries indicated.
Along the Neimark-Hopf bifurcation boundary, the types of periodic resonance (when θ = p/q) and quasiperiodic resonance (when θ is irrational) are determined by ρ = 2cos(2πθ)∈[-1,-2]. Note that, by solving (11), (α,w) is related to (β 1 ,β 2 ) by α = -[a 1 β 1 +a 2 β 2 ] and  Table 1 sets up the corresponding parameter values of (w,α) which give different types of resonances (with (p,q) = (1,2),(1,3),(2,5), (3,5), (1,5), (4,5) and one quasiperiodic orbit (with 2 θ = ). The above local bifurcation analysis and the variety of types of bifurcation along the Neimark-Hopf boundary are confirmed by our numerical simulations of the nonlinear system (8) when the parameter values are selected as indicated by Table 1. Points D,B and C in Fig. 9(a) correspond to a 1:3 and 2:5 resonances and quasi-periodic closed orbit, respectively. For the initial values near the steady state, the corresponding time series for the parameter values indicated by D,B and C in Fig. 9(a) converge to those three time series plotted in the left panel in Fig. 10. Corresponding to point D and B, (p,q) = (1,3) and (2,5) or (3,5), respectively, the periodicity of the cycles of the time series are clearly identified by the time series (on the left panel) and phase plot (on the right panel) in Fig. 10. In fact, corresponding to point B, the phase plot indicates clearly two sets of period 5 cycles.
Corresponding to point C, 2, θ = solutions with initial values near the steady state converge to the quasiperiodic time series, the bottom one on the left panel.
The quasi-periodicity of the time series is identified by the closed orbit of the phase plot, the bottom one on the right panel in Fig. 10.
As a further support to our bifurcation analysis for other combinations of the lag lengths, the case (L 1 ,L 2 ) = (2,3) is analysed in Appendix A.5 and similar results on various types of bifurcation are also found.
For fixed α = -2.494 and (L 1 ,L 2 ) = (1,2) and (1,3), bifurcation diagrams in the parameter w are given in Fig. 11. For (L 1 ,L 2 ) = (1,3) the local stability region of the steady state of the nonlinear system (8) is given in Fig. 9(b). One can see that, for (L 1 ,L 2 ) = (1,3) and the fixed α, as w increases, instability of the steady state leads to a flip type of bifurcation for a wider range of parameter of w, indicated in the upper panel of Fig. 11. However, for (L 1 ,L 2 ) = (1,2)for the fixed α, as w decreases (from w = -0.5), instability of the steady state leads to more complicated and richer dynamics, indicated by the bifurcation diagram over the range of w∈(-1,-0.6) in the lower panel of Fig. 11.

CONCLUSION
The dynamics of heterogeneous learning has been studied recently in Chiarella and He [10] when the expectations (of the first moment) of the heterogeneous agents follow various weighted average learning processes, the so-called a L process. It is found that the dynamics of the system, including stability, instability and bifurcation, are affected differently by different recursive learning processes and the heterogeneity has double edged effect on the dynamics-heterogenous learning can stabilize an otherwise unstable dynamics in some cases and destabilize an otherwise stable dynamics in other cases as well. The findings in this study provide further evidence along this line. In addition, when heterogeneous agents learn both the first moment and second moment as well, the general feature outlined in Chiarella and He [10] has been extended further so that a relatively complete picture can be drawn for the learning dynamics of the heterogeneous ALP. The interplay between the relative slopes of the supply and demand curves, the risk aversion coefficients and proportions of the heterogeneous producers as well as the lag lengths used, can lead to quite complex price behavior.

A.1. Characteristic Equation of the Heterogeneous Beliefs Model
Let x 1,t = p t , x 2,t = p t-1 , x 3,t = p t-2 ,…,x L,t = p t-(L-1) , where L = max{L 1 ,L 2 }. Then, (8) is equivalent to the following L-dimensional difference system Eq. A1: At the steady state p*, , 1  Without loss generality, it is assumed that L 1 ≤L 2 , then L = L 2 . Evaluating the function f(x t ) at the steady state, we then have: The result then follows from the Lemma in Chiarella and He [10] .

A.3. Stability Region for ALP with the Same Lag Length
For L 1 = L 2 = L, the steady state is stable for 0≤β 1 +β 2 <L and the bifurcation boundary is given by β 1 +β 2 = L, which can be written, in terms of α and w (and a 1 ,a 2 as well), as follows:  The relation α = F(w) defines a nonlinear function of w. Note that: Hence the bifurcation boundary boundary is defined by α = -L for w∈[-1,1] if a 1 = a 2 . The boundary is an increasing (decreasing) concave function of w for a 1 <(>)a 1 . In addition, α = -a 2 L for w = -1 and α = -a 1 L for w = 1. Hence, for fixed a 1 ,a 2 , the α parameter region for the local stability of the state steady is enlarged as the lag length L increases. In other words, increase of window length can stabilise an otherwise unstable steady state. A.4. Proof of Proposition 3 The characteristic equation for the ALP is given by Eq. (A.2) for general lag lengths L 1 and L 2 .
The stability region D 23 is transformed from the parameter space (β 1 , β 2 ) in Fig. 6 (d), to the parameter space (α, w) in Fig. A.1 with the corresponding flip and Neimark-Hopf boundaries indicated. Table A