#include "fun_head.h"

/***********************
rp2004 - Extendible replicator dynamics model
Esben Sloth Andersen, esa@business.aau.dk
Version: 9 April 2004

This file contains the Lsd macro code for exploring models that
are based on replicator dynamics. It emphasises structured
model development and the use of evolutionary statistics for
analysing model behaviour.
************************/

MODELBEGIN

// BASIC MODEL

EQUATION("z")
/*
Calculation of fitness characteristics
If RandWalk = 0, then fixed characteristics
If RandWalk = 1, then random walk of characteristics
*/
v[0]=V("RandWalk");
v[1]=V("InnoSize");
if (v[0]==0) v[2]=VL("z",1);
if (v[0]==1) v[2]=VL("z",1)+UNIFORM(-0.2,0.2);
RESULT(v[2])

EQUATION("N")
/*
Total number of members of the population
*/
v[0]=0;
CYCLE(cur, "Species")
 {
  v[0]=v[0]+VLS(cur,"n",1);
 }
RESULT(v[0])

EQUATION("Z")
/*
Mean characteristic
*/
v[0]=0;
v[1]=V("N");
CYCLE(cur, "Species")
 {
  v[0]=v[0]+VLS(cur,"n",1)*VS(cur,"z");
 }
RESULT(v[0]/v[1])

EQUATION("n")
/*
Replicator dynamics:
N[t]=N[t-1](1+a(f-af)/af)
*/
v[0]=VL("n",1);
v[1]=V("a");
v[2]=V("z");
v[3]=V("Z");
v[4]=v[0]*(1+v[1]*(v[2]-v[3])/v[3]);
RESULT(v[4])

// ADDITIONAL STATISTICS

EQUATION("s")
/*
Population share
*/
v[0]=V("n");
v[1]=V("N");
RESULT(v[0]/v[1])

EQUATION("InvHerf")
/*
Inverse Herfindahl index, 1/SUM(ms^2) 
   The Herfindahl is an indicator of the concentration of the market, 
computed as the square of population shares. It varies from 1 (one s=1 
and all the other equal 0) to 1/M (all species with identical shares).
   Its inverse corresponds to the number of species of the same 
size that would provide the same concentration index as the one 
measured. It varies between 1 (total concentration) and M (minimal 
concentration).
*/
v[0]=0;
CYCLE(cur,"Species")
 {
  v[1]=VS(cur,"s");
  v[0]=v[0]+v[1]*v[1];
 }
RESULT(1/v[0])

// FISSION AND FUSION OF SPECIES

EQUATION("SplitMe")
/*
If EntryExit = 0, then no fissions.
If EntryExit = 1, then fissions. Fissions of firms are
determined after the new species size is found. The 
probability of a fission depends on the parameter SplitChance 
times the population share (= lambda of a Poisson process). 
The actual splitting of a species is carried out by the
equation for Fissions.
*/
V("w"); // Ensure that reproduction coefficient is calculated
v[0] = V("EntryExit");
v[1] = V("SplitChance");
v[2] = V("s");
if (v[0]==1 && poisson(v[1]*v[2])>0)
  v[3]=1;
else
  v[3]=0;  
RESULT(v[3])

EQUATION("Fissions")
/*
This equation checks whether one or more species have set 
SplitMe = 1. If this is the case, the members of this/these 
population(s) are split according to the parameter MotherShare.
The new firm, is the smallest, and its fitness characteristic 
is equal to that of the mother firm.
*/
v[1]=V("MotherShare");
v[2]=0;
CYCLE(cur, "Species")
 {
  v[0]=VS(cur,"SplitMe");
  if(v[0]==1)
   {
   cur1=ADDOBJ_EX("Species",cur);
   MULTS(cur1,"n",(1-v[1]));
   MULTS(cur1,"s",(1-v[1]));
   // MULTS(cur1,"z",1.01); // The small is stronger!
   MULTS(cur,"n",v[1]);
   MULTS(cur,"s",v[1]);
   WRITELS(cur1,"SplitMe",0, t);
   }     
 v[2]++;
 }
RESULT(v[2])

EQUATION("Fusions")
/*
If EntryExit = 0, then no fusions.
If EntryExit = 1, then fusions. Here species whose population shares
are less than FusionLevel are deleted and their members are added
to the largest species.
*/
v[2]=v[5]=v[6]=v[7]=0;
v[4]=V("FusionLevel");
v[10] = V("EntryExit");
if (v[10]==1)
 {
  CYCLE_SAFE(cur, "Species")
   {
    v[0]=VS(cur,"s");
    if(v[0]<v[4])
     {
      v[2]+=VS(cur,"n"); 
      v[5]+=VS(cur,"s"); 
      DELETE(cur);
     } 
    if(v[6]<v[0])
     {
      v[6]=v[0];
      cur2=cur;
     } 
   }
  INCRS(cur2,"n",v[2]);
  INCRS(cur2,"s",v[5]);
 }
RESULT(v[2])

// PRICE'S STATISTICS
//         Task: Develop EVOSTAT!

EQUATION("Zchange")
/*
Change of mean characteristic
*/
v[0]=V("Z");
v[1]=VL("Z",1);
RESULT(v[0]-v[1])

EQUATION("w")
/*
w = n[t]/n[t-1]
   The reproduction coefficient of the species
*/
RESULT(V("n")/VL("n",1))

EQUATION("W")
/*
Weighted mean of the species' reproduction coefficients
*/
v[0]=0;
CYCLE(cur,"Species")
 {
  v[1] = VLS(cur,"s",1);
  v[2] = VS(cur,"w");
  v[0] = v[0]+v[1]*v[2];
 }
RESULT(v[0])

EQUATION("Var")
/*
The variance Var(z) = SUM[s[t-1]*(z[t-1]-Z[t-1])^2]
*/
v[0]=0;
v[1] = V("Z");
CYCLE(cur,"Species")
 {
  v[2] = VLS(cur,"s",1);
  v[3] = VLS(cur,"z",1);
  v[0] = v[0] + v[2]*(v[3]-v[1])*(v[3]-v[1]);
 }
RESULT(v[0])

EQUATION("SelectEffect")
/*
The selection effect is Cov(w,z)/W, where 
Cov(w,z) = SUM[s[t-1]*(w[t]-W[t])*(z[t-1]-Z[t-1]))]
   Covariance between species' reproduction coefficients and
their fitnesses
*/
v[0]=0;
v[1] = V("W");
v[2] = V("Z");
CYCLE(cur,"Species")
 {
  v[3] = VLS(cur,"s",1);
  v[4] = VS(cur,"w");
  v[5] = VLS(cur,"z",1);
  v[0] = v[0] + v[3]*(v[4]-v[1])*(v[5]-v[2]);
 }
RESULT(v[0]/v[1])

EQUATION("InnoEffect")
/*
E(s[t]*w[t]*(z[t]-z[t-1]))/W
The innovation effect as defined by George Price's equation.
*/
v[0]=0;
v[1] = V("W");
CYCLE(cur, "Species")
 {
  v[2] = VS(cur,"s");
  v[3] = V("w");
  v[4] = VS(cur,"z");
  v[5] = VLS(cur,"z",1);
  v[0] = v[0] + v[2]*v[3]*(v[4]-v[5]);
 }
RESULT(v[0]/v[1])


MODELEND

void close_sim(void)
{

}