Skraba, Andrej with Miroljub Kljajic, Davorin Kofjac, Matevz Bren and Mico Mrkaic, "Cobweb Model in System Dynamics Form - Hyperincursive Perspective", 2006 July 23-2006 July 27

Cobweb Model in System Dynamics Form
Hyperincursive Perspective

Andrej Skraba, Miroljub Kljajié, Davorin Kofjaé,
Matevz Bren and Miéo Mrkai¢é
University of Maribor, Faculty of Organizational Sciences
Kidriéeva cesta 55a, SI-4000 Kranj, Slovenia
e-mail:{andrej.skraba, miroljub.kljajic,
davorin.kofjac, matevz.bren, mico.mrkaic}@fov.uni-mb.si

Abstract

The cobweb model of competitive market dynamics has been examined in the form
of system dynamics model. Separation of the structure elements and introduction
of anticipative hyperincursive algorithm was used for transformation of the classical
cobweb model to the accelerator based one. The cyclical response of the system
that depends on the demand~supply parameters and eigenvalues of the characteristic
equation has been numerically examined. The concept of parameter differentiation
and time response of the system is transformed to the periodicity concept where
periodicity is the main, driven property of the model. As such this is the key attribute
in complex discrete agent-based adaptive anticipatory models. The periodic conditions
of the model has been analytically determined by the application of z-transform. The
periodicity conditions of the initial map have been preserved in the nonlinear case. By
the application of the Lyapunov exponents several stability regions of the nonlinear
model were numerically determined.

Keywords: cobweb, hyperincursion, system dynamics, anticipative system, nonlinear system, pe-
riod, Farey tree, chaos

1 Introduction

Cobweb model presents the market demand-supply adjustment. It is typically viewed as
the model of the agricultural pricing mechanism. The story behind the model is briefly ex-
plained as the common agriculture market adjustment mechanism: ”The quantity offered
for sale this year depends on what was planted at the start of the growing season, which
in turn depends on the last year’s price. Consumers look at the current prices, though,
when deciding to buy. The cobweb model also assumes that the market is perfectly com-
petitive and that supply and demand are both linear schedules.” For clear and extensive
introduction to the topic one should see [1, 2]. The fact that the studied cobweb model is
in the field of discrete dynamics [3] is rather an advantage since the systems of difference
equations are often easier to grasp. For example in his enduring scholarly value work
on the studies of Dynamic Systems Luenberger [4] on the first place addresses difference
and later on differential equations. The model in question has all the characteristics of
classical System Dynamics (SD) models: equilibrium, competitiveness, human perception,
Demand~Supply
Market Adjustment

: Ge : on
4

Figure 1: The main elements of the proposed cobweb SD structure

delay and adjustment but somehow it is avoiding to be settled in the common SD model
bank of each SD modeler. The main reason for elusiveness of the cobweb model is in it’s
original form which is not suitable for the straight transformation to the common elements
such as Level L and Rate R. The functions of demand Qa(k) and supply Q,(k) can be
specified in the form:

Qalk) = a+bP(k) (1)
Qs(k) c+dP(k-1) (2)

where a,b,c and d are parameters specific to individual markets. P(k) and Q,(k) should
be restricted to the positive values. In the cobweb model it is assumed that in any one
time period producers supply a given amount (determined by the previous time period’s
price) and then price adjusts so that all the product supplied are bought by customers. If
we write this in the form of equation then Qu(k) = Qs(k) which enables us to state that
the price is:

ne
; (3)

P(k) = SP(k—1) +

Equations 1, 2 and 3 are not quite in the proper form in order to perform the
transformation to the SD model. One of the things is the time argument k — 1. The
other is the missing Rate R elements and corresponding At. One should expect that the
transformation will provide the known equations in the familiar form for structure shown
in Fig. 1. The developed model should enable us to examine the properties of the cobweb
model and also to consider it’s structural and incursive perspective.

As the Wiener’s cybernetics [19] principle stands firm in the systems theory the cob-
web model principle stands as the basic linearized principle construct for the systems
interaction dependance and will probably be further more the basic starting tool for the
quantitative analysis of complex systems.

2 Transformation to SD form

Fig. 1 shows, that the price P and quantity Q of the goods should be stated as the
Level elements depending on the Rates which determine the change in price and quantity
supplied. The theory behind the cobweb model states that quantity supplied in present
depends on the ruling price in the past. Therefore the price and the quantity supplied
should be dependant variables as illustrated in Fig. 1. Restating the above Equations
eliminating the time argument k — 1 gives us the following set of equations:
Qalk+1) = a+bP(k+1) (4)
Q(k+1) = c+dP(k) (5)

Pik+1) = SPqy+ 5" (6)

The Eq. 4, 5 and 6 will enable the determination of the Rates elements. Let us
determine the Rate element for the change of Price Rp(k). As the equations are in the
difference form the Rate will be determined as R(k) = L(k +1) — L(k):

dP(k) — Qs(k
Rp(h) = P(k +1) — PR) = SFP) — Go) 7)
Little bit more work will be needed for the Rg(k) since special time consideration had
to be taken. We will apply the time arguments of k +1 and k + 2 in order to loose the
argument k — 1 which is present in the Eq. 2:

e+dP(k*) — Qs(k*
Rak) = Qa(h-+ 2) — Qa(k-+1) = aot PE) — Ca) (8)
Since the time argument k* with consideration of the Eq. 1 and 2 actually represents
the past i.e. argument of k — 1 we should state the equations for P(k — 1) and Qs(k — 1).
The Eq. 3 will enable us to state P(k — 1):

p(r—1) = PW aera (9)

