Inequality and growth:
Why fertility differential matters
David de la Croix - Matthias Doepke
American Economic Review, 93, 1091-1113, 2003

Dynamic analysis

Procedures

SetOptions[ListPlot,PlotRange->All,Joined->True];

Defining the variables

defvar:=Module[{a},
Clear[cN,cNf,cn,cx,cxf,ce,pN,pn,px,pe,growth,cap,popgrowth,income,ginix,fd,cfd,birth];
Clear[Nh,Nhf,n,x,xf,e,plN,pln,plx,ple,gag,k,nag,inc,gix,ct,pct,cbr];
cN=Array[Nh,dim];cNf=Array[Nhf,dim];
cn=Array[n,dim];cx=Array[x,dim];cxf=Array[xf,dim]; ce=Array[e,dim];
income=Array[inc,dim];
pN=Array[plN,steps];pn=Array[pln,steps];ct=Array[pct,steps];
px=Array[plx,steps];pe=Array[ple,steps];
birth=Array[cbr,steps];
growth=Array[gag,steps];popgrowth=Array[nag,steps];cap=Array[k,steps+1];
ginix=Array[gix,steps];cfd=Array[fd,steps];]

Setting the initial conditions

initcond:=Module[{scal,cxtemp,glt,Nlt},
Do[Nh[i]=1/dim,{i,1,dim}];
gdist=LogNormalDistribution[0,σ];
cdfunction=CDF[gdist,y];
pdfunction=PDF[gdist,y];
inversecdf[z_]=y/.Solve[cdfunction==z,y][[1]];
Do[x[i]=inversecdf[(i-.5)/dim] ,{i,1,dim}];
scal=Sum[x[i] Nh[i],{i,1,dim}];
Do[x[i]=x[i]/scal ,{i,1,dim}];
Nlt=(((1-η)γ)/((ϕ-θ)(1+β+γ)));
glt=If[κ<1-τ,1+ϱ, B(((((η(ϕ-θ))/(1-η))))^η)];
k[1]=((((β/(1+β))(1/(glt Nlt))(1-α)))^(1/(1-α)))/capinifact;
hh[1]=If[κ<1-τ,(B/glt)(((((η(ϕ-θ))/(1-η))))^η),1];]

Dynamic simulation

simul:=Module[{a,q},
Do[
Print["Start simulation ",Date[]];Do[n[i]=If[x[i]>θ/(η ϕ),(((1-η) γ x[i])/((ϕ x[i]-θ)(1+β+γ))),((γ )/(ϕ (1+β+γ)))],{i,1,dim}];

nsort=Sort[Array[n,dim]];
fd[j]=-Mean[Take[nsort,dim/5]]+Mean[Take[nsort,-dim/5]];
nag[j]=Sum[n[i] Nh[i],{i,1,dim}];
cbr[j]=If[j>1,nag[j]nag[j-1]/(1+nag[j-1]),nag[j]nag[j]/(1+nag[j])];
Do[Nhf[i]=((n[i])/(nag[j])) Nh[i],{i,1,dim}];
Do[e[i]=Max[((η ϕ x[i]-θ)/(1-η)),0],{i,1,dim}];
gag[j]=Sum[ Nhf[i] B (((θ+e[i]))^η)((x[i])^τ)((hh[j])^(κ+τ-1)),{i,1,dim}];
hh[j+1]=hh[j] gag[j]/(1+ϱ);
Do[inc[i]=(x[i](1-ϕ n[i])-e[i]n[i])(1-α)((k[j])^α),{i,1,dim}];
sd=Sort[income];
gix[j]=N[((Length[sd])/((Length[sd]-1)))((2Apply[Plus,Table[(2q-Length[sd]-1)sd[[q]],{q,1,Length[sd]}]])/(2 ((Length[sd])^2) Mean[sd]))];
k[j+1]=(β/(1+β))(1/(gag[j] nag[j]))(1-α)((k[j])^α);
pct[j]=Sum[If[x[i]<θ/(η ϕ),Nh[i],0] ,{i,1,dim}];
Do[xf[i]=B (((θ+e[i]))^η)((x[i])^τ) ((hh[j])^(κ+τ-1))/ gag[j],{i,1,dim}];
Print["Start sorting households ",Date[]];
Clear[tp];
tp=Dot[Transpose[{Array[xf,dim]}],{{1,0}}]+Dot[Transpose[{Array[Nhf,dim]}],{{0,1}}];
tp2=Transpose[Sort[tp]];
Do[Nhf[i]=tp2[[2]][[i]];xf[i]=tp2[[1]][[i]],{i,1,dim}];
Print["Start computing the cumulative distribution ",Date[]];
Clear[CuNhf];
CuNhf[1]=Nhf[1];
Do[CuNhf[i]=CuNhf[i-1]+Nhf[i],{i,2,dim}];
Print["Start drawing random numbers ",Date[]];
Do[zNh[i]=1/dim;
Clear[z1,tp];tp=Random[];
z1=Ordering[Abs[Table[tp,{dim}]-Array[CuNhf,dim]],1][[1]];
zx[i]=xf[z1],{i,1,dim}];
Print["Start computing the error of approximation ",Date[]];
corr=Sum[zx[i]zNh[i],{i,1,dim}];
Print["step", j];
Do[x[i]=zx[i]/corr;Nh[i]=zNh[i],{i,1,dim}],{j,1,steps}]]

Results

result2:=Module[{p1,p2,p3,p4},
Clear[p5,p6,p7,p8,p9,p10];
p5=ListPlot[growth,AxesLabel->{t,g},DisplayFunction->Identity,AxesOrigin->{1,1}];
p6=ListPlot[popgrowth,AxesLabel->{t,n},DisplayFunction->Identity,AxesOrigin->{1,1}];
p7=ListPlot[cap,AxesLabel->{t,"k"},DisplayFunction->Identity];
p8=ListPlot[ginix,AxesLabel->{t,"gini"},DisplayFunction->Identity];
p9=ListPlot[ct,AxesLabel->{t,"% const"},DisplayFunction->Identity,AxesOrigin->{1,0}];
p10=ListPlot[cfd,AxesLabel->{t,"dif fert"},DisplayFunction->Identity];
Show[GraphicsGrid[{{p5,p6},{p7,p8},{p9,p10}}]]]

Main program

ϕ=.15/2;η=.634522;β=((.99)^120);α=1/3;θ=0.01187433277439054;γ=0.2712867133790373;
B=7.3667967472312402219159866854;ϱ=(1.02^30)-1;

Basis Simulations

Endogenous versus exogenous growth

dim=1000: trial    dim=10000 for best results

dim=10000;steps=8;capinifact=1;

σ=1;

τ=.2;κ=1-τ;defvar;initcond;simul;result2

τ=.05;κ=.1;defvar;initcond;simul;result2

τ=.2;κ=.1;defvar;initcond;simul;result2

τ=.2;κ=.65;defvar;initcond;simul;result2