##############################################################
#  LNW model and plotting/saving of variables for Maple 6
#  file: "LNW.txt"
#  last update: 01/06/18 13:21:33 
#  author: Esben Sloth Andersen, esa@business.auc.dk
##############################################################

# The LNW version of the Nelson--Winter model.  Compare with 
# Andersen's paper for the DRUID NW2001 Conference and the 
# website <http://business.auc.dk/evolution/esapapers/esa01/MNW.html>.

# The whole of the present file should be pasted into an empty 
# Maple 6 worksheet (one small correction is needed if Maple 
# V.5 is used).  Run the program with the command 
# LNW(n,T,Variable,Seed,Device); Here n is number of firms; T 
# is number of periods; Variable of the variable to be 
# plotted; Seed is an integer that defines a particular series 
# of random numbers; Device is where to send the results.

# The program can be run in two major ways.  With the command 
# LNW(n,T,Variable,Seed,screen); the last argument (screen) 
# sends the output to the screen and thus allows us directly 
# to inspect the output and make experiments until we have a 
# final result.  When a final result has been obtained, we run 
# the program with the command 
# LNW(n,T,Variable,Seed,filename); The last argument is 
# <filename> and this forces the program to write the time 
# series of the variable to a file named filename.dat.  The 
# file is placed in a predifined directory.  To change 
# directory, correct the PathName in the beginning of the LNW 
# program.  If the PathName is the empty string (""), your 
# system will probably place the output in the Maple 
# directory.  The format of the file is designed for 
# processing within the LaTeX/MetaPost system (see LaTeX # 
# Graphics, pp.  72 ff.  and Andersen's MetaPost examples).

# For fixed the setup was LNW(4,1000,A,7,LNWfixed_A);
# alpha := 0.2; 
# A_init := 0.16;
# d := 1.0;
# Ltot := 100;
# r[0] := 0.0;
# rho :=  1;
# sigma := 0.01;

# For random the setup was LNW(4,1000,A,7,LNWrandom_A);
# alpha := 0.2; 
# A_init := 0.16;
# d := 1.0;
# Ltot := 100;
# r[0] := 0.05;
# rho :=  1;
# sigma := 0.01;

# For normal the setup was LNW(4,200,s,20,screen);
# > alpha := 0.2; 
# > A_init := 0.16;
# > d := 1.0;
# > Ltot := 100;
# > r[0] := 0.05;
# > rho :=  1;
# > sigma := 0.03;

# For split the setup was LNW(4,150,s,20,screen);
# > alpha := 0.2; 
# > A_init := 0.16;
# > d := 1.0;
# > Ltot := 100;
# > r[0] := 0.05;
# > rho :=  1;
# > sigma := 0.01;

# For normal 0.5 the setup was LNW(4,150,s,20,screen);
# > alpha := 0.2; 
# > A_init := 0.16;
# > d := 1.0;
# > Ltot := 100;
# > r[0] := 0.05;
# > rho :=  0.5;
# > sigma := 0.03;

LNW := proc(n,T,Variable,Seed,Device)

local f, filepath, i, j, st, t;

global _seed, A, A_init, A_max, A_mean, A_res, DrawVariable, 
L, Legend, Ltot, MergeLists, Method, P, PathName, Q, Qtot, 
S, alpha, d, lambda, pi, r, rho, s, sigma;

# The fact that variables and constants are "global" means 
# that for debugguing purposes they can be inspected 
# interactively after the simulation

option `ES Andersen, Aalborg University, 26 May 01`;

# FIRST INFORMATION TO PROGRAM USER
print(cat(`LNW is running with `,n,` firms for `,T,` periods.`));
print(cat(`Data will be sent to `,Device,`.`));
print(cat(`Pathname is: `,PathName));
st := time();

# PATHNAME FOR SAVING DATA FILES
PathName := `Macintosh HD:ENWandJEE:ENWpapertex:mpfigures:`;

# PARAMETERS
alpha := 0.2; 
A_init := 0.16;
d := 1.0;
Ltot := 100;
r[0] := 0.05;
rho :=  1;
sigma := 0.01;
_seed := Seed;
Method := random;
# If Method = random, then xxx.
# If Method = normal, then xxx.
# If Method = split, half of the firms have no R&D.
# If Method = fixed, there is no productivity change.

# MAIN VARIABLES
A := table();
A_mean := table();
L := table();
Q := table();
s := table();
pi := table();
        
# INITIALISATION
with(stats); 
with(statevalf);
for i from 1 to n do
  if Method = fixed then
    A[i,1] := exp(stats[random,normald[ln(A_init),sigma]]());
  else
    A[i,1] := A_init;
  fi;
  L[i,1] := Ltot/n;
  r[i] := r[0];
od;
if Method = split then
  for i from 1 to n/2 do
    r[i] := r[0];
    r[n/2+i] := 0;
  od;
fi;
A_mean[0] := 0.16;