Qs(k—-1)=a+P(k—1)b (10)
The Eq. 10 is set by the fact that Qu(k) = Qs(k) and Eq. 4. The consideration of
the time argument k—1 is necessary in order to perform calculation in the model. At each

time step the previous values are needed in order to perform the calculation. By putting
the P(k—-1)= bPth)-cta in the Eq. 10 we get:

b?P(k) — be + ab
d

By putting the Eq. 11 and 9 into Eq. 8 we get the simplified form for the rate
equation:

Q.(k-1) =a4+ (11)

Ro(k) = -a+¢-(b—d)P(k) (12)

As the result of the above algebraic manipulation the cobweb model could be stated
in the standard SD form:

P(k+1) = P(k)+AtRp(k) (13)
_ c+ dP(k) — Qs(k)

Rp(k) = (14)

Q(k+1) = Q(k)+AtRo(k) (15)

Ro,(k) = —a+c-(b-d)P(k) (16)

3
yo
yo
yO
te
Yo

[|

Figure 2: System Dynamics structure of the cobweb model

with the starting conditions P(0) = 45% and Q,(0) = 2 where x represents the starting
perturbation of the model. In the above set of equations the At is introduced which is not
present in the classical cobweb model formulation. If At = 1 then the model is equivalent
to the classical cobweb.

Fig. 2 shows the SD structure of the cobweb model corresponding to the Eqs. 13, 14,
15 and 16. There are two levels represented, P and Q, and two rate elements Rp and
Rg,. The model behavior is determined by the input parameters a, b,c and d as well as the
perturbation parameter p. The element Po represents the initial value of the level element
P. The initial value of the level element Q, is equal to the arbitrary value of perturbation
p.

The response of classical cobweb model developed by SD methodology is shown in Fig. 3
and Fig. 4. The applied parameter values with the description of the system response are
shown in Table 1.

a) Stable b) Unstable
20 > 60
ak pee ee oe dn ee
16]? : 20 ert
14 y 6 {|
Po t | Deal P |
28,
10 t 20 Sip
é = a
Bogs
6

-60
0 50 100 150 200 250 -300 -100 0 100 300 500 600
Q Qa

Figure 3: Response of the SD cobweb model: a) Stable, b) Unstable

There are three possibilities: a) Stable system, where the supply and demand converge,
b) Unstable system where the supply and demand diverge and c) Dynamically stable shown
in Fig. 4, where price and demand neither converge nor diverge.

Dynamically stable response indicates the periodical solution which will be of interest
in further examination of the model. In general a solution yp is periodic if Yim = Yn for
some fixed integer m and all n. The smallest integer for m is called period of the solution.
In our case the solution in Fig. 4 is a two-cycle solution.
In general the following definition will be applied [1]:

Definition 2.1 If a sequence {y} has e.g. two repeating values y; and yo, then y, and
yo are called period points, and set {y1,y2} is called a periodic orbit.

Periodical response of the system is important because real agricultural systems depend
on the cyclic behavior and could be controlled only by regarding the periodicity of such
systems. Examples from real cases could easily be found in crops as well as in the stock.

Dynamically Stable

é
op
60 wee os

0 25 50 75 100 125 150

Figure 4: Response of the SD cobweb model: Dynamically stable

2.1 Separation of the structure elements

The structure of the model in Fig. 2 shows that the Price and Quantity are related.
However the structure can be represented in a different way. By transforming the cobweb
model to SD form the model could become non-autonomous depending on the variable
At. The following two equations represent the different formulation of the cobweb model:

Qak+1) = c+ a@—@ a7)
P(k+1) = _— as)

This reformulation represents Q, and P as the non-related quantities. The only bound
that exists are the coefficients. In order to formulate the complete SD model the rate
elements which determine initial conditions shonld he determined:

Table 1: Parameter values

a b c d p Response

400 —20 -50 15 20 Stable
200 —8.1 -—43 12 90 Unstable
160 -—2 -20 2 55 Dyn. Stable

Rp(k) = P(k+2)— P(k+1) = he ae _

7 (P+ 1) — P(k)) (19)

Ra, (k) = Qs(k + 2) — Qs(k +1)

crate, 4Ql—e_

= FQslk +1) ~ Qath)) (20)

In order to meet initial conditions of the model the Q,(k — 1) should be determined:

Qelk=1) = a+ (Qk) —0) (21)

Equations for P and Q, in standard SD form are the following:

P(k+1) = P(k) + AtRp(k) (22)
Rok) = $(Py- PO e+") (23)
Qulk+1) = Qslh) + AtRo,(h) (24)
Ro,(k) = §(Q4(k) - (a+ 4(Q.(8) -9)) (25)

Eqs. 22, 23, 24 and 25 represents the cobweb model in SD form separated as shown in the
Fig. 5. Note that the terms for P and Q, are related only to the coefficients a,b, c,d and
p. P(k +1) is dependent only on the value of P(k) and coefficients a,b,c, d and p, not on
the Qs. Respectively for the Q,(k + 1).

abcdp
naan

8 P LEP), +

Re Ras

1 1

Figure 5: Cobweb model in SD form ~ separated elements

2.2 Anticipative formulation

Comparison of the structures shown in Fig. 2 and Fig. 5 indicates, that P and Q, depend
only on the parameter values a,b, c,d and p i.e. the initial conditions. Eqs. 22, 23, 24 and
25 enable the determination of entire anticipative (future event) chain while equation:

6
Plk-l)= bP(k) —c+a

