1_9788490866863 by Editorial Tirant Lo Blanch

An introduction to economic dynamics

Clara Calvo Prof. Doctora de la Facultad de EconomĂa Universidad de Valencia Carlos Ivorra Prof. Titular de la Facultad de EconomĂa Universidad de Valencia

Copyright ® 2015 Todos los derechos reservados. Ni la totalidad ni parte de este libro puede reproducirse o transmitirse por ningún procedimiento electrónico o mecánico, incluyendo fotocopia, grabación magnética, o cualquier almacenamiento de información y sistema de recuperación sin permiso escrito de los autores y del editor. En caso de erratas y actualizaciones, la Editorial Tirant lo Blanch publicará la pertinente corrección en la página web www.tirant.com (http://www.tirant.com).

© TIRANT LO BLANCH EDITA: TIRANT LO BLANCH VALENCIA TELFS.: 96/361 00 48 - 50 Email: tlb@tirant.com http://www.tirant.com Librería Virtual: http://www.tirant.es DEPOSITO LEGAL: V-711-2015 ISBN 978-84-9086-686-3 MAQUETA: Si tiene alguna queja o sugerencia, envíenos un mail a: atencioncliente@tirant.com. En caso de no ser atendida su sugerencia, por favor, lea nuestro Procedimiento de quejas en: www.tirant.net/index.php/empresa/politicas-de-empresa

Chapter I Introduction to Dynamical Analysis Dynamical analysis is a branch of mathematics which studies how one or more quantities evolve over time. Besides its purely mathematical interest, it has a wide range of applications, from theoretical physics to social sciences. Here we are, of course, particularly interested in its applications to Economics, i.e., in what is called Economic Dynamics. In this introductory chapter we will use some illustrative examples to present an overview of the economic situations where it can be applied and the kind of conclusions that they allow us to draw.

Ÿ Continuous and Discrete Dynamics The main underlying idea of dynamical analysis is that there are many situations in which the current value of one or several economic magnitudes determines in a somewhat simple way their variation in the short term, and then an adequate mathematical analysis allows us to draw conclusions of what we can expect in the mid or long term. There are essentially two ways of expressing the expected value of a magnitude in the short term as a function of its current value. The simplest one can be applied when we can assume that the quantities we are concerned with can change only after a given fixed amount of time (each month, each year, etc.). In this case we can consider the value x0 that a magnitude takes in the initial time t = 0, the value x1 it takes at the next period t = 1, the value x2 that it takes at the next period t = 2 and so on. Then we have what is called a time series, i.e., a sequence x0 , x1 , x2 , x3 , ... , xt , xt+1 , xt+2 , ... representing the values of a certain magnitude at uniform time intervals (days, months, years). Many time series interesting from an economic point of view are so complex that they must be handled by means of statistical techniques (consider, for instance, the series of daily closing values of the Dow Jones Industrial Average), but this is not the case we are going to consider. Here we will be interested in those time series for which the value xt+1 is determined by the previous value xt by a more or less simple (deterministic) mathematical law. This kind of laws determine what it is called a discrete dynamical system. (This is the simplest case. Later we will consider a more general setting.) 1. Example: Compound interest Perhaps the most simple discrete dynamical system in Economics is the compound interest. Namely, assume that a saver deposits an amount of money in an interest-bearing account providing an annual nominal interest rate j compounded monthly. If we call St the balance at the t-th month after the deposit was made, this means that St+1 - St =

j St

(1)

12 This is our first example of what is called a difference equation, i.e., an equation about a quantity St (the balance of the account at the time t) that does not provide the value of St for each time t, but just relates each value St with the difference DSt = St+1 - St . (This is not the general definition of “difference equation”, we will state it later.) Solving this equation means obtaining an explicit expression for the value St in terms of t. In this case the solution is very simple: (1) can be transformed into St+1 = St H1 + j • 12Land hence, if we call the principal (initial investment) P = S0 , it is clear that S1 = P H1 + j • 12L, S2 = S1 H1 + j • 12L = P H1 + j • 12L2 , S3 = P H1 + j • 12L3 , ... and, in general: St = P H1 + j • 12Lt

(2)

Econ om ic D yn a m ics

In spite of its simplicity, by comparing the difference equation (1) and its solution (2), we can observe a fact that, as we will see, holds in a very general setting: The difference equation does not have a unique solution, i.e., there is not a unique time series satisfying it. On the contrary, there is an infinite number of such time series, one for each possible value of the principal P. More specifically, for j = 0.01, here we have two different solutions of the same difference equation 10200.00 10302.00 10405.02 10509.07 10614.16 10720.30 ... 3700.00 3737.00 3774.37 3812.11 3850.23 3888.74 ...

corresponding to two different principals P = 10 200 and P = 3700, respectively. We say that (2) is the general solution to the given difference equation, and that each time series obtained by fixing a value for the principal is a particular solution to the difference equation. Notice that if we change the value of j we also obtain a different time series, but, since j appears in the difference equation, this is not a solution to the same equation, but to another one. On the contrary, P does not appear in the equation, and so, when we change the value of P, we go from one solution to another solution for the same equation. In general, we will see that difference equations have general solutions involving one or more arbitrary constants (like P) that, when particularized to concrete values, provide an infinite number of particular solutions. 2. Example: Continuous interest The GDP (Gross Domestic Product) per capita of USA in 2010 was 46 612 $ and until 2012 it experimented an average growth of 3.148 %. What does it mean? Certainly, it does not mean that the GDP per capita in 2012 was GDP2012 = 46 612 H1 + 0.03148L2 = 49 592.88 $ In fact, the GDP per capita in 2012 was 49 641 $. If the average growth had to be calculated by the compound interest formula, it would be the solution of the equation 49 641 = 46 612 H1 + i1 L2 , namely, i1 = 3.198 % instead of 3.148 %. The difference is small, but it has its own significance. The growth of the GDP is not measured this way, since then we would be assuming that the growth is compounded once a year, and this is not true: the GDP of a country is varying continuously, and hence we cannot speak just about a time series GDP0 , GDP1 , GDP2 , ... with a single value for each year, but instead we have a function GDPHtL providing the GDP per capita of the country at each time t, where t can take the value t = 0 (for Jan, 1, 2010), t = 1 (for Jan, 1, 2011), but also t = 0.8 (for Oct, 19, 2010). The variation of this function is expressed by its derivative: dGDP

= 0.03148 GDP.

dt Here the left hand side term expresses the rate of growth of GDPHtL, and we are saying that this rate has been the 3.148 % of the current GDP HtL at any time in the period 2010-- 2012. More accurately, we are talking about the average growth rate, and this means that the rate might have been different at each instant, but the growth of the GDP per capita in the whole interval has been the same that we would have got if the growth rate had always been equal to this mean value. This a simple example of a differential equation, a particular case of the following one: dS

=iS

(3)

dt It is a differential equation because it relates an unknown function SHtL to its derivative dS â€˘ dt. If we know SHt) for a certain time t, we know how S is growing (the speed of the growth) at that time t. (As in the case of difference equations, we will provide a more general definition later.) Solving this differential equation means obtaining the explicit function providing the value for each time .

Ch a pt e r I : I n t r odu ct ion t o D yn a m ica l An a lysis

case of difference equations, we will provide a more general definition later.) Solving this differential equation means obtaining the explicit function SHtL providing the value S for each time t. The functions SHtL satisfying (3) are said to be growing at a constant rate i of continuous interest. Notice the perfect analogy between the continuous and the discrete case, where the derivative appears as the continuous analogue to the discrete difference: d S(t)

= i S(t) , DSt =

j St 12

The factor 12 appears simply because in the discrete case we were expressing the time in months and in the continuous one we are working in years. Solving (3) is not as obvious as solving its discrete counterpart (1), but it is still easy : dS S

= i dt • à

= à i dt • ln(S)-ln(S0 L = i(t-t0 ) • ln(S) = ln HS0 L + i Ht - t0 L t

• S = eln(S0 )+i(t-t0 ) = eln(S0 ) ei(t-t0 ) = s0 ei(t-t0 ) , where S0 = S(t0 ) . Taking t0 = 0 for simplicity and calling P = S(0) , we obtain the formula S = P eit ,

(4)

which is the general solution of the differential equation (3). The particular solutions appear when the constant P is given a fixed value. The number e = 2.718281828459045 ..., which is ubiquitous in many mathematical areas, as well as in physics and other sciences, receives its name from Euler, but it was in fact discovered by the mathematician Jacob Bernoulli by studying a question on compound interest. If a principal P = 1 is invested at a 100 % of interest per year compounded at the end of the period, the final balance is S = 2, but if the interest is compounded n times in the year, the final balance will be: 1 n . SnL = 1 + n This value increases with n, but it does not tend to infinity. On the contrary, the maximum final balance that can be reached by increasing the number of times the interest is compounded is S¦ = lim 1 + n¯¦

, n

and this is the definition of e. From this definition we see that for a principal P and a nominal annual interest rate j the maximum final balance that can be reached after t years as the number of compounding periods grows is tnj

S¦ = lim PJ1 + n¯¦

j n

= P lim 1 + n¯¦

1 n j

= P lim 1 + n¯¦

= P etj ,

n j

and this fact provides an alternative interpretation of continuous compounding as the limit of the compound interest when the number of times the interest is compounded in a fixed period tends to infinity.

All in all, previous examples illustrate the two kinds of dynamical systems we will consider here: è Discrete dynamical systems in which the time variable t is restricted to integer values 0, 1, 2, 3, ... , (negative values are also alowed) and in which the dynamics are expressed by means of difference equations. è Continuous dynamical systems in which the time variable t can take arbitrary real values, and in which the dynamics are expressed by means of differential equations. Since the techniques for solving difference and differential equations are useful in a quite reduced number of cases, we will mainly concentrate ourselves in analyzing the solutions obtained with the aid of computer packages such as Mathematica.

Econ om ic D yn a m ics

computer packages such as Mathematica.

Ÿ Introduction to Mathematica Mathematica is a computational software program that can perform a wide range of mathematical computations. Here we will consider just those features related with the contents of this course. Let us recall first the very basic syntactical notions. The content of a Mathematica notebook is divided in cells. Each cell can contain text (like this one), one or more mathematical orders (input cells), or the answers provided by Mathematica to some previous orders (output cells). Here is an example of a simple input cell. By placing the cursor on it and pressing Shift-Enter (or Enter on the numeric key pad), the corresponding output will be generated: 2 + 3 5

Mathematica is case sensitive. All reserved words start with an uppercase letter (Sin, Plot, Pi, etc.). Because of this, even mathematical constants such as e or i, which are usually written in lowercase, are represented in Mathematica as E or I, respectively. Evaluate, for instance, the following cell: E^2 ä2

When this is evaluated, Mathematica writes ä2 , which is the output formated form corresponding to the plain text entry E^2. To force Mathematica to provide a numerical value of any expression we use the function N: N@E ^ 2D 7.38906

Alternatively, you can force Mathematica to work with decimal numbers by providing a decimal number in the input. E ^ 2. 7.38906

The symbol = defines assignations, whereas == defines an equation. For instance: x = 2 H*This defines the variable x to be iqual to *L 2 x ‹ 2 H*This is the statement x = 2, which is true*L y ‹ 2 H*This is the statement y = 2, whose truth value is unknown *L 2

True

y‹2

Parentheses are used to indicate as usual the priority of operations:

Ch a pt e r I : I n t r odu ct ion t o D yn a m ica l An a lysis

H2 + 4L 3 H*Comments can be inserted everywhere this way *L 18

The arguments of functions or commands are always enclosed in square brackets: Cos@PiD N@PiD Exp@2D -1

3.14159

채2

The usual brackets define arrays. The next cell contains the empty array, a 3-component array and finally a 3 x 3 matrix, which is represented just as an array of arrays A1 = 8< A2 = 83, 6, 1< A3 = 881, 2, 3<, 83, 2, 1<, 81, 0, - 1<< 8< 83, 6, 1< 881, 2, 3<, 83, 2, 1<, 81, 0, - 1<<

A matrix can be visualized in the usual way by means of the command MatrixForm: MatrixForm@A3D 1 2 3 3 2 1 1 0 -1

In order to make reference to a component of an array or matrix use double square brackets: A2@@1DD A3@@3, 2DD A3@@3DD@@2DD H*This is equivalent to previous expression *L 3

Mathematica recalls the value assigned to each variable within a session. In order to clear those values we can use the command Clear:

Econ om ic D yn a m ics

x=3 x Clear@xD x

H*This H*This H*This H*This

assigns x the value 3*L prints the value 3*L makes Mathematica forget the value of *L x just prints x since now x has no known value *L

Besides the default built-in functions, new ones can be defined this way: F = Function@8u, v<, u ^ 2 vD F@2, 3D FunctionA8u, v<, u2 vE 12

The discrete equivalent of Function is Table, which defines arrays of objects satisfying a given condition: Table@i ^ 2, 8i, 1, 5<D Table@8i, i + 1<, 8i, 4, 10<D 81, 4, 9, 16, 25< 884, 5<, 85, 6<, 86, 7<, 87, 8<, 88, 9<, 89, 10<, 810, 11<<

We will also need to know a kind of structures called “rules”. A rule is an expression like: t ^ 2 ¯ y + z H*The arrow is typed as ->*L t2 ¯ y + z

Rules can be named, as any other object: R = t^2 ¯ y + z t2 ¯ y + z

The command /. is used to apply a rule, and its effect is the replacement of the left hand side by the right hand side: 8t, t ^ 2, t ^ 3< •. R H*Replace t2 with y+z*L Hx + 2 yL ^ 3 •. y ¯ Sin@yD H*Replace y with Sin@yD*L 9t, y + z, t3 = Hx + 2 Sin@yDL3

Ch a pt e r I : I n t r odu ct ion t o D yn a m ica l An a lysis

A set of rules can also be applied at once: p + q + r •. 8p ¯ x ^ 2, q ¯ y ^ 3, r ¯ z ^ 5< x2 + y3 + z5

For instance, Mathematica renders the solutions of equations or systems of equations by means of rules: SolveA7 + 6 x - 8 x2 + x3 ‹ 0, xE :8x ¯ 7<, :x ¯

y SolveB: 1 2

J1 -

+ 3 ‹ y,

x ::x ¯ -

y 2x

5 N>, :x ¯

1 x

1 2

J1 +

5 N>>

‹ 1>, 8x, y<F

, y ¯ 1>, 8x ¯ 2, y ¯ 6<>

A final semicolon causes a line to be evaluated, but the result of the evaluation is not shown. The % symbol takes the value of the last result: 881, 2, 3<, 81, 1, 1<, 83, 1, 4<<; H*This enters a matrix, but no ouput is generated *L Inverse@%D; H*Calculates the inverse of the matrix defined before *L MatrixForm@%D H*Shows the last result in matrixform mode *L %% H*This refers to the penultimate result, in this case the inverse matrix in array mode *L 3 5 1 5 2 5

::-

1 1 -1

1 , 1,

1 5 2 5 1 5

1 2 2 1 >, : , 1, - >, : , - 1, >> 5 5 5 5 5

Econ om ic D yn a m ics

Ĺ¸ Graphing functions The basic graphing commands are Plot and Plot3D: Row@8Plot@2 x H1 - xL, 8x, 0, 1<D,

Plot3D@Sin@x yD, 8x, 0, Pi<, 8y, 0, 2 Pi<D<D

0.5 0.4 0.3 0.2 0.1

0.2

0.4

0.6

0.8

1.0

The command Plot admits several functions at once: 3

PlotB:

x ,

x >, 8x, 0, 1<F

1.0

0.8

0.6

0.4

0.2

0.4

0.6

0.8

1.0

Option PlotRange selects the range shown in the graph: Row@8Plot@20 + x ^ 3, 8x, 0, 3<D,

Plot@20 + x ^ 3, 8x, 0, 3<, PlotRange ÂŻ 80, 50<D<D

50 45 40

35 30

10 0.5

1.0

1.5

ListPlot plots arrays:

2.0

2.5

3.0 0.0

0.5

1.0

1.5

2.0

2.5

3.0

Ch a pt e r I : I n t r odu ct ion t o D yn a m ica l An a lysis

Row@8ListPlot@82, 5, 3, 6, 9, 3<D,

ListPlot@82, 5, 3, 6, 9, 3<, Joined ¯ TrueD<D

Ÿ Differential equations with Mathematica In order to solve equation (3) in Mathematica we use the command DSolve: DSolve@S '@tD ‹ i S@tD, S@tD, tD H*Notice the double equal symbol *L 99S@tD ¯ äi t C@1D==

Here C[1] represents an arbitray constant. The answer means S(t)= Ceit . A particular solution is obtained this way: DSolve@8S '@tD ‹ 0.03148 * S@tD, S@0D ‹ 46 612<, S@tD, tD 99S@tD ¯ 46 612. ä0.03148 t ==

A technical issue is that the solution appears enclosed in double brackets. This is because DSolve can be used to solve systems of differential equations with multiple solutions. Then each solution is an array of functions and Mathematica renders an array of solutions. In our case we have an array comprising a single solution, which in turn is an array comprising a single function. Notice also that each function is rendered in the form of a rule which substitutes S[t] by the corresponding solution. A slight variant is DSolve@8S '@tD ‹ 0.03148 * S@tD, S@0D ‹ 46 612<, S, tD H*Notice the S instead of S@tD in the second argument*L 99S ¯ FunctionA8t<, 46 612. ä0.03148 t E==

Now the rule provides a function instead of the expression defining the function. Then we can give a name for the function itself: F = S •. %@@1DD; F@2D 49 641.

Let us analize the syntax in some detail:

Econ om ic D yn a m ics

Solution = DSolve@8S '@tD ‹ 0.03148 * S@tD, S@0D ‹ 46 612<, S, tD; H*Call Solution to the solution found *L Solution H*This is an array of an array of rules *L Solution@@1DD H*This is the first Hand uniqueL array of rules of Solution*L S •. Solution@@1DD H*Here we apply the array of rules —a single one, in fact— to replace S with the function *L F = S •. Solution@@1DDH*Here we assign the name F to that function *L F@2D H*Now we can use F to evaluate the solution at any point *L 99S ¯ FunctionA8t<, 46 612. ä0.03148 t E== 9S ¯ FunctionA8t<, 46 612. ä0.03148 t E= FunctionA8t<, 46 612. ä0.03148 t E FunctionA8t<, 46 612. ä0.03148 t E 49 641.

Now we can plot the solution : Plot@F@tD, 8t, 0, 2<, PlotRange ¯ 80, 50 000<D 50 000 40 000 30 000 20 000 10 000

0.0

0.5

1.0

1.5

2.0

There are many differential equations that cannot be solved symbolically. Try for instance: DSolve@8y '@xD ‹ y@xD Cos@x + y@xDD, y@0D ‹ 1<, y, xD DSolve@8y£ @xD ‹ Cos@x + y@xDD y@xD, y@0D ‹ 1<, y, xD

Such equations can be solved numerically, with NSolve. In this case the initial condicion is necessary, and we must specify the interval for the variable in which the solution must be calculated (here we have specified [0, 30]):

Ch a pt e r I : I n t r odu ct ion t o D yn a m ica l An a lysis

NDSolve@8y '@xD ‹ y@xD Cos@x + y@xDD, y@0D ‹ 1<, y, 8x, 0, 30<D G = y •. %@@1DD; Plot@G@xD, 8x, 0, 30<D 88y ¯ InterpolatingFunction@880., 30.<<, <>D<<

0.5 0.4 0.3 0.2 0.1 5

As a recapitulation, we have seen two commands for solving a differential equation, NDSolve (see the previous cell) and DSolve (as in the following cell): DSolve@8S '@tD ‹ 0.03148 * S@tD, S@0D ‹ 46 612<, S, tD; F = S •. %@@1DD; Plot@F@tD, 8t, 0, 2<, PlotRange ¯ 80, 50 000<D

Ÿ Mathematical modeling The two previous examples of dynamical systems are in some sense “true”, i.e., the compound interest law is the true formula that regulates the growing of a principal invested under some contractually specified financial conditions and, on the other hand, i = 3.14 % is the true average growth of the GDP per capita of USA in the specified period since it has been calculated to satisfy the formula GDP2012 = GDP2010 ei2 from the known values of the GDP. However, in most cases the dynamical systems appearing in Economics are attempts of guessing the unknown rules determining the evolution in time of certain economic magnitudes, and we cannot take them as “true” a priori. On the contrary, most mathematical models must be confronted with the empirical data and, even if they “work well”, we cannot expect that the theoretical results predicted by them fit the real behavior of the corresponding magnitudes exactly, because they are often convenient simplifications of the reality, in which only the most influential variables are being considered. Such simplifications are interesting since they allow us to understand and explain the essential behavior of the magnitudes under consideration. The divergence between a good model and the observed data can be interpreted as perturbations due to minor factors not reflected in the model. If these perturbations are too significant, a finer model including them must be devised. 3. Example: Malthusian vs. logistic growth In 1798, the British economist Thomas Malthus published his book An Essay on the Principle of Population, in which a theory about the growth of the human population in industrialized countries was proposed. In mathematical terms, Malthus stated that the population P HtL of laborers at time t grows in proportion to itself, namely: dPHtL

= r PHtL,

(5)

dt for certain demographical constant r mainly depending on the average number of children of each family. This is the same equation as (3) and we know that its solution is P(t) = P0 ert , where P0 is the population at t = 0. More specifically, Malthus estimated that each 25 years the population doubled its size, and this leads to the value r = 0.0277 (you can check it as an exercise). On the other hand, Malthus argued that the rate of increment of food production is constant, i.e.:

Econ om ic D yn a m ics

tion doubled its size, and this leads to the value r = 0.0277 (you can check it as an exercise). On the other hand, Malthus argued that the rate of increment of food production is constant, i.e.: dFHtL

= k F0 ,

where F0 is the quantity of food available at t = 0, and it is easy to see that the solution of this differential equation is F HtL = F0 H1 + ktL. Thus the quantity of food per capita is f(t) =

F(t)

F0 H1 + ktL P0 ert

P(t)

= f0 H1 + ktL e-rt .

A reasonable estimation for k is k = 0.05. The graph on the left represents f HtL for this value of k, whereas the graph on the right corresponds to k = 5: Row@8Plot@H1 + 0.05 tL Exp@-0.0277 tD, 8t, 0, 400<D,

Plot@H1 + 5 tL Exp@-0.0277 tD, 8t, 0, 400<D<D

1.0

0.8

40 0.6 30 0.4

0.2

10 100

200

300

400

100

200

300

400

We see that, even if the growth of the production of food in the world increased one hundred times faster than the real rate, the amount of food per capita would become zero after less than four centuries, according to the Malthusian model. Today, more than two centuries after the publication of Malthus' treatise, their predictions remain very far from the real data, and hence we can conclude that his model is not a realistic one or, in other words, that it misses some relevant aspects of the population dynamics. It must also be pointed out, however, that although the Malthusian model overestimates the population growth in the long run, it is a reasonable approximation which is widely used in many economic models for describing population growth as well as other indicators of economic growth, such as the Gross Domestic Profit, the national income, etc. In 1838, after having read Malthus' treatise, the belgian mathematician Pierre-François Verhulst published an alternative growing model which can be thought of as a correction of the differential equation (5). Namely, Verhulst assumes that there exists a carrying capacity K , which is the maximum population size that the environment can sustain indefinitely, but this does not mean that the population grows until reaching this maximum size to perish then of starvation. On the contrary, Verhulst model establishes that the rate of growth of the population is reduced as long as the number of individuals approaches to the carrying capacity. This condition is expressed in the differential equation d PHtL

= r PHtL 1 -

PHtL

(6)

where r is the same demographical constant appearing in Malthus' model. Notice that when P HtL approaches K , the second factor of the right hand side of (6) becomes close to 0 and the rate of growth of P (its derivative) becomes small, as we have indicated. This equation can be simplified by dividing both sides by K : d PHtL • K dt and calling x = P • K . Hence it becomes:

PHtL K

Ch a pt e r I : I n t r odu ct ion t o D yn a m ica l An a lysis

= r xH1 - xL

(7)

dt where now 0 ¤ xHtL ¤ 1 is the relative population. For instance, x = 0.25 means that the number of individuals is the 25 % of the carrying capacity. We can solve it with Mathematica and the result is: Clear@x, x0, rD; DSolve@8x '@tD == r * x@tD H1 - x@tDL, x@0D ‹ x0<, x, tD; x = x •. %@@1DD; x@tD x0 = 0.3; r = 1.2; Plot@x@tD, 8t, 0, 10<D är t x0 1 - x0 + är t x0 1.0 0.9 0.8 0.7 0.6 0.5 0.4 2

Equivalently: x(t) =

x0 ert

1 + x0 (ert - 1L

(8)

where x0 = x(0) is the initial relative population. Verhulst called this function the logistic function and (7) is called the logistic equation (although nobody knows why he considered that name adequate for such function). The figure shows the logistic function for an initial (relative) population x0 = 0.3 and a demographical constant r = 1.2. The Verhulstʼs model does not predict an indefinite growth leading to starvation. Instead, it is easy to show that lim t¯¦

x0 ert

1 + x0 (ert -1)

for each x0 , r > 0. This means that, unless the population is sterile (r = 0) it always grows asymptotically to the carrying capacity, but it does not exceed this limiting value. It is not clear that the growth of human population follows a logistic pattern, but the logistic growth model has found a great number of applications in biology, medicine (in the study of the growth of tumors), in phisics and many other sciences. In economics it has been applied among other contexts to describe the diffusion of innovations. The main conclusion that one should extract from the previous examples is that a mathematical model cannot be considered “true” a priori, but it must be tested. The exponential (Malthusian) growth and the logistic growth are two models of growth, and each one can describe some growth processes properly, but they can also fail to describe some others.

Econ om ic D yn a m ics

The Newtonʼs laws of motion, together with Newtonʼs gravitation law provides a system of differential equations that can predict the position of a planet in terms of its mass, the mass of the Sun and its initial position in an initial instant t0 . The predictions of such model are very accurate but, however, they are not exact, since they neglects the gravitatorial influence of the other planets. In 1821, the French astronomer Alexis Bouvard published a table of observed positions of the planet Uranus, and noticed that there were significant deviations from the values predicted by Newtonʼs laws, even taking into account the effect of the other planets. He suggested than they could be due to the presence of an unknown planet. In 1846, the French mathematician Urbain Le Verrier analyzed those differences and calculated what the orbit of the unknown planet should be in order to cause such perturbations in the orbit of Uranus. On 23 September, Le Verrier wrote to an astronomer of Berlin Observatory urging him to watch certain position in the sky and look for a moving object that could not be a star. That very evening the planet Neptune was discovered less than 1ë far from the position predicted by Le Verrier. Notice that, in order to describe the movement of a planet, we have a simple model which takes into account only the planet and the Sun, and other much more complicated models taking also into account the other planets. Those are more accurate, but this does not mean at all that the first simpler model is useless. In fact, this less accurate model is preferable in order to understand the movement of the planets, since it establishes that the orbit of a planet is essentially a consequence of the gravity of the Sun, whereas the gravity of the other planets just introduce small perturbations into the effect of the Sun. In general, a not too accurate model relying on fewer hypotheses (as the presence of the Sun) can be epistemologicaly better than a more accurate model relying on more hypotheses (as the presence of all planets), since the first one reveals that its few assumptions are the main relevant causes of the situation observed, whereas the other hypotheses are secondary ones that can be negligible in many contexts. When the objective is making predictions, the best model is the most complex and accurate, but when the objective is understanding, the best model is the simplest one providing enough accuracy.

Ÿ Difference equations with Mathematica Let us see now how Mathematica can be used to analyze the solutions of a difference equation. We consider an example which illustrates that choosing a continuous or a discrete model for describing a given situation can lead to dramatically different conclusions, and hence it may be a crucial decision that cannot be set rashly. Namely, let us study the discrete version of the logistic growth. 4. Example: Discrete logistic growth The continuous growth models can be adequate for describing populations that are growing continuously i.e., such that new individuals can appear at any moment. However, there are many biological species having a fixed annual breeding period, and so a discrete model that only considers the population Pt in the year t = 0, 1, 2, ... could be more accurate. The discrete analog of the Malthusian growth is described by the same equation as the compound interest, namely: Pt+1 = H1 + iL Pt = r Pt , where Pt is the population at time t. From this basic model, we can incorporate a logistic correction to the model by introducing a carrying capacity K : Pt+1 = r Pt 1 -

and dividing both sides by K and defining xt = Pt • K , this equation becomes the so-called logistic difference equation: xt+1 = r xt (1-xt ).

(9)

If we call f HxL = r x(1-x) , the equation has the form xt+1 = f(xt ) . Hence, given an initial (relative) population x0 , the corresponding time series is x0 , f(x0 ) , f(f(x0 )) , f(f(f(x0 ))) , ...

(10)

Ch a pt e r I : I n t r odu ct ion t o D yn a m ica l An a lysis

So it is convenient to define the n-th iterate of a function f , denoted by f nL , as the result of applying n times f to x , i.e., f 0L (x) = x, f 1L (x) = f(x) , f 2L (x) = f(f(x)) , f 3L (x) = f(f(f(x))) , ... In these terms, the time series (also called the orbit) corresponding to an initial state x0 is determined as xn = f nL Hx0 L. Mathematica can easily calculate these orbits with the command NestList: Clear@rD f = Function@8x<, r x H1 - xLD; r = 0.25; NestList@f, 0.9, 10D H*Calculate 10 iterations of f starting in 0.9 *L ListPlot@%, Joined ¯ True, PlotRange¯ 80, 1<D 90.9, 0.0225, 0.00549844, 0.00136705, 0.000341296, 0.0000852948,

0.0000213219, 5.33036 µ 10-6 , 1.33258 µ 10-6 , 3.33145 µ 10-7 , 8.32862 µ 10-8 = 1.0 0.8 0.6 0.4 0.2

We observe that starting near the carrying capacity (x0 = 0.9) the population tends to 0, whereas in the continuous case the limit of every orbit was always 1, i.e., the population always approached its maximum possible value. Let us change the demographic constant r: Clear@rD f = Function@8x<, r x H1 - xLD; r = 1.25; NestList@f, 0.9, 30D; ListPlot@%, Joined ¯ True, PlotRange¯ 80, 1<D 1.0

0.8

0.6

0.4

0.2

We can see that, for r = 1.25, the population stabilizes, but not in 1, as in the continuous case, but it tends to a value less than 0.2. The next panel shows the effects of changing the initial value and the demographic constant interactively:

Econ om ic D yn a m ics

r Initial value

1.0

0.8

0.6

0.4

0.2

We can see that the behavior of the dynamical system depends heavily on the value of r , and for values close to 4 the behavior is chaotic. The following graph shows the 100 first values for r = 3.96: Clear@rD f = Function@8x<, r x H1 - xLD; r = 3.96; NestList@f, 0.9, 100D; ListPlot@%, Joined -> True, PlotRange -> 80, 1<D 1.0

0.8

0.6

0.4

0.2

100

This chaotic behavior of the discrete logistic growth (for certain values of the demographical constant) does not appear in the continuous case, where the population always tends to the carrying capacity. However, notice that equation (9) is not exactly the discrete equivalent of (7), which would be D xt = r xt H1 - xt L â€” xt+1 = xt + r xt H1 - xt L. This equation makes sense for 0 < r < 1, and it is left as an exercise to check that its behavior is not chaotic at all. On the contrary, it perfectly parallels the behaviour of the continuous logistic model. On the other hand, we will see later that there are also chaotic continuous dynamical systems.

Ĺ¸ Introduction to Stochastic Dynamics There is no universally accepted mathematical definition of when the behavior of a dynamical system is chaotic. The underlying idea is that a dynamical system is chaotic when it does not follow any recognizable pattern, as in the case of the previous example (the discrete logistic growth for r = 3.96, as oposed to the non chaotic case r = 1.25, in which the orbit tends to stabilize). In some sense, we can say that a chaotic dynamical system takes random values, but this must be precised, since there is a very different

Ch a pt e r I : I n t r odu ct ion t o D yn a m ica l An a lysis

to the non chaotic case r = 1.25, in which the orbit tends to stabilize). In some sense, we can say that a chaotic dynamical system takes random values, but this must be precised, since there is a very different kind of random dynamical systems that should not be confused with the chaotic ones. These are the stochastic dynamical systems, but before introducing them let us try to especify in what sense a chaotic dynamical system can be considered as random. ã Chaos and randomness Let us consider a finite sample 8xn , ... , xm < of the orbit of the starting point x0 = 0.9 for the logistic map corresponding to the parameter r = 3.96. To follow the reasoning take n = 100, m = 120, but then we will consider larger intervals: r = 3.96; a = 0.9; n = 100; m = 120; f = Function@8x<, r x H1 - xLD; Sample = Drop@NestList@f, a, mD, n + 1D H*Drop deletes the first n+1 terms*L 80.561028, 0.975251, 0.0955799, 0.34232, 0.891542, 0.382911, 0.935709, 0.238224, 0.718634, 0.800708, 0.631916, 0.921089, 0.287829, 0.811734, 0.605174, 0.946196, 0.2016, 0.637392, 0.915249, 0.307169<

Now order the sample: Sample = Sort@SampleD 80.0955799, 0.2016, 0.238224, 0.287829, 0.307169, 0.34232, 0.382911, 0.561028, 0.605174, 0.631916, 0.637392, 0.718634, 0.800708, 0.811734, 0.891542, 0.915249, 0.921089, 0.935709, 0.946196, 0.975251<

Consider, for instance, the fifth term of the sequence, 0.307169. Since it is ordered, we can say that, when taking one of its terms at random, the probability of obtaining a number less than or equal to 0.307169 is 5 • Hm - nL = 0.25. Hence, the array: i

F>, 8i, 1, m - n<F m-n ListPlot@%, AspectRatio ¯ 1D

TableB:Sample@@iDD, NB

880.0955799, 0.05<, 80.2016, 0.1<, 80.238224, 0.15<, 80.287829, 0.2<, 80.307169, 0.25<, 80.34232, 0.3<, 80.382911, 0.35<, 80.561028, 0.4<, 80.605174, 0.45<, 80.631916, 0.5<, 80.637392, 0.55<, 80.718634, 0.6<, 80.800708, 0.65<, 80.811734, 0.7<, 80.891542, 0.75<, 80.915249, 0.8<, 80.921089, 0.85<, 80.935709, 0.9<, 80.946196, 0.95<, 80.975251, 1.<< 1.0

0.8

0.6

0.4

0.2

0.4

0.6

0.8

Econ om ic D yn a m ics

Represents a finite approximation of the function PHX ¤ xL giving the probability that a member of the sample is less than or equal to x. Now we do all these calculation at once with a larger sample: r = 3.96; a = 0.9; n = 100; m = 4100; f = Function@8x<, r x H1 - xLD; Sample = Sort@ Drop@NestList@f, a, mD, n + 1DD; ListPlot@Table@8Sample@@iDD, i • Hm - nL<, 8i, 1, m - n<D, AspectRatio ¯ 1, PlotStyle ¯ 8PointSize@0.004D<D 1.0

0.8

0.6

0.4

0.2

0.4

0.6

0.8

1.0

Changing the values of n and m we can see that the graph is essentially the same for large intervals. Thus we conclude that there is a function F (defined as the limit of the functions we are obtaining as the length of the sample tends to infinity) that allows us to calculate PHu ¤ X ¤ v L = F HvL - F HuL, i.e., the probability that a random point in the orbit of the fixed starting point is in the interval @u, v D. In other words, the sample can be considered as a sample of a random variable X defined by the cumulative distribution function F . This is not true for all values of r . The next panel allows us to compare the cumulative dis-tribution function for any 0 ¤ r ¤ 4:

3.3

3.5

3.555

0.2

0.4

3.585

3.6

3.69

3.74

3.84

3.96

1.0

0.8

0.6

0.4

0.2

0.0

0.6

0.8

1.0

We observe that for r = 2 the orbit stabilizes in a single point, so that there is just one possible event that happens with probability 1. For r = 3.3 there are two possible results that occur with probability 1 • 2. For there are four equiprobable results. For there are eight. The chaotic behavior can