# MAIN PROGRAM
for t from 1 to T do

   # SHORT RUN
   for i from 1 to n do 
      Q[i,t] := A[i,t]*(1-r[i])*L[i,t];
   od;  
   Qtot := sum(Q[k,t], k = 1..n);
   P := Ltot/Qtot;
   for i from 1 to n do 
      pi[i,t] := P*Q[i,t] - L[i,t];
      s[i,t] := Q[i,t]/Qtot;
      S[i,t] := L[i,t]/Ltot;
   od;

   # NEW TECHNO
   for i from 1 to n do 
     A_res[i] := 0;
   od;
   A_max := max(seq(A[i,t], i = 1..n));
   A_mean[t] :=  A_mean[t-1]*(1 + phi);
   for i from 1 to n do 
     if Method = fixed then
       A[i,t+1] := A[i,t];
     fi;
     if Method = random then
       A_res[i] := exp(stats[random,normald[ln(A[i,t]),sigma]]());
       if A_res[i] > A[i,t] then
         A[i,t+1] := A_res[i];
       else
         A[i,t+1] := A[i,t];
       fi;
     fi;
     if Method = normal or Method = split then
       # Research
       lambda[i] := d*r[i]*L[i,t];
       if stats[random,poisson[lambda[i]]]() > 0 then
         if evalf(rand()/1000000000000) < rho then
           A_res[i] := 
             exp(stats[random,normald[ln(A[i,t]),sigma]]()); 
         else A_res[i] := A_max; 
         fi;
       fi;

       # Technochoice
       A[i,t+1] := max(A[i,t], A_res[i]);
     fi;
   od;


   # NEW LABOUR
   for i from 1 to n do

      L[i,t+1] := L[i,t] + pi[i,t];

   od;

od;

# PROCEDURES FOR PLOTTING OF OUTPUT TO TERMINAL

# Function that helps the creation of a PLOT data structure

MergeLists := proc(list1, list2)
local LookUpItem,i;
   LookUpItem := proc(list1, list2, seriesnumber)
   local listnumber;
      if type(seriesnumber, even) then
         listnumber := seriesnumber/2;
         RETURN(list1[listnumber]);
      else
         listnumber := (seriesnumber-1)/2;
         RETURN(list2[listnumber]);
      fi;
    end;       
if not nops(list1) = nops(list2) then
    ERROR(`the number of elements in lists must be the same`);
fi;
[seq(LookUpItem(list1, list2, i), 
     i = 2..nops(list1)*2 + 1)];     
end;

# Legend for plots.  For Maple V.5 the concatenation operator 
# "||" must be replaced by a dot ".".

Legend := proc(N)
global periods,plotdata,text;
local colour,i,n,a,b;
colour := seq('COLOUR'(HUE,n/N), n=1..N);
a:=NULL; b := NULL;
for n from 1 to N do
    a := a, 'CURVES'([[0,N-n],[10,N-n]],
    colour[n],THICKNESS(2));
    b := b, 'TEXT'([11,N-n],`Firm `||n,FONT(TIMES,ROMAN,12));
od;
text := cat(`Legend for `,N,` firms`);
plotdata := 'PLOT'(a, b, 'TITLE'(text), AXESSTYLE(NONE));
plotdata;

end;

# Procedure that plots a time series for a firm-level 
# variable.  This procedure draws time series data for the 
# firms of the industry.  In order to allow simple comparisons 
# across time series for different variables, each firm is 
# given a specific colour that is the same in all plots.

DrawVariable := proc(var,N,T)
global periods,plotdata,text;
local colour,tt,t,n,a;
periods := seq(t, t = 1..T);
colour := seq('COLOUR'(HUE,n/N), n=1..N);
a:=NULL;
for n from 1 to N do
        a := a, 'CURVES'([seq(MergeLists([periods[tt]],[var[n,tt]]),tt=1..T)],
    colour[n],THICKNESS(2));
od;
text := cat(`Variable `,var,` for `, N, ` firms and `, T, ` periods.`);
plotdata := 'PLOT'(a, 'TITLE'(text), AXESSTYLE(NORMAL), AXESTICKS(5,5));
plotdata;
end;

# OUTPUT FROM THE LNW PROGRAM
for i from 1 to n do
for t from 1 to T do
 if s[i,t] < evalf(exp(-20)) then
   s[i,t] := s[i,t-1]
 fi;
od;
od;
print(cat(`Computation time: `,time() - st,` seconds.`));
print(``);
if Device = screen then
  print(DrawVariable(Variable,n,T));
else
  filepath := cat(PathName,Device,`.dat`);
  f := fopen(filepath,WRITE,TEXT);
  for i from 1 to n do
    if n*T < 502 then
      for t from 1 to T do
        fprintf(f, "%g %g\n",t,Variable[i,t]);
      od;
      fprintf(f, "\n");
    else
      for t from 1 by 5 to T do
        fprintf(f, "%g %g\n",t,Variable[i,t]);
      od;
      fprintf(f, "\n");
    fi;
  od;
  fclose(f);
fi;
end: