# Replicator dynamics in discrete time
# RepDyn2.mws Ñ ESA, 12 Nov 98
# 1. Replicator dynamics
# 
# In formal models of economic evolution replicator dynamics plays an
# important role (presented by Hofbauer and Sigmund, Silverberg, etc.).
# This dynamics is most easy to understand in the case where we have a
# pure selection process, i.e. where a population with a given initial
# set of different routines undergoes a process of selection. 
# 
# If the population share of each of  n; competing "species" (firms,
# technologies) is s[i,t]; and their "competetiveness" or "relative
# fitness" is c[i,t];, then the so-called replicator equation or Fisher
# equation represents a simple selection mechanism. It can be
# represented in two ways. First, there is the differential equation:
> ds[i]/dt := alpha*(c[i]-C)*s[i];#      (1)
# Second, there is the difference equation:
> s[i,t+1] := s[i,t]+alpha*(c[i,t]-C[t])*s[i,t];#      (2)
# In both cases C[t] := sum(c[i,t]*s[i,t],i = 1 .. n);. This is called
# the population-weighted average fitness. Of course, we also have that
# S[t] := sum(s[i,t],i = 1 .. n); is equal to 1.
# 
# Fisher's (1930) theorem of natural sepection is that the average
# competitiveness (C[t];) increases at a rate proportional to the
# variance of the relative fitnesses (c[i,t];) of the population. As the
# population share of the most competitive "species" increases, the
# variance decreases and thus the rate of change. In the end the
# polulation has a complete dominance of one species, and there is no
# further change. Of course, things get more complicated if there is a
# process of mutation/innovation by which new "species" are introduced
# into the population.
# 
# It is a mechanism like the replicator dynamics which underlies e.g.
# Winter's (1964) pioneering critique of Alchian (1950) Friedman (1953).
# Winter agrees with Alchian and Friedman that firms with lower than
# average unit costs have higher profits to reinvest and thus they can
# expand capacity and enlarge their market shares. However, this process
# does not necessarily lead to overall use of best-practice technology.
# The reason is e.g. that an efficient selection process presuppose that
# there is some degree of "stickiness" of the productivity of the firm
# and that the alpha; parameter has a size that is high enough to remove
# the reintroduction of low-productivity behaviour in the population.
# 
# 2. Numerical simulations and their problems
# 
# Equation (1) is used for deriving the theoretical results of Fisher,
# but equation (2) is often used for evolutionary simulations. In this
# case we have problems because of the stepwise movement of the (market)
# shares. To understand and solve this problem, we of-course have to
# follow the KISS principle. At the same time we should recognise that
# we are facing a basic problem of numerical computation which probably
# has been studied extensively and for which we probably can find help
# in the Bible (i.e. in Numerical Recipies in C). 
# 
# For the moment we shall, however, appraoch the problem in three steps.
# First, we shall study pure selection. Second, we shall study the
# influence of simple random innovations. Third, we shall try to solve
# the problems we see.
# 
# 2.1. Pure selection
# 
# In this note our means of studying the replicator dynamics is Maple
# (V.5), a symbolic math package which can also be used for numerical
# expiriments. We could also have used Excel or the LSD system. In our
# experiments we shall use randomly distributed competitivenesses. They
# are produced by Maple's random number generator which in the standard
# setting gives a 12-digit random number, e.g.
# 
> _seed := 1; rand(); rand(); 

                              _seed := 1


                             427419669081


                             321110693270

# To simplify, we start by defining a Maple procedure that generates
# competitivenesses between 0 and 1:
> Rnd := proc()
>    evalf(rand()/10^12)
> end:
> _seed := 1; Rnd(); Rnd(); 

                              _seed := 1


                             .4274196691


                             .3211106933