5 (25

and Eq. 21 enable the determination of feedback (past event) chain. The representation of
the Feedback ~ Anticipative chain is shown in Fig. 6. The dynamics of interest is therefore
the chains dynamics which is dependant on the parameters a(t), b(t), c(t), d(¢) and p(t).
Both chains are actually dependant on strategy dynamics which could be formulated as
the f(a,b,c,d,p,t).

Feedback chain Anticipative chain

intr
[Ptken) James [ Pea) JL Peet) J] Po JL Pte t) J] Pema) ]---->[ Pein) ]

[rsthen) ]——--s- | Q5(6-2) JL Qsthet) J] 50) [stk] [ a.(k+2)]---- = [astern]
i
Initial Conditions
a(t) bi) c(t) a(t) pit)

Figure 6: Feedback ~ Anticipative chain

Application of hyperincursive algorithm and inspection of gained equations with Dubois’ [6]
formulation of logistic growth and previous research [11, 12] yields the following set of
equations for the hyperincursive cobweb model:

P(k+2) = 5(4-(F=**)) (27)
Q(k+2) = s(c-«-5(P-<)) (28)
with initial conditions:
Po(k+1) = ee (29)
Pk) = pag Ara (30)
Qso(k+1) = p (31)
Qu(h) = 2+ 5(Qulk-+1)-c) (32)

Terms A(k) and B(k) in Eq. 27 could be replaced by the terms P(k+1) or P(k), similarly
C(k) and D(k) in Eq. 28 by Qs(k +1) or Qs(k). This yields 16 different combinations of
system defined by Eq. 27 and Eq. 28 that should be studied. The appropriate explana-
tion of modified system structure should follow the techniques of graphical solutions for
homogenous difference equations [22].

The system combination further examined will have the following terms: A = P(k +
1), B = P(k),C = Qs(k +1) and D = Q,(k). This yields the following set of equations:
P(k +2)

(p+ - (Oey) aa)

§(Qs(k +1) -a—2(Qu(k) ~0)) (34)

Eq. 33 and Eq. 34 with initial conditions stated by Eqs. 29 ~ 32 could be modelled
as shown by the structure in the Fig. 7. The structure represents the cobweb model in
hyperincursive form modelled by classic SD elements. The Euler integration method is
applied with time-step At = 1.

Qs(k + 2)

abcdop
TTTTT

i Ti nt

y y

(+1) (ky () Ket}
Rp Rp Rp Ras Ras Ras
(k+2)| (+1) « 0) (k+1)| (k+2)}

Figure 7: Structure of the hyperincursive Cobweb model; Euler integration, At = 1

Eq. 33 and Eq. 34 could be reformulated in order to show the dependency of the
future-present-past events:

P(k) = so a +S P(b+1) (35)
Qalk) = FQsle +1) + 5Qs(k- 1) +a-% (36)

Eq. 35 and Eq. 36 state that the value of the present is dependent on the past as
well as on the future. This paradoxical statement is realizable since the formulation of
feedback anticipative chain could be stated. Fig. 7 has two delay chains, one for P and
one for Q,. One might notice, that the level and rate elements are dependant only on the
coefficients and initialization values.

3 Results

Fig. 8 represents some examples of hyperincursive cobweb model. The matrix of graphs in
Fig. 8 represents different modes of cyclic behavior named by the shape of the 2d mapping
of P(k +1), P(k + 2) (second column represents the shapes in 2d mapping). Each row
is named by second column 2d mapping shape and represents the response of the model
in four different mappings in each column different mapping. First column represents the
dynamics in the classical cobweb form. 3’¢ and 4‘* columns represents the 3d mapping
which is 2d mapping expanded with the time step k on x-axis. Fig. 9 represents cobweb and
2d mapping of the system before, in and after pentagon synchronization. Here pentagon
shape is taken only as the example. Each of the presented shapes in Fig. 8 could be
observed in the similar way. First pair of pictures in first row represents the behavior
of the system before the synchronization. Notice the arrows in the first row under 2d
mapping which show the movement of the vertices of the graph. The vertices converge to
the one point which is shown in the second row. Points cover one another and in the case
of pentagon synchronization estimated values of synchronization parameters are: a = 400,
b = —20, ¢ = —50, d = —12.3671 and p = 160. Third row in Fig. 9 represents the state
of the system after the pentagon synchronization where arrows represent the divergent
movement of the vertices points of the graph. The system is in transition to the next
full synchronization which is estimated as the quad synchronization where parameter d
is near 0. The synchronization representation in 3d might be better observed than in
the 2d mapping. The responses of the system in Fig. 8 were gained according to the
parameter values gathered in Table 2. The change was made in parameter d which yielded
the synchronization patterns as shown by the shape column. The parameter values were
gained by the simulation where the range of parameter d was set as of [—40, 40] with
Ad = 0.001. The condition for parameter values determination was set by the rule of
acceptable error between simulation steps and definition 2.1 of synchronization .

4 Periodicity of the system

The z-transform is the basis of an effective method for solution of linear constant-coefficient.
difference equations. It essentially automates the process of determining the coefficients
of the various geometric sequences that comprise a solution [4]. The application of z-
transform on the Eq. 33 and Eq. 34 with initial conditions stated by Eqs. 29 ~ 32
gives:

rps =z + yodz — yor?
re) -l+dz-2 (37)
Inverse z-transform yields the following solution:
n (a ~V-4+ é)"
Y-'2) = 2% y (d— V4 BY" —
@) w ( ) - aa
a" yd (d- VA)" »
$y (dt V4 4+ 2)
V4t@ w( )
an (d+ V-aF®)" aI yod (a +V-4F@)"
+ = (38)

Vase vate

9
Cobweb 2d mapping 3d mapping, P, Qs 3d mapping, P

