This is the program file for James B. Bullard and Steven Russell's
article "How Costly is Sustained Low Inflation for the U.S. Economy?" in 
the May/June 2004 issue of the Federal Reserve Bank of St. Louis Review.

This program creates the data in Tables 3, 4, and 7, as well as Figures 3 and 4.



(*  $Output = OpenAppend["d:\\programs\\welfare\\calga71.out",FormatType->OutputForm]; *)

Print["THIS IS THE GA PROGRAM TO CALIBRATE THE MODEL."];
Print["Choose nine deep parameters to hit nine targets."];
Print["How costly is continued low inflation for the US economy?  J. Bullard and S. Russell."];
Print["This version: August 3, 1998."];
Print[" "];

(* print detailed output if print is 1 *)

print=0;

(* ga parameters *)

numberofstrings = 30; (* number of strings must be an even integer *)
stringlength = 9; (* program is hardcoded for stringlength of 9 *)
pcross = .95;
pshuffle = .33;
parithmatic = .34;
psinglepoint = 1-pshuffle-parithmatic; (* pshuffle + parithmatic must be less than 1 *)
pmutate = 1/stringlength;
gamax = 1000;
bbb = 2;

(* min and max of parameters *)
(* this is the parameter space we are searching over *)

rhomin = 0.0; rhomax = 0.1;
gammamin = 1.1; gammamax = 40;
etamin = .075; etamax = .4;
alphamin = .25; alphamax = .4;
deltamin = .025; deltamax = .1;
costmin = .01; costmax = .04;
rrmin = .01; rrmax = .08;
tauwmin = .01; tauwmax = .4;
taukcmin = .01; taukcmax = .25;

(* Target values *)

targeticgmin = .01; targeticgmax = .03; targeticg = .015;
targetaltmin = .075; targetaltmax = .4; targetalt = .3156; (* Note--based on 14 hour day *)
targetboymin = .255; targetboymax = .686; targetboy = .47;
targetkoymin = 2.32; targetkoymax = 4.32; targetkoy = 3.32;
targetiokmin = .066; targetiokmax = .086; targetiok = .076;
targetcoymin = .05; targetcoymax = .07; targetcoy = .06;
targetmoymin = .041; targetmoymax = .078; targetmoy = .0592;
targetgoymin = .121; targetgoymax = .18; targetgoy = .151;
targettkcmin = .036; targettkcmax = .201; targettkc = .119;

(*  Parameter settings  *)
(* These are the fixed parameters in our setup *)

(*  Number of periods in agents' lives  *)
n = 55;  (*  Program is hard-coded for n=55  *)

(*  Preferences  *)
lbar = 1;
periodsofretirement = 11;

(*  Production, inputs and depreciation  *)
lambda = 1.015;
psi = 1.017;

(*  Government policy parameters  *)
pi = .04;
caprm = 1/(1+pi);
taukp = .2;

(*  Target caprb  *)
targetcaprb = 1.05;
Print[" "];
Print["Note: targetcaprb = ",targetcaprb];
Print[" "];
 
(*  Productivity profile  *)
(*  Set flag to choose endowment profile  *)
endow = 2;
j=.;
If[endow==0,
     Print["Using Auerbach-Kotlikoff."];
     pattern=Exp[4.47 + .033*(j-20) - .00067*(j-20)^2];
,
     If[endow==1,
          Print["Using Hansen 1."];
          pattern = Exp[-1.46+.070731*j-.00075*(j^2)];
     ,
          If[endow==2,
               Print["Using Hansen 2."];
               pattern=-4.34+.61312*j-.027397*(j^2)+.000630074*(j^3)-
                    .00000717078*(j^4)+.0000000314186*(j^5);
          ,
               Print["Warning:  Endowment choice failure."];
          ];
     ];
];
normalize = NIntegrate[pattern,{j,21,21+55/n}];
For[jj=1,jj<n+1,jj++,
     jmin=21 + (jj-1)*(55/n);
     jmax=jmin + 55/n;
     e[jj] = (NIntegrate[pattern,{j,jmin,jmax}]/normalize);
];

(*  End preliminary portion of program  *)
Print["End preliminary section of program."];
Print[" "];

(* Begin ga section of program. *)
Print["Begin ga section of program."];
Print[" "];

(* Create some strings *)

If[print==1,Print["Creating a set of strings."]];
For[i=1,i<=numberofstrings,i++,
     x[i,1] = Random[Real,{rhomin,rhomax}];
     x[i,2] = Random[Real,{gammamin,gammamax}];
     x[i,3] = Random[Real,{etamin,etamax}];
     x[i,4] = Random[Real,{alphamin,alphamax}];
     x[i,5] = Random[Real,{deltamin,deltamax}];
     x[i,6] = Random[Real,{costmin,costmax}];
     x[i,7] = Random[Real,{rrmin,rrmax}];
     x[i,8] = Random[Real,{tauwmin,tauwmax}];
     x[i,9] = Random[Real,{taukcmin,taukcmax}];
];

(* Begin ga loop *)

For[ga=1,ga<=gamax,ga++,
If[print==1,Print["ga=",ga]];

(* Evaluate fitness of existing strings via loop *)

For[floop=1,floop<=numberofstrings,floop++,
     
     If[print==1,Print["floop=",floop]];

     rho = x[floop,1];
     gamma = x[floop,2];
     eta = x[floop,3];
     alpha = x[floop,4];
     delta = x[floop,5];
     cost = x[floop,6];
     rr = x[floop,7];
     tauw = x[floop,8];
     taukc = x[floop,9];

     beta = 1/(1+rho);
     nu = 1-eta*(1-gamma);
     sigma = 1/nu;

     (*  Interest rate definition  *)
     (*  Begin by assuming agents will save in every period of life.  *)

     caprb = targetcaprb;
     caprd = caprb;
     caprdat = (1-taukp)*caprd + taukp*caprm;
     For[i=1,i<=n-1,i++,
          caprhat[i] = (1-taukp)*caprd + taukp*caprm;
     ];
     caprkcost = (caprb - rr*caprm)/(1-rr);
     caprkact = caprkcost + cost;
     caprknetd = (caprkact - taukc*caprm)/(1 - taukc);
     caprktotal = caprknetd + delta;

     k[1] = ((caprknetd-1+delta)/alpha)^(1/(alpha-1));
     w[1] = (1-alpha)*(k[1]^alpha);

     solution=0;
     borrow=0;
     retire=periodsofretirement;
     cornerflag=0;

     While[solution==0,

          Clear[c];
          Clear[l];

          (* Construct budget constraint equation as function of c1 alone *)

          bc = c1 + Sum[c[i]*Product[(1/caprhat[j]),{j,1,i-1}],{i,2,n}] -
               (1-tauw)*w[1]*e[1]*(lbar-l[1]) - 
               Sum[(lambda^(ii-1))*(1-tauw)*w[1]*e[ii]*(lbar-l[ii])*
               Product[1/caprhat[jj],{jj,1,ii-1}],{ii,2,n-retire}];
          l[1] = ((1-eta)/eta)*(c1/(w[1]*(1-tauw)*e[1]));
          For[i=2,i<=n-retire,i++,
               c[i] = c1*((((lambda^(i-1))*e[i])/(e[1]))^(((1-gamma)*(eta-1))/gamma))*
                    (beta^((i-1)/gamma))*Product[caprhat[j]^(1/gamma),{j,1,i-1}];
               l[i] = ((1-eta)/eta)*(c[i]/((lambda^(i-1))*w[1]*(1-tauw)*e[i]));
          ];
          If[retire >= 1,
               For[i=n-retire+1,i<=n,i++,
                    c[i] = (((c1^-gamma)*(((1-eta)/(eta*(1-tauw)*w[1]*e[1]))^((1-eta)*(1-gamma))))/
                         ((lbar^((1-eta)*(1-gamma)))*(beta^(i-1))*
                         Product[caprhat[j],{j,1,i-1}]))^(1/(eta-1-gamma*eta));
                    l[i] = lbar;
               ];
          ];

          (* Now solve for c1] *)

          c[1] = c1 /. FindRoot[bc==0,{c1,.15}];

          (* Find implied values for c[i], l[i], and s[i] *)

          l[1] = ((1-eta)/eta)*(c[1]/((1-tauw)*w[1]*e[1]));
          assets[1] = (1-tauw)*w[1]*e[1]*(lbar- l[1]);
          s[1]=assets[1]-c[1];
          u[1] = (1/(1-gamma))*(((c[1]^eta)*(l[1]^(1-eta)))^(1-gamma));
          For[i=2,i<=n-retire,i++,
               c[i] = c[1]*((((lambda^(i-1))*e[i])/e[1])^(((1-gamma)*(eta-1))/gamma))*
                         (beta^((i-1)/gamma))*Product[caprhat[j]^(1/gamma),{j,1,i-1}];
               l[i]=((1-eta)/eta)*(c[i]/((lambda^(i-1))*(1-tauw)*w[1]*e[i]));
               assets[i] = ((1-tauw)*(lambda^(i-1))*w[1]*e[i]*(lbar-l[i]))+(caprhat[i-1]*s[i-1]);
               s[i]=assets[i]-c[i];
               u[i] = (beta^i)*(1/(1-gamma))*(((c[i]^eta)*(l[i]^(1-eta)))^(1-gamma));
          ];
          If[retire>=1,
               For[i=n-retire+1,i<=n,i++,
                    c[i] = (((c[1]^-gamma)*(((1-eta)/(eta*(1-tauw)*w[1]*e[1]))^((1-eta)*
                         (1-gamma))))/((lbar^((1-eta)*(1-gamma)))*
                         (beta^(i-1))*(Product[caprhat[j],{j,1,i-1}])))^(1/(eta-1-gamma*eta));
                    l[i] = lbar;
                    assets[i] = caprhat[i-1]*s[i-1];
                    s[i]=assets[i] - c[i];
                    u[i] = (beta^i)*(1/(1-gamma))*(((c[i]^eta)*(l[i]^(1-eta)))^(1-gamma));
               ];
          ];


          If[cornerflag==0 && borrow>=1 && s[borrow] > 0,
               If[print==1,Print["Corner solution detected, approximation proceeding."]];
               caprhat[borrow] = caprdat +
                    (-borrowatcaprdat/(-borrowatcaprdat+s[borrow]))*(caprkact - caprdat);
               If[print==1,Print["borrowatcaprdat = ",borrowatcaprdat]];
               If[print==1,Print["saveatcaprkact = ",s[borrow]]];
               If[print==1,Print["caprhat[",borrow,"] has been set to ",caprhat[borrow]," ..."]];
               cornerflag = 1;
          ,
               (*  Now check for retirement and borrowing, with no corner problem. *)
               If[cornerflag==1,
                    If[print==1,Print["... and savings at the corner turns out to be ",s[borrow]]];
               ];
               If[l[n-retire]>lbar,
                    If[s[borrow+1]<0,
                         Print["Warning: both retirement and borrowing need correction."];
                    ];
                    (* Correcting for retirement, assuming any retirement occurs late in life *)
                    If[print==1,Print["Adjusting for retirement in period ",n-retire]];
                    retire=retire+1;
               ,
                    (* Check and correct for borrowing *)
                    (*  Assuming any borrowing occurs early in life  *)
                    If[s[borrow+1]<0,
                         borrow=borrow+1;
                         cornerflag=0;
                         borrowatcaprdat = s[borrow];
                         If[print==1,Print["Adjusting for borrowing in period ",borrow]];
                         caprhat[borrow] = caprkact;
                    ,
                         For[i=1,i<=n-retire,i++,
                              If[l[i]>lbar,
                                   Print["Warning: additional retirement occurring."]; 
                              ];
                         ];
                         For[i=borrow+1,i<=n-1,i++,
                              If[s[i]<0,
                                   Print["Warning: additional borrowing occurring."];
                              ];
                         ];
                         solution = 1;

                    ]; (* end if on borrowing *)

               ]; (* end if on retirement *)

          ]; (* end if checking for corner conditions *)

     ]; (* End while no solution loop *)

     (*  We now have c[i], l[i], s[i] for this value of caprk *)
     (*  Use these results to define aggregate savings and other aggregate variables  *)

     aggs[1] = s[1];
     For[i=2,i<=n,i++,
          aggs[i] = ((lambda*psi)^(1-i))*s[i];
     ];
     caps = Sum[aggs[i],{i,1,n-1}];
     capsplus = 0;
     capsminus = 0;
     For[i=1,i<=n-1,i++,
          If[aggs[i]<0,
               capsminus = capsminus - aggs[i];
          ,
               capsplus = capsplus + aggs[i];
          ];
     ];

     (*  Aggregate values at this value of caprk  *)

     capl = Sum[(psi^(1-i))*e[i]*(lbar-l[i]),{i,1,n}];
     capkt = k[1]*capl;
     capktp1 = lambda*psi*capkt;
     capm = (rr/(1-rr))*(capsminus+capktp1);
     capb = capsplus - (1/(1-rr))*(capsminus + capktp1);
     capy = (capl^(1-alpha))*(capkt^alpha);
     capg = tauw*w[1]*capl + taukp*(caprd-caprm)*(capsplus/(lambda*psi)) + 
          (1-(caprm/(lambda*psi)))*capm + (1-(caprb/(lambda*psi)))*capb +
          taukc*(caprknetd-caprm)*capkt;
     indcongrth = (1/(n-1))*Sum[1+((c[i]-c[i-1])/c[i-1]),{i,2,n}];
     avglabortime=(1/n)*Sum[(lbar-l[i])/lbar,{i,1,n}];

     (* Evaluate this steady state for fitness *)

     eqicg = indcongrth -1;
     If[eqicg<targeticgmin,
          icgfit[floop] = (Abs[(eqicg - targeticg)/targeticg] + 
               ((eqicg - targeticgmin)/targeticgmin)^2)/
               (Abs[(targeticgmin - targeticg)/targeticg]);
     ,
          If[eqicg>targeticgmax,
               icgfit[floop] = (Abs[(eqicg - targeticg)/targeticg] +
                    ((eqicg - targeticgmax)/targeticgmax)^2)/
                    (Abs[(targeticgmax - targeticg)/targeticg]);
          ,
               If[eqicg>targeticg,
                    icgfit[floop] = Abs[(eqicg - targeticg)/targeticgmax]/
                         Abs[(targeticgmax - targeticg)/targeticgmax];
               ,
                    icgfit[floop] = Abs[(eqicg - targeticg)/targeticgmin]/
                         Abs[(targeticgmin - targeticg)/targeticgmin];
               ];
          ];
     ];
     eqalt = avglabortime;
     If[eqalt<targetaltmin,
          altfit[floop] = (Abs[(eqalt - targetalt)/targetalt] + 
               ((eqalt - targetaltmin)/targetaltmin)^2)/
               (Abs[(targetaltmin - targetalt)/targetalt]);
     ,
          If[eqalt>targetaltmax,
               altfit[floop] = (Abs[(eqalt - targetalt)/targetalt] +
                    ((eqalt - targetaltmax)/targetaltmax)^2)/
                    (Abs[(targetaltmax - targetalt)/targetalt]);
          ,
               If[eqalt>targetalt,
                    altfit[floop] = Abs[(eqalt - targetalt)/targetaltmax]/
                         Abs[(targetaltmax - targetalt)/targetaltmax];
               ,
                    altfit[floop] = Abs[(eqalt - targetalt)/targetaltmin]/
                         Abs[(targetaltmin - targetalt)/targetaltmin];
               ];
          ];
     ];
     eqboy = capb/capy;
     If[eqboy<targetboymin,
          boyfit[floop] = (Abs[(eqboy - targetboy)/targetboy] + 
               ((eqboy - targetboymin)/targetboymin)^2)/
               (Abs[(targetboymin - targetboy)/targetboy]);
     ,
          If[eqboy>targetboymax,
               boyfit[floop] = (Abs[(eqboy - targetboy)/targetboy] +
                    ((eqboy - targetboymax)/targetboymax)^2)/
                    (Abs[(targetboymax - targetboy)/targetboy]);
          ,
               If[eqboy>targetboy,
                    boyfit[floop] = Abs[(eqboy - targetboy)/targetboymax]/
                         Abs[(targetboymax - targetboy)/targetboymax];
               ,
                    boyfit[floop] = Abs[(eqboy - targetboy)/targetboymin]/
                         Abs[(targetboymin - targetboy)/targetboymin];
               ];
          ];
     ];
     eqkoy = capkt/capy;
     If[eqkoy<targetkoymin,
          koyfit[floop] = (Abs[(eqkoy - targetkoy)/targetkoy] + 
               ((eqkoy - targetkoymin)/targetkoymin)^2)/
               (Abs[(targetkoymin - targetkoy)/targetkoy]);
     ,
          If[eqkoy>targetkoymax,
               koyfit[floop] = (Abs[(eqkoy - targetkoy)/targetkoy] +
                    ((eqkoy - targetkoymax)/targetkoymax)^2)/
                    (Abs[(targetkoymax - targetkoy)/targetkoy]);
          ,
               If[eqkoy>targetkoy,
                    koyfit[floop] = Abs[(eqkoy - targetkoy)/targetkoymax]/
                         Abs[(targetkoymax - targetkoy)/targetkoymax];
               ,
                    koyfit[floop] = Abs[(eqkoy - targetkoy)/targetkoymin]/
                         Abs[(targetkoymin - targetkoy)/targetkoymin];
               ];
          ];
     ];
     eqiok = (capktp1 - capkt + delta*capkt)/capkt;
     If[eqiok<targetiokmin,
          iokfit[floop] = (Abs[(eqiok - targetiok)/targetiok] + 
               ((eqiok - targetiokmin)/targetiokmin)^2)/
               (Abs[(targetiokmin - targetiok)/targetiok]);
     ,
          If[eqiok>targetiokmax,
               iokfit[floop] = (Abs[(eqiok - targetiok)/targetiok] +
                    ((eqiok - targetiokmax)/targetiokmax)^2)/
                    (Abs[(targetiokmax - targetiok)/targetiok]);
          ,
               If[eqiok>targetiok,
                    iokfit[floop] = Abs[(eqiok - targetiok)/targetiokmax]/
                         Abs[(targetiokmax - targetiok)/targetiokmax];
               ,
                    iokfit[floop] = Abs[(eqiok - targetiok)/targetiokmin]/
                         Abs[(targetiokmin - targetiok)/targetiokmin];
               ];
          ];
     ];
     eqcoy  = (cost*(capkt + capsminus))/capy;
     If[eqcoy<targetcoymin,
          coyfit[floop] = (Abs[(eqcoy - targetcoy)/targetcoy] + 
               ((eqcoy - targetcoymin)/targetcoymin)^2)/
               (Abs[(targetcoymin - targetcoy)/targetcoy]);
     ,
          If[eqcoy>targetcoymax,
               coyfit[floop] = (Abs[(eqcoy - targetcoy)/targetcoy] +
                    ((eqcoy - targetcoymax)/targetcoymax)^2)/
                    (Abs[(targetcoymax - targetcoy)/targetcoy]);
          ,
               If[eqcoy>targetcoy,
                    coyfit[floop] = Abs[(eqcoy - targetcoy)/targetcoymax]/
                         Abs[(targetcoymax - targetcoy)/targetcoymax];
               ,
                    coyfit[floop] = Abs[(eqcoy - targetcoy)/targetcoymin]/
                         Abs[(targetcoymin - targetcoy)/targetcoymin];
               ];
          ];
     ];
     eqmoy = capm/capy;
     If[eqmoy<targetmoymin,
          moyfit[floop] = (Abs[(eqmoy - targetmoy)/targetmoy] + 
               ((eqmoy - targetmoymin)/targetmoymin)^2)/
               (Abs[(targetmoymin - targetmoy)/targetmoy]);
     ,
          If[eqmoy>targetmoymax,
               moyfit[floop] = (Abs[(eqmoy - targetmoy)/targetmoy] +
                    ((eqmoy - targetmoymax)/targetmoymax)^2)/
                    (Abs[(targetmoymax - targetmoy)/targetmoy]);
          ,
               If[eqmoy>targetmoy,
                    moyfit[floop] = Abs[(eqmoy - targetmoy)/targetmoymax]/
                         Abs[(targetmoymax - targetmoy)/targetmoymax];
               ,
                    moyfit[floop] = Abs[(eqmoy - targetmoy)/targetmoymin]/
                         Abs[(targetmoymin - targetmoy)/targetmoymin];
               ];
          ];
     ];
     eqgoy = capg/capy;
     If[eqgoy<targetgoymin,
          goyfit[floop] = (Abs[(eqgoy - targetgoy)/targetgoy] + 
               ((eqgoy - targetgoymin)/targetgoymin)^2)/
               (Abs[(targetgoymin - targetgoy)/targetgoy]);
     ,
          If[eqgoy>targetgoymax,
               goyfit[floop] = (Abs[(eqgoy - targetgoy)/targetgoy] +
                    ((eqgoy - targetgoymax)/targetgoymax)^2)/
                    (Abs[(targetgoymax - targetgoy)/targetgoy]);
          ,
               If[eqgoy>targetgoy,
                    goyfit[floop] = Abs[(eqgoy - targetgoy)/targetgoymax]/
                         Abs[(targetgoymax - targetgoy)/targetgoymax];
               ,
                    goyfit[floop] = Abs[(eqgoy - targetgoy)/targetgoymin]/
                         Abs[(targetgoymin - targetgoy)/targetgoymin];
               ];
          ];
     ];
     eqtkc =  (taukc*(caprknetd-caprm)*capkt)/capg;
     If[eqtkc<targettkcmin,
          tkcfit[floop] = (Abs[(eqtkc - targettkc)/targettkc] + 
               ((eqtkc - targettkcmin)/targettkcmin)^2)/
               (Abs[(targettkcmin - targettkc)/targettkc]);
     ,
          If[eqtkc>targettkcmax,
               tkcfit[floop] = (Abs[(eqtkc - targettkc)/targettkc] +
                    ((eqtkc - targettkcmax)/targettkcmax)^2)/
                    (Abs[(targettkcmax - targettkc)/targettkc]);
          ,
               If[eqtkc>targettkc,
                    tkcfit[floop] = Abs[(eqtkc - targettkc)/targettkcmax]/
                         Abs[(targettkcmax - targettkc)/targettkcmax];
               ,
                    tkcfit[floop] = Abs[(eqtkc - targettkc)/targettkcmin]/
                         Abs[(targettkcmin - targettkc)/targettkcmin];
               ];
          ];
     ];
     fitness[floop] = icgfit[floop] + altfit[floop] + boyfit[floop] + koyfit[floop] + iokfit[floop] +
          coyfit[floop] + moyfit[floop] + goyfit[floop] + tkcfit[floop];

(* End fitness loop, the output is one fitness value for each string *)

];

(* Print a status line about average values *)

avgrho = Sum[(1/numberofstrings)*x[i,1],{i,1,numberofstrings}];
avggamma = Sum[(1/numberofstrings)*x[i,2],{i,1,numberofstrings}];
avgeta = Sum[(1/numberofstrings)*x[i,3],{i,1,numberofstrings}];
avgalpha = Sum[(1/numberofstrings)*x[i,4],{i,1,numberofstrings}];
avgdelta = Sum[(1/numberofstrings)*x[i,5],{i,1,numberofstrings}];
avgcost = Sum[(1/numberofstrings)*x[i,6],{i,1,numberofstrings}];
avgrr = Sum[(1/numberofstrings)*x[i,7],{i,1,numberofstrings}];
avgtauw = Sum[(1/numberofstrings)*x[i,8],{i,1,numberofstrings}];
avgtaukc = Sum[(1/numberofstrings)*x[i,9],{i,1,numberofstrings}];
Print[" "];
Print["Iteration ",ga];
Print["Average fitness is ",
     AccountingForm[Sum[(1/numberofstrings)*(fitness[i]),{i,1,numberofstrings}],5],
     " with average parameter values:"];
Print["     rh=",AccountingForm[avgrho,3],
     " gm=",AccountingForm[avggamma,3],
     " et=",AccountingForm[avgeta,3],
     " al=",AccountingForm[avgalpha,3],
     " dl=",AccountingForm[avgdelta,3]];
Print["     c=",AccountingForm[avgcost,3],
     " rr=",AccountingForm[avgrr,3],
     " tw=",AccountingForm[avgtauw,3],
     " tc=",AccountingForm[avgtaukc,3]];
avgicgfit = Sum[(1/numberofstrings)*icgfit[i],{i,1,numberofstrings}];
avgaltfit = Sum[(1/numberofstrings)*altfit[i],{i,1,numberofstrings}];
avgboyfit = Sum[(1/numberofstrings)*boyfit[i],{i,1,numberofstrings}];
avgkoyfit = Sum[(1/numberofstrings)*koyfit[i],{i,1,numberofstrings}];
avgiokfit = Sum[(1/numberofstrings)*iokfit[i],{i,1,numberofstrings}];
avgcoyfit = Sum[(1/numberofstrings)*coyfit[i],{i,1,numberofstrings}];
avgmoyfit = Sum[(1/numberofstrings)*moyfit[i],{i,1,numberofstrings}];
avggoyfit = Sum[(1/numberofstrings)*goyfit[i],{i,1,numberofstrings}];
avgtkcfit = Sum[(1/numberofstrings)*tkcfit[i],{i,1,numberofstrings}];
Print["Average fitness values:"];
Print["     icg=",AccountingForm[avgicgfit,3],
     " alt=",AccountingForm[avgaltfit,3],
     " boy=",AccountingForm[avgboyfit,3],
     " koy=",AccountingForm[avgkoyfit,3],
     " iok=",AccountingForm[avgiokfit,3]];
Print["     coy=",AccountingForm[avgcoyfit,3],
     " moy=",AccountingForm[avgmoyfit,3],
     " goy=",AccountingForm[avggoyfit,3],
     " tkc=",AccountingForm[avgtkcfit,3]];

(* Report the string with the highest fitness *)
(* Also keeps track of second best string  *)
(*
nfit1 = 1;
nfit2 = 1;
bestfit = 99999999;
For[i=1,i<=numberofstrings,i++,
     If[fitness[i]<bestfit, (* means string is better *)
          nfit2 = nfit1;
          nfit1 = i;
          bestfit = fitness[i];
     ];
];
Print["The best string is number ",nfit1,
     ", with fitness ",AccountingForm[fitness[nfit1],5],", and parameter values:"];
Print["     rh=",AccountingForm[x[nfit1,1],3],
     " gm=",AccountingForm[x[nfit1,2],3],
     " et=",AccountingForm[x[nfit1,3],3],
     " al=",AccountingForm[x[nfit1,4],3],
     " dl=",AccountingForm[x[nfit1,5],3]];
Print["     c=",AccountingForm[x[nfit1,6],3],
     " rr=",AccountingForm[x[nfit1,7],3],
     " tw=",AccountingForm[x[nfit1,8],3],
     " tc=",AccountingForm[x[nfit1,9],3]];
Print["Fitness values by target:"];
Print["     icg=",AccountingForm[icgfit[nfit1],3],
     " alt=",AccountingForm[altfit[nfit1],3],
     " boy=",AccountingForm[boyfit[nfit1],3],
     " koy=",AccountingForm[koyfit[nfit1],3],
     " iok=",AccountingForm[iokfit[nfit1],3]];
Print["     coy=",AccountingForm[coyfit[nfit1],3],
     " moy=",AccountingForm[moyfit[nfit1],3],
     " goy=",AccountingForm[goyfit[nfit1],3],
     " tkc=",AccountingForm[tkcfit[nfit1],3]];
Print[" "];
*)

(* Create a new generation of strings *)

(* Reproduction *)

If[print==1,Print["Reproduction."]];

(* Now use tournament selection *)

For[i=1,i<=numberofstrings,i++,
     string1 = Random[Integer,{1,numberofstrings}];
     string2 = Random[Integer,{1,numberofstrings}];
     If[fitness[string1] < fitness[string2], (*means string1 is better *)
          For[j=1,j<=stringlength,j++,
               nb[i,j] = x[string1,j];
          ];
     ,
          For[j=1,j<=stringlength,j++,
               nb[i,j] = x[string2,j];
          ];
     ];
];

(* Crossover *)

If[print==1,Print["Crossover."]];

(* no crossover for last two strings *)

For[i=1,i<=numberofstrings,i+=2,
          If[Random[]<pcross,

               typecross = Random[];
               
               (* Three types of crossover with given probability *)

               If[typecross < pshuffle,
                    (* Shuffle crossover *)
                    For[j=1,j<=stringlength,j++,
                         If[Random[]<.5,
                              x[i,j] = nb[i+1,j];
                              x[i+1,j] = nb[i,j];
                         ,
                              x[i,j] = nb[i,j];
                              x[i+1,j] = nb[i+1,j];
                         ];
                    ];
               ,
                    If[typecross < pshuffle+parithmatic,
                         (* Arithmatic crossover *)
                         aaa = Random[];
                         For[j=1,j<=stringlength,j++,
                              x[i,j] = aaa*nb[i,j] + (1-aaa)*nb[i+1,j];
                              x[i+1,j] = (1-aaa)*nb[i,j] + aaa*nb[i+1,j];
                         ];
                    ,
                         (* Single point crossover *)
                         crosspoint = Random[Integer,{1,stringlength}];
                         For[j=1,j<=stringlength,j++,
                              If[j<=crosspoint,
                                   x[i,j] = nb[i,j];
                                   x[i+1,j] = nb[i+1,j];
                              ,
                                   x[i,j] = nb[i+1,j];
                                   x[i+1,j] = nb[i,j];
                              ];
                         ];
                    ];
               ];
          ,
               (* No crossover at all *)
               For[j=1,j<=stringlength,j++,
                    x[i,j] = nb[i,j];
                    x[i+1,j] = nb[i+1,j];
               ];

          ];
];
     
(* Mutation *)

If[print==1,Print["Mutation."]];

For[i=1,i<=numberofstrings,i++,
     If[Random[]<pmutate,
          temp = x[i,1];
          rrr = Random[];
          If[Random[]<.5,
               x[i,1] = temp + (rhomax - temp)*(1 - rrr^((1-(ga/gamax))^bbb));
          ,
               x[i,1] = temp - (temp - rhomin)*(1 - rrr^((1-(ga/gamax))^bbb));
          ];
     ];
     If[Random[]<pmutate,
          temp = x[i,2];
          rrr = Random[];
          If[Random[]<.5,
               x[i,2] = temp + (gammamax - temp)*(1 - rrr^((1-(ga/gamax))^bbb));
          ,
               x[i,2] = temp - (temp - gammamin)*(1 - rrr^((1-(ga/gamax))^bbb));
          ];
     ];
     If[Random[]<pmutate,
          temp = x[i,3];
          rrr = Random[];
          If[Random[]<.5,
               x[i,3] = temp + (etamax - temp)*(1 - rrr^((1-(ga/gamax))^bbb));
          ,
               x[i,3] = temp - (temp - etamin)*(1 - rrr^((1-(ga/gamax))^bbb));
          ];
     ];
     If[Random[]<pmutate,
          temp = x[i,4];
          rrr = Random[];
          If[Random[]<.5,
               x[i,4] = temp + (alphamax - temp)*(1 - rrr^((1-(ga/gamax))^bbb));
          ,
               x[i,4] = temp - (temp - alphamin)*(1 - rrr^((1-(ga/gamax))^bbb));
          ];
     ];
     If[Random[]<pmutate,
          temp = x[i,5];
          rrr = Random[];
          If[Random[]<.5,
               x[i,5] = temp + (deltamax - temp)*(1 - rrr^((1-(ga/gamax))^bbb));
          ,
               x[i,5] = temp - (temp - deltamin)*(1 - rrr^((1-(ga/gamax))^bbb));
          ];
     ];
     If[Random[]<pmutate,
          temp = x[i,6];
          rrr = Random[];
          If[Random[]<.5,
               x[i,6] = temp + (costmax - temp)*(1 - rrr^((1-(ga/gamax))^bbb));
          ,
               x[i,6] = temp - (temp - costmin)*(1 - rrr^((1-(ga/gamax))^bbb));
          ];
     ];
     If[Random[]<pmutate,
          temp = x[i,7];
          rrr = Random[];
          If[Random[]<.5,
               x[i,7] = temp + (rrmax - temp)*(1 - rrr^((1-(ga/gamax))^bbb));
          ,
               x[i,7] = temp - (temp - rrmin)*(1 - rrr^((1-(ga/gamax))^bbb));
          ];
     ];
     If[Random[]<pmutate,
          temp = x[i,8];
          rrr = Random[];
          If[Random[]<.5,
               x[i,8] = temp + (tauwmax - temp)*(1 - rrr^((1-(ga/gamax))^bbb));
          ,
               x[i,8] = temp - (temp - tauwmin)*(1 - rrr^((1-(ga/gamax))^bbb));
          ];
     ];
     If[Random[]<pmutate,
          temp = x[i,9];
          rrr = Random[];
          If[Random[]<.5,
               x[i,9] = temp + (taukcmax - temp)*(1 - rrr^((1-(ga/gamax))^bbb));
          ,
               x[i,9] = temp - (temp - taukcmin)*(1 - rrr^((1-(ga/gamax))^bbb));
          ];
     ];
];

(* Do next iteration of ga loop *)

];


(* End program *)

Print[" "];
Print["End of program."];