# 
# Using this procedure, we write a simple Maple program (SimpleRepDyn) 
# that reflects the following situation. We have N; firms with fixed and
# randomly distributed competitivesses between 0 and 1. All firms have
# equal market share of size 1/N;. In each period market shares are
# updated using equation (2). 
# 
> SimpleRepDyn := proc(N,T,alpha,Dig)
> global c,C,s,S;
> local i,t;
> 
>    Digits := Dig;
> 
>    c := array(1..N); for i from 1 to N do c[i] := Rnd() od;
>    s := array(1..N,0..T); for i from 1 to N do s[i,0] := 1/N od;
>    S := array(0..T); S[0] := sum(s['j',0],'j'=1..N);
>    C := array(0..T): C[0] := sum(c['j']*s['j',0],'j'=1..N);
> 
>    for t from 1 to T do
>       for i from 1 to N do
>          s[i,t] := s[i,t-1] + alpha*(c[i]-C[t-1])*s[i,t-1]
>       od;
>       C[t] := sum(c['j']*s['j',t],'j'=1..N);
>       S[t] := sum(s['j',t],'j'=1..N)
>    od;
>    print(`OK`)
> end:
# 
# We now test the results of SimpleRepDyn:
# 
> SimpleRepDyn(5,300,0.3,10);
                                            OK
> plot(S[t],t=0..300);

> plot(C[t],t=0..300);

> plot([s[1,t],s[2,t],s[3,t],s[4,t],s[5,t]],t=1..300);

# 
# The calculations are influenced by the precision of our calculations.
# If we fix the number of digits to 2 rather than to 10, we get a more
# irregular behaviour:
# 
> SimpleRepDyn(5,300,0.3,3);
                                            OK
> plot(S[t],t=0..300);

> plot(C[t],t=0..300);

> plot([s[1,t],s[2,t],s[3,t],s[4,t],s[5,t]],t=1..300);

# 2.2. Introducing innovations
# 
# Fisher's law of simple selection processes is quite simple: in the
# long run the firm with the best competitiveness takes over the whole
# market. But what happens if innovations are introduced? This can be
# studied by including a simple change process whereby firms in each
# period get a random increase in their competitiveness. More
# specifically, the competitiveness of each firm is increased
# incrementally according to the the following equation:
> c[i,t] := c[i,t-1]+beta*Rnd;#      (3)
# where Rnd is a random number between 0 and 1.
# 
#  Let ud try to modify SimpleRepDyn to InnoRepDyn to reflect this
# process
# .
> InnoRepDyn := proc(N,T,alpha,beta,Dig)
> global c,C,s,S;
> local i,t;
> 
>    Digits := Dig;
> 
>    c := array(1..N,0..T); for i from 1 to N do c[i,0] := Rnd() od;
>    s := array(1..N,0..T); for i from 1 to N do s[i,0] := 1/N od;
>    S := array(0..T); S[0] := sum(s['j',0],'j'=1..N);
>    C := array(0..T): C[0] := sum(c['j',0]*s['j',0],'j'=1..N);
> 
>    for t from 1 to T do
>       for i from 1 to N do
>          s[i,t] := s[i,t-1] + alpha*(c[i,t-1]-C[t-1])*s[i,t-1];
>          c[i,t] := c[i,t-1] + beta*Rnd()
>       od;
>       C[t] := sum(c['j',t]*s['j',t],'j'=1..N);
>       S[t] := sum(s['j',t],'j'=1..N)
>    od;
>    print(`OK`)
> end:
> InnoRepDyn(5,30,0.3,0.3,10);
                                            OK
> plot(S[t],t=0..30);

> plot(C[t],t=0..30);

> plot([s[1,t],s[2,t],s[3,t],s[4,t],s[5,t]],t=1..30);

# This looks OK, but if we make further computations, we find out that
# something strange is happening. In a particular run we see that around
# the famous 68 generation something strange is happening: the sum of
# the market shares starts to occillate around 1 and finally it moves in
# the famous period 73 towards 0, i.e. the system collapses. 
# 
> InnoRepDyn(5,74,0.3,0.3,10);
                                            OK
> plot(S[t],t=0..74);

> plot(C[t],t=0..74);

> plot([s[1,t],s[2,t],s[3,t],s[4,t],s[5,t]],t=1..74);

# 
# Here is, from another run with another set of random numbers, another
# kind of strange behaviour:
# 
> plot([s[1,t],s[2,t],s[3,t],s[4,t],s[5,t]],t=1..74);

# 
# 2.3. Correcting for the "border" problems
# 
# It is obvious that the problem is that the discrete form of the
# replicator equation (2) gives problems when we introduce innovation in
# the form of equation (3). One solution is that we normalise so that
# there are never market shares below 0 or above 1 that enters into
# equation (2). This problem seems to be solved in the procedure
# NewInnoRepDyn.
# 
> NewInnoRepDyn := proc(N,T,alpha,beta,Dig)
>    # To be written
> end:
>