200 ———- | _Osike2) P(k+2)
ee 10
100 5

Triangular
Qs (k+2)
o 8
P(k+2)_,
oa 8

°
oa
3
°
oa
3

i

Quad
Qs (k+2)
°
Ro
N

2 oO 2 oO
PUK#2) x 40° Ptk#t) x 10°

Pentagon
4 Qs(k#2)
oe 08

q
eo 8s

Pentagram
Qs (k+2)
= 8

o BES
oud

5 0 51015 5 0 5
P(k+2) Plk+1)

3
a

z
£500) 7
6 0

Hexagon
Bes
888

eo 88

ag

0 20 40 60 20 40 60 P(k+2)” O k Pike) Ok
P(k+2) P(k+1)

°

000

Nonagram

4 Qs (k#2) _,

3

oe 8

=

8 ose
8
°
8
&
8

-20 0 20 40 60
P(k+2) P(kH1)

Hexagram
Qs (k+2)
3

eo 8

a
m=
PkpZ

onda

0 5 10 15 0 5 10 15
P(k+2) Plk+1) P(k+2) Pik+1)

Figure 8: Results of the hyperincursive cobweb model; examples of cyclic behavior; notice
rows and columns naming

1
Cobweb 2d mapping

Qfk+2)
2

Before synchronization

20 0 2 40 20 0 20 40
P(k#2) P(k+1)

‘Synchronization

Q{k+2)
°

After synchronization

20 o 20 40
P(k+2)

Figure 9: Cobweb and 2d mapping of the system before, in and after pentagon synchro-
nization

In order to gain conditions for the periodic response of the system the following equa-
tion should he solved:

Y(z) = 9 (39)

Let us compute a numerical example of periodic solution applying the z-transform. The
period examined will be the period of 9 i.e. n = 9. In Eq. 39 one should put the condition
n=9. One of the possible solutions for the initial condition worth of examination is the
following:

= ——— +501 +iv3))? (40)

i
(4(-1 + iv3))*
The term (—1+ iV3)3 (let us denote the term as 2*) could be expressed by three different
imaginary values in polar form:

2 2
Z = ¥2( cos 2 +isin =) (41)
B= ¥2( cos a +isin =) (42)
3 14 14
3 = ¥2( cos >= +isin +) (43)

am
By putting Eq. 41, Eq. 42 and Eq. 43 into Eq. 40 and performing trigonometric
reduction one gets the following solutions:

Qn 4a 80
dy = 2cos dy = 2cos dg = 2cos (44)

By inspecting the Eq. 40 and considering the equation for the roots of complex
numbers [15]:

2k
Yz= vr( wae 7 + isin sao (45)

the general form of the solution for the parameter d could be assumed:
d = 2cos Pua (46)

where n is the period and m = 1,2,3,...,n — 1. Similar procedure could be performed for
the arbitrary period n. More general solutions which regards also the parameter b which
was fixed for the purpose of determination of solutions is:

d = 2bcos fern (47)

The solutions could in some cases be expressed in alternative algebraic or trigonometric
form. Tab. 2 shows some of the solutions for the parameter d. Alternative solutions could
be expressed as the roots of the polynomial. The table incorporates the Shape symbols
which are important at the study of the response of dynamical systems. Especially this is
the case at the examination of complex nonlinear dynamical systems [16] or their represen-
tation in SD form [29]. Mappings of the system and visualization of the periodic solution
is therefore important for the analyzing of periodic or chaotic solutions of differential and
discrete difference equations. Shape description is taken from the vocabulary of proper
shapes although the response of the system is primarily not proper. Numerical values
of the solutions for parameter d are important since this values also confirm the findings
of Sonis [24, 5] about the domain of attraction for 2D dynamics by n-dimensional linear
bifurcation analysis. One of the important conditions gained by the proposed inspection
is the value of the period n = 10 which is in close relation to the period n = 5. The value
of parameter dis d= 3(1 + 75) with numerical value 1.61803... This solution represents
the ” Golden Ratio” (#). Some of the different representations of solution for parameter d
value at period n = 10 are:

1
dip = = 2eos = = 5 (1+ V5) = 1.61803308874989481820... (48)

The first solution of parameter d at period n = 10 connects the considered discrete system
with the Fibonacci numbers given by the infinite series:

d al ( yer 49
gp ate (49)

u=t

The fact, that the periodicity conditions of the examined discrete system incorporates
the golden ratio number ¢ could be observed in other studies [8] of complex nonlinear

12
expansions of the basic cob-web systems e.g. Brock and Hommes ” Almost Homoclinic
Tangency Lemma”. One should expect that the symmetric response in n — mapping
should follow the pattern with the match in certain point of solution with the @ condition.
The source of the mentioned condition is presented by the preceding procedure. (The
value of parameter d for mentioned period n = 5 is d= BES = 0.61803... often called the
” Golden Mean”.) (*)Since the period 2 is on the boundary of the solution the periodic
response of the system depends on the initial conditions. Example of numerical values of
the period-2 response: a = 1,b = 1,c = -1,d = —2,p = 1. (**) The value for tetragon
is taken in the limit since the system of equation returns the undefined value when d = 0
therefore one should consider the tetragon period condition as the value approaching to
zero i.e. d — 0. In this case the system response is undetermined (2) in it’s critical point.

Table 2: Synchronization parameter values

Shape Description a b c d p
4 Triangular 400 —20 —50 20.0000 160
° Quad 400 —20 —50 -—0.0010 160
Q Pentagon 400 —20 —50 —12.3671 160
* Pentagram 400 -—20 -—50 32.3620 160
OQ Hexagon 400 -—20 -—50 -—20.0000 160
OQ  Nonagram 400 —20 -50 -6.9450 160
.° 3 Hexagram 400 —20 —50 15.3070 160

Discrete 2d — map should be analyzed according to the variation of parameter d and
the determinant A = p? — 4q. The range of the cyclical behavior is determined by the
classical imaginary solution of the dynamical system which is in our case defined by the
characteristic equation:

—2b+ d+ V—4b? + d?
A= (50)
2b

Stability result corresponds to the polynomial 4? = tr.J\ — det J where periodic solutions
will be considered. One should consider e.g. [5] for details and research in the field of
discrete dynamics cobweb models [28, 13]. The ordering of the equilibriums is determined
by the general Eq. 47. The rational fraction , which is in our case transformed by
the Eq. 47 to the value of the parameter d, corresponds to the Farey sequence which
could be represented by the Farey tree. Important to notice is the sequences of Fibonacci
numbers [17] {F;/Fii1,i € 8} and {F;/Fi,2,7 € 8}. Since the observed system is linear
the change in parameter b will not influence the periodicity condition. The emergence of
the system periodic stability in the shape of n-sided polygon could be observed not only
in economical systems [21, 20]; the n-sided polygon and the Farey tree organization of the
equilibria could be observed in the technical systems as for example in laser control as the
paradigin of the chaotic system [23].

18
4.1 Bifurcation Analysis - Extension to Nonlinear systems

Periodicity conditions in previous section are somehow general and could be transferred
from linear systems to nonlinear, see for example [5, 24]. In order to observe the pseudo-
bifurcation response of the initial 2-d discrete map the bifurcation was performed for
the change in parameter d in the range d € |—2,2] which covers the whole periodicity
area of the studied model. Fig. 10 represents pseudo bifurcation diagram for the linear
system. One could observe the strong periodicity points which are marked with the vertical
lines and correspond to the previously studied Farey tree sequence of polygons. X-azis
represents the variation in the parameter d and on the Y-axis the values of P;, are shown.
The graph shows the periodicity response as major gaps in the bifurcation with the pole
at the origin of X-azis.

rit

5
* ~12
3

5

Figure 10: Bifurcation diagram - initial linear 2d discrete map
The observation developed so far leads to the following proposition:

Proposition 4.1 Periodicity conditions in the linear 2-d map of the cobweb anticipative
system exist in the nonlinear expansion of the system.

If we consider the proposed socio-economic model by introducing the nonlinearity there
are several implications [27, 30, 31]: a) neccessity of the nonlinear equation application
in order to simulate the system, b) evolutionary character of the socio-economic systems
could only be revealed via computer simulation and c) nonexistence of analytical solutions.
The obvious approach to model complex dynamical systems is therefore the mixture of
several approaches, continuous simulation, discrete-event simulation and application of
nonlinear dynamics models. Important fact that should be considered at the modeling
of complex systems is the impossibility of prediction [27] as the immanent characteristic
which is present in nonlinear chaotic model representations. Consider generic alteration

14
of the initial anticipative cobweb model:

Pe(k+1) = P+ Prpilk)-

1
- (Fx)+ seam)

Peei(kt+1) = Prpo(k)
Prp(kt+1) = § (Pao - Pet) e+e)
P2o(k+1) = Polk) + Px(k)Pxpi(k) -
— vPz(k) (51)

Slight modification of initial Hicks’ model [22] with applied accelerator principle [34] gives
the interesting response. The system can be represented in three dimensions which reveals
the periodicity of the system for which the previously determined conditions of Farey tree
generally still holds. Fig. 11 shows the 3d bifurcation diagram for the altered model.
One can see the four attractors which are simultaneous and represent the four possible
equilibrium states for the trade dynamics. The 4-cycle characteristic is preserved at the
alteration of the parameter v which could be observed in the Fig. 12. The four dots on the
center-right side of the figure represents the four-cycle characteristic of the response. The
larger orbits indicate the steep change in the modus of the system. However the analytical
proof of the periodicity would be hard. The underlying Farey sequence define the adapted
nonlinear 2-d discrete map. Such evidences are also found in other works in nonlinear
system analysis for example [8] or in the recent works of Swedish economist T. Puu.

05 E
S Pi 05

Pk

Figure 11: left) Emergence of two Synchronous Attractors in the nonlinear case where
d = 0.3833 and b; = 0.33 and right) Emergence of Four Synchronous Attractors where
d = 0.160793 and 6; = 0.18

In order to analyze the preservation of the periodic solutions the most. significant
periodic solution i.e. the period 6 has been applied to the system Eq. 51 which is
restated in the following form:

d
Pepi = 7PKpi—Px+e—a
1
Pe = Pea-——
K KPL— Dope

15
Pkpt

1000 50

Pz ~500  -50 Pk

Figure 12: Preservation of Four Synchronous Attractors in the nonlinear case where d =
0.16404706 and b; = 0.18

Pz = Pz(l—v)+PrpiPK (52)

The following condition has been applied to the Eq. 52:

dPpy , d(-at+e+ Pe — Px)
a b

—Pxpi — +

d PK. 1
-a+e- Pep) +— +

1
(Pri — pep;) (Pepi Px + (1-0) Pz)
d PRP

—(Ppi Px) + (-at+ e+ —[— — Px)

1
(Prpi- Py Py? —(l-v)Pzt+

(1—v) (Pxpi Px + (1—v) Pz) =0

which yields the periodicity solutions that are related to the system attractor. For ex-
ample, for period 10 the initial values are:Px p) = —1, Px = —1.61803, Pz = 1.61803, d =
1.618033989, v = 1, for period 7: d = cos ay = 1, Pep; = —2.16007, Pe = —0.533493, Pz =
1.15238.

For the nonlinear mapping symmetry could be expected as shown with numerical
examination of the system for the parameter d € |0.3833, 0.4538

As the illustration of the presented concept and it’s relation to the real world [39] the
following example will be considered: Fig. 14 shows the real world data for the Number
of births, daily, Quebec, January 1, 1977 to December 31, 1990 [25]. Other series closer to

16
14

i
a5 -1 05 0 05 1 15

Figure 13: Symmetry of the nonlinear map for d € [0.3833, 0.4538]

the Supply~Demand system could also be considered here. This time series shows nothing
in particular. However, if we consider the data by the Poincaré first-return map [35] as
shown in Fig.15 one could observe the strong quad (>) polygon showing the basic periodic
characteristic of considered system.

500

450

Number of births
yp 8 @ Be
3.3 8 8

x
Ss

500 1000 1500 2000 2500 3000 3500 4000 4500 5000
[day] from 1.1. 1977

Figure 14: Real world system response - number of births from 1.1.1977

The periodically driven system is identified in the Fig.15 where the strong 4-cycle
resonance is shown. One of the common periodicity conditions is 6-cycle as the cycle of
possible optimal response.

If we consider the proposed nonlinear system by the examination of Lyapunov expo-

17
0.15,

ot

0.05)

b(k+1)

0.05

04 0.05 0 0.05 of 0.415

Figure 15: Poincaré first-return map

nents defined by Eq. 53 the chaotic property is indicated [36]

= (inf (eo)] + In| f'(ea)| +... +In|f"Hn)]) (53)

Eq. 53 tells us that the Lypunov exponent is the rate of divergence of the two trajecto-
ries. Stated differently, we have here the average of the natural logarithm of the absolute
value of the derivatives of the map function evaluated at the trajectory points. If the ap-
plication of the map function to two nearby points leads to two points further apart then
the absolute value of the derivative of the map function is greater than 1 when evaluated
at those trajectory points. If the absolute value is greater than 1 then the logarithm is
positive. If the trajectory points continue to diverge on the average, then the average of
the logarithm of the derivatives absolute values is positive. In our case the Gram-Schmidt
orthogonalization procedure [37, 38, 18] which gives the following results for the three
Lyapunov exponents, \1, Az and Ag shown in Fig.16.

From the Fig.16 the stable regions where all three the Lyapunov exponents 41, Az and
A3 are less than 0 could be observed. For the studied system the regions of stability are
important. As the example of the application of Lyapunov exponents one should consider
parameter bifurcation in the parameter space of interest in order to identify regions of
stability. As an example, the computation of the Lyapunov exponents for parameters d €
[0.1, 4] and v € [0.5, 0.6] was performed indicating the region of stability shown in the Fig.
17. The blue region shows the chaotic parameter space of the system while the white region
indicates the stable region. Such stability regions are usually periodic regions obeying
previously mentioned Farey sequence rule. At the analysis of nonlinear system the regions
of stability are of particular interest [33, 32, 10, 9, 14]. This regions are dependant on the
parameter values and initial conditions; here we have the function f(A, Xo, Yo, Zo) where
A represents set of parameter values, Xo, Yo and Zp represents the set of initial values.
Analytical determination of boundary conditions of periodicity regions is difficult task and

18
02

Lyapunov Exponent

of 0.15 0.2 0.25 03 0.35 04
bifurcation parameter

Figure 16: Lyapunov exponents for 3-d discrete nonlinear map for 1, A2 and A3

in many cases impossible. For the determination of the stability regions, the parameter
space search is performed as the aggregate Lyapunov estimator dependant on the system
response for the function f(A, Xo, Yo, Zo). However, one should be aware, that a full
understanding of a nonlinear system might require computing solutions for essentially all
initial conditions {26]. Since this is infeasible, we would like some simple ways to summarize
the possible behaviors as proposed by the Farey tree sequence. Here the application of
powerful parallel computing provides the basic platform for systematic [40, 41] nonlinear
system analysis. But one should be warned that only the analytical explanations of the
nonlinear systems’ behavior provide methodologically correct characterization

4.2 Equilibrium Analysis

Equilibrium condition for the P segment of the hyperincursive cobweb system could be
stated as:

(2-8 _ “a eh 8) _?
b\ 6 d ~

(54)

The equilibrium values of the parameters for the P segment of the system are: a =c
and h=d.

Equilibrium condition for the Q, segment of the hyperincursive cobweb model could
be stated as:

Dp

§(p—a-5(a+ 4-0) -0)) = pat 0) (55)

Q, segment of the system has no solution for the equilibrium values of the parameters.
When the equilibrium conditions for the P segment of the system are considered in all the

19
Figure 17: Lyapunov exponents stability region (white) for 3-d discrete nonlinear map for
parameter space d € [0.1,4] and v € (0.5, 0.6]

cases, the Q, segment of the system could not be in stable state. Graphical presentation
of the equilibrium conditions a = p and b = d is shown in Fig. 18. The hexagon
shape is known as the possible optimal shape in the context of spatial economics [21].
The existence of hexagonal shapes is known - in the space economy and explained by the
structural stability [21]

Qs (kH1)

Qs (k)

Figure 18: Response of the Q, segment of the system while in equilibrium

Proposition 4.2 Equilibrium condition for the P segment of the anticipative cobweb sys-
tem defined by the equations from Eq. 29 to Eq. 84 is: a=c=pandb=d. At this
conditions, the Q, segment of the system has the response of hexagon shape with vertices
{(a,0), (a, a), (0, a), (—a, 0), (—a, —a), (0, -—a)} in Qs(k), Qs(k +1) mapping.

While the response of the system for the P segment is in equilibrium, the Q, segment
of the system has the hexagon-like shape with significant dimension of parameter a value
and edges dimension of a and ay2.

20
Proposition 4.3 Triangular (A\) i.e. three-period response in 2-dimensional mapping is
determined by the condition b = —d.

In order to gain the term for the period in the P values one should apply Prop. 2.1. The
values for times 1,2,3,4 should be symbolically expressed. By inserting the Eq. 29 and
Eq. 30 into Eq. 33 the following term is gained:

b(p-c)
p-a a-—ct ——)
b d (56)

P(k+2)=$(

By repetition of similar procedure the equation for P(k + 3) considering the period 4
condition i.e.: P(k +3) = P(k) one should get the following equation:

_ _ b(p—c) _ _
GC @ a-c+ 5 2 :\ fo (67)

b\bo\ 6b d d
with solution b = —d € R/{0}.

Proposition 4.4 Hexagon (©) i.e. sir-period response in 2-dimensional mapping is
determined by the condition b = d.

The symbolic solution is gained by Def.2.1. The values for times 1,...,6 should be
symbolically expressed. Determination of P(k + 2) is based on expressions for P(k + 1)
and P(k), determination of P(k + 3) is based on expressions for P(k + 2) and P(k + 1)
etc. Periodic condition is expressed as P(k + 6) — P(k) = 0 with solution b = d € R/{0}

The general periodicity solution has been determined by the z-transform analysis. For
parameter d the solution is:

dees (58)

n

where n is the period and m = 1,2,3,...,.n —1. Here, as mentioned, the above procedure
could be performed for the arbitrary period n.

5 Discussion

The story revealed from the developed hyperincursive model raises the following questions:
a) Does a change in the strategy changes the structure or does it change only the relations
between the elements of the structure? b) Does changing of the strategy change the future
as well as the past?

A change in the strategy would mean a new and different future and should also mean a
different past if the change in the stratezy would occur earlier. The hyperincursive cobweb
model enables us to change the future as well as the past chain of events. However, differ-
ent examination of the system dynamics is proposed where change in the key parameters
is performed while observing the change in complete future and past chain rather than
observing the classical time response of the system.

The following procedure proposition emerges which enabled the anticipative formulation
of the classical dynamic system: since the hyperincursive systems are hard to determine
[7, 6], the developed anticipatory mechanisms should be applied, therefore the model
should be a) transformed in the separated form b) provide past-future chain property and
c) apply the hyperincursive structure to the studied model.

21
The developed model shows that by the statement of general rule of the system, the
synchronization of entire feedback-anticipative chain could be gained by setting the ap-
propriate strategy in the form of parameters value set which should be time dependant.
The idea for the simulation proposed in the paper is quite different from the common
paradigm. The structure of the model should yield the entire feedback-anticipative chain
and the observation of the entire system response should be made. This provides new
and quite challenging responses which should initiate further interest and examination of
proposed model.

One of the interesting responses from the model is synchronization of the model at the
certain time steps. Entire feedback-anticipative chain i.e. all point set, is synchronized
according to the period of the system. The solution of the periodicity conditions for the
2-d discrete linear cobweb map provided the means to determine the periodicity condi-
tions. Analytical approach with z-transformation provides the proper way to determine
the periodic solutions. The emergence of Farey tree as the rational fraction representa-
tion yields the organization of the periodicity solutions. The developed model shows that
by the statement of general rule of the system, the synchronization of entire feedback-
anticipative chain could be gained by setting the appropriate strategy in the form of
parameters value set which should be time dependant. The bifurcation experiment with
the nonlinear mapping provided the example of periodicity transposition to the systems
of higher complexities. Period 6 has been determined as one of the most stable periodic
solution, which has been explicitly shown by the analysis of system’s attractors

Acknowledgement
This research was supported by the Ministry of Higher Education, Science and Technology of the Republic of Slove-
nia (Program No. UNI-MB-0986-P5-0018).

References

[1] Shone. R. (1997) Economic Dynamics. Cambridge University Press

[2] Rosser. M. (2003) Basic Mathematics for Economists, 2"4ed. Routledge

[3] Azariadis C., (1993). Intertemporal Macroeconomics. Blackwell Publishers. Cambridge MA.

[4] Luenberger, D. G. (1979) Introduction to Dynamics Systems \ Theory, Models & Applications. John & Wiley
Sons Inc.

[5] Sonis, M. (1999). Critical Bifurcation Surfaces of 3D Discrete Dynamics. Discrete Dynamics in Nature and
Society, Vol. 4, pp. 333-343. OPA Overseas Publishers Association N.V,

[6] Dubois D., Resconi G. (1992) Hyperincurisvity — a new mathematical theory, 1°ted. Presses Universitaires de
Liége

[7] Rosen R. (1985) Anticipatory Systems. Pergamon Press

[8] Brock, W. A., Hommes, C. H. (1997) A rational route to randomness. Econometrica. Vol. 65, No. 5, 1059 -
1095

[9] Gallas, J., Nusse, H. E. (1996) Periodicity versus chaos in the dynamics of cobweb models. Journal Of Economic
Behavior & Organization. Vol. 29, 447 — 464.

[10] Hommes, C. H. (1998) On the consistency of backward-looking expectations: The case of the cobweb. Journal
Of Economic Behavior & Organization. Vol. 33, 333 — 362.

[11] Kijajié M. (1998) Modeling and Understanding the Complex System Within Cybernetics. Ramaekers M. J
(Ed.), 15th International Congress on Cybernetics. Association International de Cybernetique.
Namur. 864 — 869

[12] Kljajié M. (2001) Contribution to the meaning and understanding of anticipatory systems. Dubois D. M. (Ed.)
Computing anticipatory systems, (AIP conference proceedings, 573). Melville (New York): American Institute
of Physics, 400 — 411

[13] Kljajié, M., Skraba, A., Kofjaé, D. Bren, M., (2005) Discrete Cobweb Model Dynamics with Anticipative
Structure. WSEAS Tranactions on Systems. Vol. 4, Issue 5.

[14] Matsumoto, A. (1997) Ergodic cobweb chaos. Discrete Dynamics in Nature and Society . Vol. 1, 135 — 146

22
15] Kreyszig, E. (1993) Advanced Engineering Mathematics. 7" ed, John & Wiley Sons Inc.
(15) (1993) Advanced Math *" ed. John & Wiley Sons I

[16] Strogatz S. H. (1994). Nonlinear dynamics and chaos: with applications to physics, biology, chemistry and
engineering. Addison-Wesley Publishing Company.

[17] Kocic L. M., Stefanovska L. R., (2004). Golden Route to Chaos. Facta Universitatis. Series: Mechanics, Auto-
matic Control and Robotics. Vol. 4, No. 16., pp. 491-496.

[18] Steeb W-H. (1999). The Nonlinear Workbook. World Scientific Publishing. Singapore.
[19] Wiener N. (1961). Cybernetics or Control and Communication in the Animal. MIT Press: Cambridge, MA

20] Puu T., Sushko I. (2004). A business cycle model with cubic nonlinearity. Chaos, Solitons and Fractals 19 pp.
AQ7612.

[21] Puu T. (2005). On the Genesis of Hexagonal Shapes Tnu Puu. Networks and Spatial Economics. Vol.5, Iss. 1
pp. 5-20.

[22] Puu T. (1963). A Graphical Solution of Second-Order Homogeneous Difference Equations. Oxford Economic
Papers. Mar 1963; 15, 1. pp. 53-58.

[23]. Kociuba G.,Heckenberg N. R., White A. G. (2001). Transforming chaos to periodic oscillations. Physical Review
E, Vol. 64, No. 5. DOL: 10.1103/PhysRevE.64.056220.

[24] Sonis M. (1996). Linear Bifurcation Analysis with Application to Relative Socio-Spatial Dynamics. Discrete
Dynamics in Nature and Society. Vol. 1. pp. 45-56

[25] Hyndman, R.J. (n.d.) Time Series Data Library, http://www.robhyndman.info/TSDL/. Accessed on
27!" Feb. 2006.

[26] Cull P., Flahive M., Robson R. (2005). Difference Equations - From Rabbits to Chaos. Springer.

[27] (1990) Radzicki, M.J. Institutional Dynamics, Deterministic Chaos and Self-Organizing Systems, Journal of
Economic Issues, Vol. XXIV No.1 March 1990, pp. 57-102

[28] Skraba, A., Kljajié, M., Kofjaé, D. Bren, M., Mrkaié M. (2005). Periodic Cycles in Discrete Cobweb Model.
WSEAS Transactions on Mathematics. Issue 3, Vol. 4, July 2005. pp. 196-203.

[29] Pavlov O. V., Radzicki M, and Saeed K. (2005). Stability in a Superpower-Dominated Global Economic System,
Journal of Economic Issues. Vol. XXXIX, No. 2, June 2005. pp. 491-500.

[30] Radzicki M. J. (1993). Methodologia oeconomiae et systematis dynamis. System Dynamics Review 6 No. 2.
Summer. pp. 123-147.

[31] Sprott J. C. (2004). Dynamical Models of Love. Nonlinear Dynamics, Psychology and Life Sciences, Vol. 8, No
3. Society of Chaos Theory in Psychology & Life Sciences.

[32] Guerra J. H., Rodriguez F.F. (1997). A qualitative study of the disagregated long-wave model. System Dynamics
Review. Vol 13, No.1. pp. 87-96.

[33] Sice P., Mosekilde B., Moscardini A., Lawler K. and French I. (2000). Using system dynamics to analyse
interactions in duopoly competition. System Dynamics Review. Vol. 16, No. 2. pp. 113-133.

[34] Moscardini A.O., Lawler K. and Loufti M. (1998) Systemic Insight into Economic Principles. 16th International
Conference of the System Dynamics Society.

[35] Berry B.J.L. and Kim H. (2000). Long Waves 1790-1990in Chaos Theory in The Social Sciences. Kiel L.D and
Elliott E. (Eds.)

[36] Hilborn R.C. (2000). Chaos and Nonlinear Dynamics. Oxford University Press
[37] Alligood K.T., Sauer T.D. and Yorke J.A. (1996). Chaos - An introduction to Dynamical Systems. Springer.

[38] Govorukhin V.N. MATDS — MATLAB-based program for dynamical system _ investigation.
http://www.math.rsu.ru/mexmat/kvm/matds/ Accessed on 28'” Feb.2006.

[39] Peterson D.W. and Eberlein R.L. (1994). Reality Check: a bridge between systems thinking and system
dynamics. System Dynamics Review. Vol 10, no. 2-3. pp. 159-174.

[40] Sprott J.C. (1994). Some simple chaotic flows. Physical Review E. Vol. 50, No. 2. APS.

[41] Chechin G.M. and Ryabov D.S. (2004). Three-dimensional chaotic flows with discrete symmetries. Physical
Review F. Vol. 69. APS.

23