Rich Nations, Poor Nations: How much can multiple equilibria explain? (CALIBRATION MATERIALS)
Bryan Graham, University of California - Berkeley
Jon Temple, Bristol University

This Mathematica notebook includes all the code required to replicate the calibration results and figures in Graham and Temple (2004). The code was executed in Mathematica 3.0. Please report any errors or suggestions for improvement to me by e-mail at:

bgraham@econ.berkeley.edu
or by mail at:

Department of Economics
University of California - Berkeley
549 Evans Hall #3880
Berkeley, CA 94720

I regret that I cannot provide any support to potential users. However, feel free to modify and use the code as you wish for any noncommerical purposes, just please cite our paper.

Thanks,

Bryan Graham

Loading of libraries and data

SetDirectory["c:\BSG_WORK\Research\RichPoorPaper"]\(<< Graphics`Graphics`
<< Graphics`MultipleListPlot`
<< 
        Statistics`ContinuousDistributions`\)<< Statistics`DescriptiveStatistics`<< Statistics`MultiDescriptiveStatistics`<< Statistics`LinearRegression`\(<< \ Statistics`NonlinearFit`
<< NumericalMath`ListIntegrate`\)

Off[General::spell1]

$TextStyle  = 
       {FontFamily  ->  "Courier",  FontSize  ->  12}

Nonparametric statistical modules

This section includes code for all the nonparametric statistical techniques used in Graham and Temple (2003). The coding is heavily commented and variable names are usually named in a way consistent with the corresponding theoretical literature. Procedures below include a module for univariate density estimatiom via the Rosenblatt-Parzen Kernel Estimator (Pagan and Ullah 1999:9), a module for graphing the nonparametric density estimate, and modules for optimal bandwidth selection via the Sheather and Jones (1991) "plug in" method and the simple "rule of thumb" method suggested by Silverman (1986). Some basic routines for the inequality decompositions reported in the paper are also included.

Nonparametric kernel density estimation using the Rosenblatt-Parzen kernel estimator as decribed by Pagan and Ullah (1999, Chapter 2).

KDensity[ed_, n_, win_, rescale_, EvalMin_, EvalMax_, NumberEval_, 
          SigL_]:=Module[{v, h, fHat, KernelResult, rd,  Minimum,  Maximum, 
              IntegralOfKernelSquared, fLowConfidence, fHighConfidence,  
              Zcrit},
  
  
    		  
  
  		  		  
gauss[v_] := (1/(((2  π))^(1/2))) Exp[-(v^2)/2];
   
Sort[ed]; 
Minimum := ed[[1]];
Maximum := ed[[n]];
  
rd = If[rescale == 1, 
                  Table[(((ed[[i]] - Minimum)/(Maximum - Minimum)))
                      , {i, 1, n, 1}], ed]; 
  
h = If[win == 0, 
                  0.9*Min[InterquartileRange[rd]/1.34, StandardDeviation[rd]]*
                    (n^(-1/5)), win]; 
  
fHat = 
                Table[(1/(n*h))*
                    (∑_(j = 1)^n)(
                        gauss[((EvalMin + 
                              i*(((EvalMax - EvalMin)/NumberEval)) - 
                              rd[[j]])/h)]), {i, 1, NumberEval, 1}]
                  ; 
  
                  		Zcrit = Quantile[NormalDistribution[0, 1], 1 - SigL/2];
IntegralOfKernelSquared = 
                NIntegrate[
                  (gauss[x])^2, {x, -∞, ∞}];
fLowConfidence = 
                Table[fHat[[i]] - 
                    Zcrit*(((n*h))^(-1/2))*
                      (fHat[[i]]*IntegralOfKernelSquared)^(1/2), {
                    i, 1, NumberEval, 1}];
fHighConfidence = 
                Table[fHat[[i]] + 
                    Zcrit*(((n*h))^(-1/2))*
                      (fHat[[i]]*IntegralOfKernelSquared)^(1/2), {
                    i, 1, NumberEval, 1}];	  
KernelResult = 
                Table[0, {i, 1, NumberEval, 1}, {j, 1, 4, 1}];
Do[KernelResult[[i, 1]] = 
                  EvalMin + i*(((EvalMax - EvalMin)/NumberEval));  
                KernelResult[[i, 2]] = fHat[[i]]; 
                
KernelResult[[i, 3]] = 
                  fLowConfidence[[i]]; 
                KernelResult[[i, 4]] = fHighConfidence[[i]]; 
                
TabulatedResult = CalibrationResult, {i, 
                  1, NumberEval, 1}];
Sort[KernelResult]
];

Bandwidth selection for univariate density estimation using the Sheather and Jones (1991) plug-in method as described in their article.

KDenBandPISheather[ed_, n_, win_, rescale_, numiter_, precision_]:=Module[{v, h, rd,  Minimum,  Maximum,  RK, σ2, IQR, a, b, Sa, 
              Tb, x, hstar},
  
  
  		  		  
                
  		  
  
                
gauss[v_] := (1/(((2*3.14))^(1/2))) Exp[-(v^2)/2];
  		RK = NIntegrate[
                  (gauss[x])^2, {x, -∞, ∞}];
σ2 = 
                NIntegrate[
                  (x^2)*(gauss[x]), {x, -∞, ∞}];
Sort[ed];
Minimum := ed[[1]];
Maximum := ed[[n]];
rd = If[rescale == 1, 
                  Table[(((ed[[i]] - Minimum)/(Maximum - Minimum)))
                      , {i, 1, n, 1}], ed];
IQR = InterquartileRange[rd];
  		a = 0.920*IQR*(n^(-1/7));
b = 0.912*IQR*(n^(-1/9));
Tb = (((n*(n*1)))^(-1))*(b^(-7)) 
                  (∑_(i = 1)^n)
                      (∑_(j = 1)^n)(
                        gauss[((rd[[i]] - rd[[j]])/b)]);
Sa = (((n*(n*1)))^(-1))*(a^(-5)) 
                  (∑_(i = 1)^n)
                      (∑_(j = 1)^n)(
                        gauss[((rd[[i]] - rd[[j]])/a)]);
α2[h_] := 1.357*(((Sa/Tb))^(1/7))*(h^(5/7));
Sα[h_] := 
                (((n*(n*1)))^(-1))*α2[h]^(-5) 
                  (∑_(i = 1)^n)
                      (∑_(j = 1)^n)(
                        gauss[((rd[[i]] - rd[[j]])/α)2[h]])
                        ;
h = If[win == 0, 
                  0.9*Min[InterquartileRange[rd]/1.34, StandardDeviation[rd]]*
                    (n^(-1/5)), win];
hstar = 
                x /. First[
                    FindRoot[
                      (((RK/(σ(2^2)*Sα[x])))^(1/5))*
                          (n^(-1/5)) - x, {x, h}, 
                      WorkingPrecision -> precision, 
                      MaxIterations -> numiter]]
];

Bandwidth selection for univariate density estimation using Silverman's (1986) "rule-of-thumb" as decribed in Pagan and Ullah (1999, Chapter 2).

KDenBandPISilverman[ed_, n_, rescale_]:=Module[{v, h, rd,  Minimum,  Maximum},		  	Sort[ed]; 
            
Minimum := ed[[1]]; 
            
Maximum := ed[[n]]; 
            
rd = 
              If[rescale == 1, 
                Table[(((ed[[i]] - Minimum)/(Maximum - Minimum)))
                    , {i, 1, n, 1}], ed]; 
            
h = 
              0.9*Min[InterquartileRange[rd]/1.34, StandardDeviation[rd]]*
                (n^(-1/5))
];

Graphing routine for nonparametric estimators

GraphKernelEstimate[KernelEstimate_, data_, n_, Confidence_, 
          IncludeData_, xAxisLabel_, yAxisLabel_, title_, darkness_, xMin_, 
          xMax_,  yMin_, yMax_]:=Module[{fHat, fLowConfidence, fHighConfidence, FinalPlot},
  		  		  
  
                
  
  
  
  		  
fHat = 
              Table[If[j == 1, KernelEstimate[[i, 1]], 
                  KernelEstimate[[i, 2]]], {i, 1, n, 1}, {j, 1, 2, 1}]; 
            
fLowConfidence = 
              Table[If[j == 1, KernelEstimate[[i, 1]], 
                  KernelEstimate[[i, 3]]], {i, 1, n, 1}, {j, 1, 2, 1}]; 
            fHighConfidence = 
              Table[If[j == 1, KernelEstimate[[i, 1]], 
                  KernelEstimate[[i, 4]]], {i, 1, n, 1}, {j, 1, 2, 1}]; 
            
FinalPlot = 
              If[Confidence == 1, 
                If[IncludeData == 1, 
                  MultipleListPlot[fHat, fLowConfidence, fHighConfidence, 
                    data, PlotJoined -> {True, True, True, False}, 
                    Frame -> True, 
                    FrameLabel 
                      -> {xAxisLabel, yAxisLabel, title, ""}, 
                    Axes -> False, 
                    SymbolShape 
                      -> {None, None, None, PlotSymbol[Diamond]}, 
                    PlotStyle 
                      -> {{GrayLevel[darkness], Thickness[0.0075]}, {
                          Dashing[{0.025, 0.0125}], GrayLevel[darkness], 
                          Thickness[0.005]}, {Dashing[{0.025, 0.0125}], 
                          GrayLevel[darkness], Thickness[0.005]}, {
                          GrayLevel[0], Thickness[0.005]}}, 
                    PlotRange -> {{xMin, xMax}, {yMin, yMax}}], 
                  MultipleListPlot[fHat, fLowConfidence, fHighConfidence, 
                    PlotJoined -> True, Frame -> True, 
                    FrameLabel 
                      -> {xAxisLabel, yAxisLabel, title, ""}, 
                    Axes -> False, 
                    SymbolShape -> {None, None, None}, 
                    PlotStyle 
                      -> {{GrayLevel[darkness], Thickness[0.0075]}, {
                          Dashing[{0.025, 0.0125}], GrayLevel[darkness], 
                          Thickness[0.005]}, {Dashing[{0.025, 0.0125}], 
                          GrayLevel[darkness], Thickness[0.005]}}, 
                    PlotRange -> {{xMin, xMax}, {yMin, yMax}}]], 
                If[IncludeData == 1, 
                  MultipleListPlot[fHat, data, 
                    PlotJoined -> {True, False}, Frame -> True, 
                    FrameLabel 
                      -> {xAxisLabel, yAxisLabel, title, ""}, 
                    Axes -> False, 
                    SymbolShape -> {None, PlotSymbol[Diamond]}, 
                    PlotStyle 
                      -> {{GrayLevel[darkness]}, {GrayLevel[0]}}, 
                    PlotRange -> {{xMin, xMax}, {yMin, yMax}}], 
                  ListPlot[fHat, PlotJoined -> True, Frame -> True, 
                    FrameLabel 
                      -> {xAxisLabel, yAxisLabel, title, ""}, 
                    Axes -> False, 
                    PlotStyle 
                      -> {GrayLevel[darkness], Thickness[0.0075]}, 
                    PlotRange -> {{xMin, xMax}, {yMin, yMax}}]]]
];

A simple routine to generate the empirical cummulative density function of a vector of data

(* This function takes a list of data, Density_, and its length, n_, and returns an n_ by two matrix of the cummulative distribution function of the data *)

CummDenF[Density_, n_] := 
      Table[If[j == 1, Density[[i]], i/n], {i, 1, n, 1}, {j, 1, 2, 1}];

Routines for the calculation of Atkinson's (1970) inequality measures, the Gini coefficient, and the mean logarithmic deviation. Note that one half times the coefficient of variation squared is equal to the generalized entropy measure of inequality when c = 2.

\(AtkinsonMeasure[y_, n_, ϵ_] := 
      1 - \((Sum[
              \((y[\([i]\)]\/Mean[y])\)\^\(1 - ϵ\)*1\/n, {i, 1, n, 
                1}])\)\^\(1/\((1 - ϵ)\)\); 

    Gini[y_, n_] := 
      \(1\/\(2*n\^2*Mean[y]\)\) 
        \(∑\+\(i = 1\)\%n
            \(∑\+\(j = 1\)\%n Abs[y[\([i]\)] - y[\([j]\)]]\)\); 

    CoefOfVariation[y_, n_] := 
      \((\(1\/\(n*Mean[y]\^2\)\) 
            \(∑\+\(i = 1\)\%n\((y[\([i]\)] - Mean[y])\)\^2\))\)\^\(1/2
          \); 
MeanLogDev[y_, n_] := 
      \((\(1\/n\) \(∑\+\(i = 1\)\%n\((Log[Mean[y]/y[\([i]\)]])\)\))\); 

    RelativeMeanDeviation[y_, n_] := 
      \((\(1\/n\) \(∑\+\(i = 1\)\%n Abs[y[\([i]\)]/Mean[y] - 1]\))\); 

    Entropy[y_, n_, c_] := 
      \((1\/n)\)*\((1\/\(c*\((c - 1)\)\))\)*
        \((Sum[\((y[\([i]\)]\/Mean[y])\)\^c - 1, {i, 1, n, 1}])\); 

    Theil[y_, n_] := 
      \((1\/n)\)*
        Sum[\((y[\([i]\)]\/Mean[y])\)*Log[y[\([i]\)]\/Mean[y]], {i, 1, n, 
            1}]; \)

The following functions decompose the sources of inequality into between and within components as in Shorrocks (1980). The procedure works for the case of two groups, see Mookerjee and Shorrocks (1982) for an example on UK income data.

DCMeanLogDev[y_, n_, d_]:=Module[{v0, L0, v1, L1, I0, w0, w1, b, decomp},		  
  
  		  		v0 = N[Sum[If[d[[i]] == 0, 1, 0], {i, 1, n, 1}]/n, 10];
v1 = N[1 - v0, 10];
L0 = (
                    Sum[If[d[[i]] == 0, y[[i]], 0], {i, 1, n, 1}]/
                      (v0*n))/Mean[y];		L1 = (
                    Sum[If[d[[i]] == 1, y[[i]], 0], {i, 1, n, 1}]/
                      (v1*n))/Mean[y];
I0 = MeanLogDev[y, n];
w0  =  
                v0*((1/(v0*n))*
                      Sum[If[d[[i]] == 0, 
                          (Log[(L0*Mean[y])/y[[i]]]), 0], {i, 1, 
                          n, 1}]);		w1  = 
                v1*((1/(v1*n))*
                      Sum[If[d[[i]] == 1, 
                          (Log[(L1*Mean[y])/y[[i]]]), 0], {i, 1, 
                          n, 1}]);
b = v0*Log[1/L0] + v1*Log[1/L1];
  	decomp = {I0, w0/I0, w1/I0, b/I0, v0, L0*Mean[y], v1, 
                  L1*Mean[y]}
];

DCEntropy[y_, n_, d_, c_]:=Module[{v0, L0, v1, L1, I0, w0, w1, b, decomp},		  
  
  		  			v0 = N[Sum[If[d[[i]] == 0, 1, 0], {i, 1, n, 1}]/n, 10];
v1 = N[1 - v0, 10];
L0 = (
                    Sum[If[d[[i]] == 0, y[[i]], 0], {i, 1, n, 1}]/
                      (v0*n))/Mean[y];	L1 = (
                    Sum[If[d[[i]] == 1, y[[i]], 0], {i, 1, n, 1}]/
                      (v1*n))/Mean[y];
I0 = Entropy[y, n, c];
w0  =  
                v0*(L0^c)*
                  (((1/(v0*n)))*((1/(c*(c - 1))))*
                      (Sum[If[d[[i]] == 0, 
                            (((y[[i]]/(L0*Mean[y])))^c) - 1, 0], {i, 
                            1, n, 1}]));		w1  = 
                v1*(L1^c)*
                  (((1/(v1*n)))*((1/(c*(c - 1))))*
                      (Sum[If[d[[i]] == 1, 
                            (((y[[i]]/(L1*Mean[y])))^c) - 1, 0], {i, 
                            1, n, 1}]));		b = (1/(c*(c - 1)))*
                  (v0*((L0^c) - 1) + v1*((L1^c) - 1));		  
decomp = {I0, w0/I0, w1/I0, b/I0, v0, L0*Mean[y], v1, 
                  L1*Mean[y]}
	];

DCTheil[y_, n_, d_]:=Module[{v0, L0, v1, L1, I0, w0, w1, b, decomp},		  
  
  
  		v0 = N[Sum[If[d[[i]] == 0, 1, 0], {i, 1, n, 1}]/n, 10];
v1 = N[1 - v0, 10];
L0 = (
                    Sum[If[d[[i]] == 0, y[[i]], 0], {i, 1, n, 1}]/
                      (v0*n))/Mean[y];		L1 = (
                    Sum[If[d[[i]] == 1, y[[i]], 0], {i, 1, n, 1}]/
                      (v1*n))/Mean[y];
I0 = Theil[y, n];
w0  =  
                v0*L0*(
                    ((1/(v0*n)))*
                      (Sum[If[d[[i]] == 0, 
                            (y[[i]]/(L0*Mean[y]))*
                              Log[(y[[i]]/(L0*Mean[y]))], 0], {i, 
                            1, n, 1}]));		w1  = 
                v1*L1*(
                    ((1/(v1*n)))*
                      (Sum[If[d[[i]] == 1, 
                            (y[[i]]/(L1*Mean[y]))*
                              Log[(y[[i]]/(L1*Mean[y]))], 0], {i, 
                            1, n, 1}]));		b = v0*L0*Log[L0] + v1*L1*Log[L1];		  
decomp = {I0, w0/I0, w1/I0, b/I0, v0, L0*Mean[y], v1, 
                  L1*Mean[y]}
];

Calibration subroutines

The following modules implement the calibration strategy outlined in the appendix of the paper. FullCalibration provides a complete set of calibration results, while WelfareCalibration and TFPCalibration provide just welfare and TFP results respectively. Also included are a shortcut alternative equilibrium assignment subroutine and some code for the aggregated growth accounting exercise done in Section 7 of the paper.

FullCalibration[rdlist_, nc_, γ_, α_, β_, 
          λ_, cθ_]:=Module[{dlist, cslist, θlist, blista, blistb, blist, zlist, b, 
              CalibrationResult, ylist, traplist, TFPm, TFPa, TabulatedResult}
              ,
  
		  		  
  
  
  
                  
		  
dlist = rdlist/100; 
cslist = 
                Table[(((1 - dlist[[i, 2]])/dlist)[[i, 2]]) 
                      (dlist[[i, 1]]/(1 - dlist[[i,  1]])) 
                      (((1 - γ)/(1 - α - β))) - 
                    1, {i, 1, nc, 1}];
θlist = 
                Table[(
                    (((dlist[[i, 1]]/(1 - dlist[[i, 1]])) 
                        (((1 - γ)/(1 - α - β))) 
                        ((α/γ)) 
                        (((1 + cθ)/(1 + cslist[[i]])))
                          )/(1 + 
                        (dlist[[i, 1]]/(1 - dlist[[i, 1]])) 
                          (((1 - γ)/(1 - α - β))) 
                          ((α/γ)) 
                          (((1 + cθ)/(1 + cslist[[i]]))))))
                    , {i, 1, nc, 1}];
blista = 
                Table[b  /. 
                    First[FindRoot[
                        1 == (((dlist[[i,  1]]/b))^β) 
                            (((((1 - dlist[[i, 1]])/(1 - b)))
                                  )^λ) 
                            ((((((1 - 
                                      dlist[[i,  1]] 
                                        (1 - 
                                        (((1 - γ)/(1 - α - 
                                        β))) ((α/γ)) 
                                        (((1 + cθ)/(1 + 
                                        cslist[[i]]))))))/((
                                    1 - b 
                                        (1 - 
                                        (((1 - γ)/(1 - α - 
                                        β))) ((α/γ)) 
                                        (((1 + cθ)/(1 + 
                                        cslist[[i]]))))))))
                                )^(α - 
                                γ (1 + λ))), {b, 
                          0.9999999999},  AccuracyGoal -> 10,  
                        WorkingPrecision -> 40, MaxIterations -> 100]], {i, 
                    1, nc, 1}];blistb = 
                Table[b  /. 
                    First[FindRoot[
                        1 == (((dlist[[i,  1]]/b))^β) 
                            (((((1 - dlist[[i, 1]])/(1 - b)))
                                  )^λ) 
                            ((((((1 - 
                                      dlist[[i,  1]] 
                                        (1 - 
                                        (((1 - γ)/(1 - α - 
                                        β))) ((α/γ)) 
                                        (((1 + cθ)/(1 + 
                                        cslist[[i]]))))))/((
                                    1 - b 
                                        (1 - 
                                        (((1 - γ)/(1 - α - 
                                        β))) ((α/γ)) 
                                        (((1 + cθ)/(1 + 
                                        cslist[[i]]))))))))
                                )^(α - 
                                γ (1 + λ))), {b, 
                          0.0000000001},  AccuracyGoal -> 10,  
                        WorkingPrecision -> 40, MaxIterations -> 100]], {i, 
                    1, nc, 1}];blist = 
                Table[If[
                    Abs[blista[[i]] - dlist[[i, 1]]] < 
                      Abs[blistb[[i]] - dlist[[i, 1]]], 
                    blistb[[i]], blista[[i]]], {i, 1, nc, 1}];		  
traplist = 
                Table[If[blist[[i]] < dlist[[i, 1]], 1, 0], {i, 1, 
                    nc, 1}];
  
zlist = 
                Table[(
                    (((blist[[i]]/(1 - blist[[i]])) 
                        (((1 - γ)/(1 - α - β))) 
                        ((α/γ)) 
                        (((1 + cθ)/(1 + cslist[[i]])))
                          )/(1 + 
                        (blist[[i]]/(1 - blist[[i]])) 
                          (((1 - γ)/(1 - α - β))) 
                          ((α/γ)) 
                          (((1 + cθ)/(1 + cslist[[i]]))))))
                    , {i, 1, nc, 1}];
  
ylist = 
                Table[Log[
                      (((1 + 
                              (blist[[i]]/(1 - blist[[i]])) 
                                (((1 - γ)/(1 - α - 
                                      β))) ((α/γ)) 
                                (((1 + cθ)/(1 + cslist[[i]])))
                                  )/(1 + 
                              (dlist[[i, 1]]/(1 - dlist[[i, 1]]))
                                   (
                                  ((1 - γ)/(1 - α - 
                                      β))) ((α/γ)) 
                                (((1 + cθ)/(1 + cslist[[i]])))
                                  ))) 
                        (((((1 - blist[[i]])/(1 - 
                                  dlist[[i, 1]])))
                              )^((1 - γ) 
                              (1 + λ))) 
                        (((((1 - zlist[[i]])/(1 - 
                                θlist[[i]]))))^(γ 
                            (1 + λ)))] + rdlist[[i, 3]], {i, 
                    1, nc, 1}];
  
TFPm = 
                Table[((1/((1 - γ)*(1 + λ))))*
                      (Log[1 - rdlist[[i, 2]]/100] + 
                          rdlist[[i, 3]] + rdlist[[i, 6]]) - 
                    (γ/(1 - γ))*
                      (Log[1 - θlist[[i]]] + 
                          rdlist[[i, 4]] + rdlist[[i, 6]]) - 
                    rdlist[[i, 5]] - 
                    (Log[1 - rdlist[[i, 1]]/100] + rdlist[[i, 6]])
                      , {i, 1, nc, 1}];	
  
                  	TFPa = 
                Table[(1/(1 - α - β)) 
                      (Log[rdlist[[i, 2]]/100] + rdlist[[i, 3]] + 
                          rdlist[[i, 6]]) - 
                    ((α )/(1 - α - β)) 
                      (Log[θlist[[i]]] + rdlist[[i, 4]] + 
                          rdlist[[i, 6]]) - 
                    (β/(1 - α - β)) 
                      (rdlist[[i, 7]]) -  rdlist[[i, 5]] - 
                    (Log[rdlist[[i, 1]]/100] + rdlist[[i, 6]]), {
                    i, 1, nc, 1}];
  
CalibrationResult = Table[0, {i, 1, nc, 1}, {j, 1, 15, 1}];
Do[CalibrationResult[[t, 1]] = traplist[[t]];  
                CalibrationResult[[t, 2]] = 
                  If[traplist[[t]] == 1, rdlist[[t, 3]], 
                    ylist[[t]]]; 
                
CalibrationResult[[t, 3]] = 
                  If[traplist[[t]] == 0, rdlist[[t, 3]], 
                    ylist[[t]]]; 
                CalibrationResult[[t, 4]] = dlist[[t, 1]]; 
                CalibrationResult[[t, 5]] = blista[[t]]; 
                CalibrationResult[[t, 6]] = blistb[[t]]; 
                CalibrationResult[[t, 7]] = θlist[[t]]; 
                CalibrationResult[[t, 8]] = 
                  If[traplist[[t]] == 1, θlist[[t]], 
                    zlist[[t]]]; 
                CalibrationResult[[t, 9]] = 
                  If[traplist[[t]] == 0, θlist[[t]], 
                    zlist[[t]]]; 
                CalibrationResult[[t, 10]] = cslist[[t]];  
                CalibrationResult[[t, 11]] = TFPm[[t]];  
                CalibrationResult[[t, 12]] = TFPa[[t]];  
                CalibrationResult[[t, 13]] = λ; 
                CalibrationResult[[t, 14]] = rdlist[[t, 22]]; 
                CalibrationResult[[t, 15]] = rdlist[[t, 8]]; 
                , {t, 1, nc, 1}];
TabulatedResult = CalibrationResult];

\(\(WelfareCalibration[rdlist_, nc_, γ_, α_, β_, 
        λ_, cθ_, wbcode_] := 
      Module[{dlist, cslist, θlist, blista, blistb, blist, zlist, 
          wlist, b}, 
dlist = rdlist/100; \ 
        
cslist = 
          Table[\((\(1 - dlist[\([i, 2]\)]\)\/dlist[\([i, 2]\)])\) 
                \((dlist[\([i, 1]\)]\/\(1 - dlist[\([i, \ 1]\)]\))\) 
                \((\(1 - γ\)\/\(1 - α - β\))\) - 1, {i, 
              1, nc, 1}]; 
        
θlist = 
          Table[\((
              \(\((dlist[\([i, 1]\)]\/\(1 - dlist[\([i, 1]\)]\))\) 
                  \((\(1 - γ\)\/\(1 - α - β\))\) 
                  \((α\/γ)\) 
                  \((\(1 + cθ\)\/\(1 + cslist[\([i]\)]\))\)\)\/\(1 + 
                  \((dlist[\([i, 1]\)]\/\(1 - dlist[\([i, 1]\)]\))\) 
                    \((\(1 - γ\)\/\(1 - α - β\))\) 
                    \((α\/γ)\) 
                    \((\(1 + cθ\)\/\(1 + cslist[\([i]\)]\))\)\))\), {
              i, 1, nc, 1}]; 
        
blista = 
          Table[b\  /. 
              First[FindRoot[
                  1 == \(\((dlist[\([i, \ 1]\)]\/b)\)\^β\) 
                      \(\((\(1 - dlist[\([i, 1]\)]\)\/\(1 - b\))
                            \)\^λ\) 
                      \((\((1 - 
                                dlist[\([i, \ 1]\)] 
                                  \((1 - 
                                      \((\(1 - γ\)\/\(1 - α - 
                                        β\))\) \((α\/γ)\) 
                                        \((\(1 + cθ\)\/\(1 + 
                                        cslist[\([i]\)]\))\))\))\)\/\((
                              1 - b \((
                                    1 - \((
                                        \(1 - γ\)\/\(1 - α - 
                                        β\))\) \((α\/γ)\) 
                                        \((\(1 + cθ\)\/\(1 + 
                                        cslist[\([i]\)]\))\))\))\))
                          \)\^\(α - γ \((1 + λ)\)\), {
                    b, 0.9999999999}, \ AccuracyGoal -> 10, \ 
                  WorkingPrecision -> 40, MaxIterations -> 100]], {i, 1, nc, 
              1}]; blistb = 
          Table[b\  /. 
              First[FindRoot[
                  1 == \(\((dlist[\([i, \ 1]\)]\/b)\)\^β\) 
                      \(\((\(1 - dlist[\([i, 1]\)]\)\/\(1 - b\))
                            \)\^λ\) 
                      \((\((1 - 
                                dlist[\([i, \ 1]\)] 
                                  \((1 - 
                                      \((\(1 - γ\)\/\(1 - α - 
                                        β\))\) \((α\/γ)\) 
                                        \((\(1 + cθ\)\/\(1 + 
                                        cslist[\([i]\)]\))\))\))\)\/\((
                              1 - b \((
                                    1 - \((
                                        \(1 - γ\)\/\(1 - α - 
                                        β\))\) \((α\/γ)\) 
                                        \((\(1 + cθ\)\/\(1 + 
                                        cslist[\([i]\)]\))\))\))\))
                          \)\^\(α - γ \((1 + λ)\)\), {
                    b, 0.0000000001}, \ AccuracyGoal -> 10, \ 
                  WorkingPrecision -> 40, MaxIterations -> 100]], {i, 1, nc, 
              1}]; blist = 
          Table[If[
              Abs[blista[\([i]\)] - dlist[\([i, 1]\)]] < 
                Abs[blistb[\([i]\)] - dlist[\([i, 1]\)]], blistb[\([i]\)], 
              blista[\([i]\)]], {i, 1, nc, 1}]; 
        
zlist = 
          Table[\((
              \(\((blist[\([i]\)]\/\(1 - blist[\([i]\)]\))\) 
                  \((\(1 - γ\)\/\(1 - α - β\))\) 
                  \((α\/γ)\) 
                  \((\(1 + cθ\)\/\(1 + cslist[\([i]\)]\))\)\)\/\(1 + 
                  \((blist[\([i]\)]\/\(1 - blist[\([i]\)]\))\) 
                    \((\(1 - γ\)\/\(1 - α - β\))\) 
                    \((α\/γ)\) 
                    \((\(1 + cθ\)\/\(1 + cslist[\([i]\)]\))\)\))\), {
              i, 1, nc, 1}]; 		
        wlist = Table[
            If[j == 1, rdlist[\([i, 3]\)], 
              If[j == 2, 
                Log[\((\(1 + 
                            \((blist[\([i]\)]\/\(1 - blist[\([i]\)]\))\) 
                              \((\(1 - γ\)\/\(1 - α - β\))
                                \) \((α\/γ)\) 
                              \((\(1 + cθ\)\/\(1 + cslist[\([i]\)]\))
                                \)\)\/\(1 + 
                            \((dlist[\([i, 1]\)]\/\(1 - dlist[\([i, 1]\)]\))
                                \) \((
                                \(1 - γ\)\/\(1 - α - β\))
                                \) \((α\/γ)\) 
                              \((\(1 + cθ\)\/\(1 + cslist[\([i]\)]\))
                                \)\))\) 
                      \(\((\(1 - blist[\([i]\)]\)\/\(1 - dlist[\([i, 1]\)]\))
                            \)\^\(\((1 - γ)\) \((1 + λ)\)\)\) 
                      \((\(1 - zlist[\([i]\)]\)\/\(1 - 
                              θlist[\([i]\)]\))\)\^\(γ 
                          \((1 + λ)\)\)] + rdlist[\([i, 3]\)], 
                rdlist[\([i, wbcode]\)]]], {i, 1, nc, 1}, {j, 1, 3, 1}]
	]; 
    \)\)

TFPCalibration[rdlist_, nc_, γ_, α_, β_, 
          λ_, cθ_]:=Module[{dlist, cslist, θlist, TFPm, TFPa, CalibrationResult, 
              TabulatedResult},
  
\(\(dlist = rdlist/100; 
                \)
\(cslist = 
                  Table[\((\(1 - dlist[\([i, 2]\)]\)\/dlist[\([i, 2]\)])\)*
                        \((dlist[\([i, 1]\)]\/\(1 - dlist[\([i, \ 1]\)]\))\)*
                        \((\(1 - γ\)\/\(1 - α - β\))\) - 
                      1, {i, 1, nc, 1}]; 
                \)
\(θlist = 
                  Table[\((
                      \(\((dlist[\([i, 1]\)]\/\(1 - dlist[\([i, 1]\)]\))\)*
                          \((\(1 - γ\)\/\(1 - α - β\))\)*
                          \((α\/γ)\)*
                          \((\(1 + cθ\)\/\(1 + cslist[\([i]\)]\))
                            \)\)\/\(1 + 
                          \((dlist[\([i, 1]\)]\/\(1 - dlist[\([i, 1]\)]\))\)*
                            \((\(1 - γ\)\/\(1 - α - β\))
                              \)*\((α\/γ)\)*
                            \((\(1 + cθ\)\/\(1 + cslist[\([i]\)]\))
                              \)\))\), {i, 1, nc, 1}]; 
                \)
TFPm = 
              Table[\((1\/\(\((1 - γ)\)*\((1 + λ)\)\))\)*
                    \((Log[1 - rdlist[\([i, 2]\)]/100] + rdlist[\([i, 3]\)] + 
                        rdlist[\([i, 6]\)])\) - 
                  γ\/\(1 - γ\)*
                    \((Log[1 - θlist[\([i]\)]] + rdlist[\([i, 4]\)] + 
                        rdlist[\([i, 6]\)])\) - rdlist[\([i, 5]\)] - 
                  \((Log[1 - rdlist[\([i, 1]\)]/100] + rdlist[\([i, 6]\)])
                    \), {i, 1, nc, 1}]; 		
		
            TFPa = Table[
                \(1\/\(1 - α - β\)\) 
                    \((Log[rdlist[\([i, 2]\)]/100] + rdlist[\([i, 3]\)] + 
                        rdlist[\([i, 6]\)])\) - 
                  \(\(α\ \)\/\(1 - α - β\)\) 
                    \((Log[θlist[\([i]\)]] + rdlist[\([i, 4]\)] + 
                        rdlist[\([i, 6]\)])\) - 
                  β\/\(1 - α - β\)\ 
                    \((rdlist[\([i, 7]\)])\) - \ rdlist[\([i, 5]\)] - 
                  \((Log[rdlist[\([i, 1]\)]/100] + rdlist[\([i, 6]\)])\), {i, 
                  1, nc, 1}]; 	
            
CalibrationResult = 
              Table[0, {i, 1, nc, 1}, {j, 1, 8, 1}]; 
            
Do[
              CalibrationResult[\([t, 1]\)] = TFPa[\([t]\)]; \ 
              CalibrationResult[\([t, 2]\)] = TFPm[\([t]\)]; 
              
CalibrationResult[\([t, 3]\)] = 
                dlist[\([t, 1]\)]; 
              CalibrationResult[\([t, 4]\)] = θlist[\([t]\)]; 
              CalibrationResult[\([t, 5]\)] = cslist[\([t]\)]; 
              CalibrationResult[\([t, 6]\)] = rdlist[\([t, 3]\)]; 
              CalibrationResult[\([t, 7]\)] = rdlist[\([t, 8]\)]; 
              \(CalibrationResult[\([t, 8]\)] = rdlist[\([t, 22]\)]; \), 
              \ {t, 1, nc, 1}]; 
            
TabulatedResult = CalibrationResult\)];

AggegrateTFP[rdlist_, nc_, γ_, α_, β_, 
          λ_, cθ_, YL_, KL_, HL_, DL_, TFP_, HJ_]:=Module[{CR, v, aA, bA, aH, bH, FD, dA, dAh, A, Ah, Aj, t},		  		CR = FullCalibration[rdlist, nc, γ, α, β, 
                  λ, cθ];
  		aA = Table[
                  (rdlist[[i, 2]]/100)*α + 
                    (1 - rdlist[[i, 2]]/100)*γ, {i, 1, nc, 1}]
                  ;bA = Table[
                  (rdlist[[i, 2]]/100)*β, {i, 1, nc, 1}];		Print["***Basic Statistics for Actual World Distribution of \
Output***"];		Print["Mean capital share : ", N[Mean[aA], 4], 
                " Standard deviation: ", N[StandardDeviation[aA], 4]];		Print["Mean land share    : ", N[Mean[bA], 4], 
                " Standard deviation: ", N[StandardDeviation[bA], 4]];			Print["Mean log Y/L       : ", N[Mean[YL], 4], 
                " Standard deviation: ", N[StandardDeviation[YL], 4]];		  		v = Table[
                  v  /. First[
                      Solve[(v/(1 - v)) == 
                          (CR[[i, 6]]/(1 - CR[[i, 6]]))*
                            (((rdlist[[i, 2]]/100)/(1 - 
                                  rdlist[[i, 2]]/100)))*
                            (((1 - 
                                  rdlist[[i, 1]]/100)/(rdlist[
                                    [i, 1]]/100))), v]], {i, 1, nc, 1}]
                  ;	aH = Table[
                  v[[i]]*α + (1 - v[[i]])*γ, {i, 1, 
                    nc, 1}];
bH = Table[v[[i]]*β, {i, 1, nc, 1}];			Print["***Basic Statistics for Counterfactual World \
Distribution of Output***"];			Print["Mean capital share : ", N[Mean[aH], 4], 
                " Standard deviation: ", N[StandardDeviation[aH], 4]];		Print["Mean land share    : ", N[Mean[bH], 4], 
                " Standard deviation: ", N[StandardDeviation[bH], 4]];		  		FD = Table[
                  If[j == 1, KL[[i]], 
                    If[j == 2, DL[[i]], 
                      If[j == 3, HL[[i]], 
                        If[j == 4, YL[[i]], 
                          If[j == 5, CR[[i, 3]], TFP[[i]]]]]]], {i, 
                    1, nc, 1}, {j, 1, 6, 1}];
Print["Mean log Y/L       : ", 
                N[Mean[Table[FD[[i, 5]], {i, 1, nc, 1}]], 4], 
                " Standard deviation: ", 
                N[StandardDeviation[Table[FD[[i, 5]], {i, 1, nc, 1}]], 4]]
                ;
		FD = Sort[FD];
  		dA = Table[
                  If[i == 1, 0, 
                    (FD[[i, 4]] - FD[[i - 1, 4]]) - 
                      0.5*(aA[[i]] + aA[[i - 1]])*
                        (FD[[i, 1]] - FD[[i - 1, 1]]) - 
                      0.5*(bA[[i]] + bA[[i - 1]])*
                        (FD[[i, 2]] - FD[[i - 1, 2]]) - 
                      0.5*
                        ((1 - aA[[i]] - bA[[i]]) + 
                            (1 - aA[[i - 1]] - bA[[i - 1]]))*
                        (FD[[i, 3]] - FD[[i - 1, 3]])], {i, 1, 
                    nc, 1}];			dAh = Table[
                  If[i == 1, 0, 
                    (FD[[i, 5]] - FD[[i - 1, 5]]) - 
                      0.5*(aH[[i]] + aH[[i - 1]])*
                        (FD[[i, 1]] - FD[[i - 1, 1]]) - 
                      0.5*(bH[[i]] + bH[[i - 1]])*
                        (FD[[i, 2]] - FD[[i - 1, 2]]) - 
                      0.5*((1 - aH[[i]] - bH[[i]]) + 
                            (1 - aH[[i - 1]] - bH[[i - 1]]))*
                        (FD[[i, 3]] - FD[[i - 1, 3]])], {i, 1, 
                    nc, 1}];		AjA = Table[
                  FD[[i, 4]] - 
                    (((1/3)/(1 - 1/3))*
                          (FD[[i, 1]] - FD[[i, 4]]) + 
                        FD[[i, 3]]), {i, 1, nc, 1}];			AjH = Table[
                  FD[[i, 5]] - 
                    (((1/3)/(1 - 1/3))*
                          (FD[[i, 1]] - FD[[i, 5]]) + 
                        FD[[i, 3]]), {i, 1, nc, 1}];			A = If[HJ == 1, 
                  Table[(∑_(j = 1)^i) dA[[j]], {i, 1, nc, 1}], 
                  Table[If[i == 1, 0, 
                      (FD[[i, 4]] - Mean[FD][[4]]) - 
                        0.5*(aA[[i]] + Mean[aA])*
                          (FD[[i, 1]] - Mean[FD][[1]]) - 
                        0.5*(bA[[i]] + Mean[bA])*
                          (FD[[i, 2]] - Mean[FD][[2]]) - 
                        0.5*(
                            (1 - aA[[i]] - bA[[i]]) + 
                              (1 - Mean[aA] - Mean[bA]))*
                          (FD[[i, 3]] - Mean[FD][[3]])], {i, 
                      1, nc, 1}]];		Ah = If[HJ == 1, 
                  Table[(∑_(j = 1)^i) dAh[[j]], {i, 1, nc, 1}], 
                  Table[If[i == 1, 0, 
                      (FD[[i, 5]] - Mean[FD][[5]]) - 
                        0.5*(aH[[i]] + Mean[aH])*
                          (FD[[i, 1]] - Mean[FD][[1]]) - 
                        0.5*(bH[[i]] + Mean[bH])*
                          (FD[[i, 2]] - Mean[FD][[2]]) - 
                        0.5*(
                            (1 - aH[[i]] - bH[[i]]) + 
                              (1 - Mean[aH] - Mean[bH]))*
                          (FD[[i, 3]] - Mean[FD][[3]])], {i, 
                      1, nc, 1}]];		AjA = Exp[AjA]/Exp[AjA[[121]]];
AjH = Exp[AjH]/Exp[AjH[[121]]];		o = Table[
                  If[j == 1, A[[i]], 
                    If[j == 2, Ah[[i]], 
                      If[j == 3, AjA[[i]], 
                        If[j == 4, AjH[[i]], 
                          If[j == 5, FD[[i, 6]], FD[[i, 4]]]]]]], {i, 
                    1, nc, 1}, {j, 1, 6, 1}]			
];

This function implements the λ* equilibrium assignment criterion discussed in the text.

GoodOrBad[rdlist_, nc_, γ_, α_, β_, 
            λ_, cθ_]:=Module[{λCrit, a, s},
a = Table[rdlist[[i, 1]]/100, {i, 1, nc, 1}];
s = Table[rdlist[[i, 2]]/100, {i, 1, nc, 1}];		λCrit=Table[(a[[i]]*(γ - α)+β*(1 - a[[i]])-(γ - α - β)*(α/γ)*(1 - a[[i]])*(((s[[i]])/(1 - s[[i]]))))*(((a[[i]]*(1 - γ)+α*(((s[[i]])/(1 - s[[i]])))*(1 + 
                                    a[[i]]*
                                      (((1 - γ)/γ)))
                                  ))^(-1)),{i, 1, nc, 1}];				Trap = 
                  Table[If[j == 1, 
                      If[λCrit[[i]] < λ, 1, 0], 
                      If[j == 2, λCrit[[i]], 
                        rdlist[[i, 8]]]], {i, 1, nc, 1}, {j, 1, 3, 1}]
];

A first look at the cross country data

Here we read in the core dataset used for analysis. The data augments the Hall and Jones (1999) dataset to include sectoral data on labor shares and value added. The data set includes the following variables: a88, s88adj, logYL, logKL, logHL, logL, logD, Region, bgcode, logA, SocInf, EngFrac, EurFrac, Latitude, EcOrg, YrsOpen, GADP, logFrankRomer, Mining, s88, pop88, and wbcode

\(rdlist = 
      ReadList["cdata.txt", 
        \ {Number, \ Number, Number, Number, Number, Number, Number, Number, 
          Number, Number, Number, Number, Number, Number, Number, Number, 
          Number, Number, Number, Number, Number, String}]; 

    Short[rdlist]\)nc = 127;

This section uses the Sheather and Jones (1991) "plug in" method to choose "optimal" bandwidth levels for density estimates of Hall and Jones' (1999) output per worker data. Note this routine is computationally intensive and takes a few minutes to complete.

\(JonesYL = Table[rdlist[\([i, 3]\)], {i, 1, nc, 1}]; 

    JonesKL = Table[rdlist[\([i, 4]\)], {i, 1, nc, 1}]; 

    JonesHL = Table[rdlist[\([i, 5]\)], {i, 1, nc, 1}]; 

    JonesL = Table[rdlist[\([i, 6]\)], {i, 1, nc, 1}]; 

    JonesTFP = Table[rdlist[\([i, 10]\)], {i, 1, nc, 1}]; 

    Land = Table[rdlist[\([i, 7]\)], {i, 1, nc, 1}]; 

    JonesDL = Land - JonesL; 

    na88\  = 1 - \ \((Table[rdlist[\([i, 1]\)], {i, 1, nc, 1}])\)/100; \)

hYL = KDenBandPISheather[JonesYL, nc, 0, 0, 15, 15]hKL = KDenBandPISheather[JonesKL, nc, 0, 0, 15, 15]hHL = KDenBandPISheather[JonesHL, nc, 0, 0, 15, 15]hTFP = KDenBandPISheather[JonesTFP, nc, 0, 0, 15, 15]

This section plots simple kernel density estimates for the TFP, GDP, and capital stock (human and physical) data used in the Hall and Jones (1999) growth accounting exercise.

\(FigYL := 
      GraphKernelEstimate[
        KDensity[JonesYL, nc, hYL, 0, Min[JonesYL] - 1, Max[JonesYL] + 1, 
          250, 0.05], {}, 250, 1, 0, "Log of output per worker", 
        "f[y] - kernel density estimate", "", \ 0.0, 6.5, 11, 0, 
        0.5]; 
FigKL := 
      GraphKernelEstimate[
        KDensity[JonesKL, nc, hKL, 0, Min[JonesKL] - 1, Max[JonesKL] + 1, 
          250, 0.05], {}, 250, 1, 0, 
        "Log of physical capital per worker", 
        "f[k] - kernel density estimate", "", \ 0.0, 5.5, 12.0, 0, 
        0.4]\)FigHL := 
      GraphKernelEstimate[
        KDensity[JonesHL, nc, hHL, 0, Min[JonesHL] - 1, Max[JonesHL] + 1, 
          250, 0.05], {}, 250, 1, 0, "Log of human capital per worker", 
        "f[h] - kernel density estimate", "",  0.0, 0, 1.2, 0, 1.75]
        FigHJTFP := 
      GraphKernelEstimate[
        KDensity[JonesTFP, nc, hTFP, 0, Min[JonesTFP] - 1, Max[JonesTFP] + 1, 
          250, 0.05], {}, 250, 1, 0, 
        "Log of Hall-Jones total factor productivity", 
        "f[A] - kernel density estimate", "",  0.0, 6.0, 9.5, 0, 0.8]

Here's a simple density plot for the distribution of labor across the agricultural and non-agricultural sectors.

hna88 = KDenBandPISheather[na88, nc, 0, 0, 15, 15]

Figna88 := 
      GraphKernelEstimate[KDensity[na88, nc, hna88, 0, 0, 1, 250, 0.05], {}, 
        250, 1, 0, "Share of workers in the non-agricultural sector", 
        "f[1-a] - kernel density estimate", "",  0.0, 0, 1, 0, 2.25]

Show[GraphicsArray[{{FigYL, FigKL}, {Figna88, FigHL}}]]

Aggregate levels accounting and multiple equilibria (Section 7, Table 3 & Figure 4)

This section explores the implications of multiple equilibria for interpretating aggregate growth accounting exercises. We use the standard Christenson, Cummings and Jorgenson (1980) method for constructing aggregate TFP estimates by country and we then try a see how these estimates are sensitive to the presence of multiple equilibria in the data. This code reproduces the material used in Section 7 of the paper. For each case the "1" to "5" suffixes correspond to the cases of λ=0.05 to λ=0.5.

\(ATFP1 = 
      AggegrateTFP[rdlist, nc, 0.35, 0.4, 0.2, 0.05, 0, JonesYL, JonesKL, 
        JonesHL, JonesDL, JonesTFP, 0]; 

    ATFP2 = AggegrateTFP[rdlist, nc, 0.35, 0.4, 0.2, 0.10, 0, JonesYL, 
        JonesKL, JonesHL, JonesDL, JonesTFP, 0]; 

    ATFP3 = AggegrateTFP[rdlist, nc, 0.35, 0.4, 0.2, 0.20, 0, JonesYL, 
        JonesKL, JonesHL, JonesDL, JonesTFP, 0]; 

    ATFP4 = AggegrateTFP[rdlist, nc, 0.35, 0.4, 0.2, 0.30, 0, JonesYL, 
        JonesKL, JonesHL, JonesDL, JonesTFP, 0]; 

    ATFP5 = AggegrateTFP[rdlist, nc, 0.35, 0.4, 0.2, 0.50, 0, JonesYL, 
        JonesKL, JonesHL, JonesDL, JonesTFP, 0]; \)

\(ALow = Table[ATFP1[\([i, 1]\)], {i, 1, nc, 1}]; 

    AHigh1 = Table[ATFP1[\([i, 2]\)], {i, 1, nc, 1}]; 

    AHigh2 = Table[ATFP2[\([i, 2]\)], {i, 1, nc, 1}]; 

    AHigh3 = Table[ATFP3[\([i, 2]\)], {i, 1, nc, 1}]; 

    AHigh4 = Table[ATFP4[\([i, 2]\)], {i, 1, nc, 1}]; 

    AHigh5 = Table[ATFP5[\([i, 2]\)], {i, 1, nc, 1}]; \)

Calculate the fraction of the variation in the aggregate TFP estimates "explained" by multiple equilibria.  This material is reported in Table 3 of the paper.

ex1 = 1 - Variance[AHigh1]/Variance[ALow]ex2 = 1 - Variance[AHigh2]/Variance[ALow]ex3 = 1 - Variance[AHigh3]/Variance[ALow]ex4 = 1 - Variance[AHigh4]/Variance[ALow]ex5 = 1 - Variance[AHigh5]/Variance[ALow]

TableForm[DispersionReport[ALow]]TableForm[DispersionReport[AHigh1]]TableForm[DispersionReport[AHigh2]]TableForm[DispersionReport[AHigh3]]TableForm[DispersionReport[AHigh4]]TableForm[DispersionReport[AHigh5]]

Plot the aggregate TFP distribution along with several counterfactuals. This generates Figure 4 of the paper.

hALow = KDenBandPISheather[ALow, nc, 0, 0, 15, 15]hAHigh1 = KDenBandPISheather[AHigh1, nc, 0, 0, 15, 15]hAHigh2 = KDenBandPISheather[AHigh2, nc, 0, 0, 15, 15]hAHigh3 = KDenBandPISheather[AHigh3, nc, 0, 0, 15, 15]hAHigh4 = KDenBandPISheather[AHigh4, nc, 0, 0, 15, 15]hAHigh5 = KDenBandPISheather[AHigh5, nc, 0, 0, 15, 15]

\(DLow = 
      KDensity[ALow, nc, hALow, 0, Min[ALow] - 0.25, Max[ALow] + 0.25, 250, 
        0.05]; 
DHigh1 = 
      KDensity[AHigh1, nc, hAHigh1, 0, Min[AHigh1] - 0.25, 
        Max[AHigh1] + 0.25, 250, 0.05]; 

    DHigh2 = KDensity[AHigh2, nc, hAHigh2, 0, Min[AHigh2] - 0.25, 
        Max[AHigh2] + 0.25, 250, 0.05]; 

    DHigh3 = KDensity[AHigh3, nc, hAHigh3, 0, Min[AHigh3] - 0.25, 
        Max[AHigh3] + 0.25, 250, 0.05]; 

    DHigh4 = KDensity[AHigh4, nc, hAHigh4, 0, Min[AHigh4] - 0.25, 
        Max[AHigh4] + 0.25, 250, 0.05]; 

    DHigh5 = KDensity[AHigh4, nc, hAHigh5, 0, Min[AHigh5] - 0.25, 
        Max[AHigh5] + 0.25, 250, 0.05]; \)

\(l1 = 
      Table[If[j == 1, DLow[\([i, 1]\)], DLow[\([i, 2]\)]], {i, 1, 250, 1}, {
          j, 1, 2, 1}]; 

    l2 = Table[
        If[j == 1, DHigh1[\([i, 1]\)], DHigh1[\([i, 2]\)]], {i, 1, 250, 1}, {
          j, 1, 2, 1}]; 

    l3 = Table[
        If[j == 1, DHigh2[\([i, 1]\)], DHigh2[\([i, 2]\)]], {i, 1, 250, 1}, {
          j, 1, 2, 1}]; 

    l4 = Table[
        If[j == 1, DHigh3[\([i, 1]\)], DHigh3[\([i, 2]\)]], {i, 1, 250, 1}, {
          j, 1, 2, 1}]; 

    l5 = Table[
        If[j == 1, DHigh4[\([i, 1]\)], DHigh4[\([i, 2]\)]], {i, 1, 250, 1}, {
          j, 1, 2, 1}]; 

    l6 = Table[
        If[j == 1, DHigh5[\([i, 1]\)], DHigh5[\([i, 2]\)]], {i, 1, 250, 1}, {
          j, 1, 2, 1}]; \)

ATFPFig1 = 
      MultipleListPlot[l1, l3, l5, PlotJoined -> True, 
        Frame -> True, 
        FrameLabel 
          -> {
            StyleForm["Log TFP minus world average", FontSize -> 16], 
            StyleForm["Estimated density function", FontSize -> 16], 
            "", ""}, Axes -> False, 
        SymbolShape -> {None, None}, 
        PlotStyle 
          -> {{GrayLevel[0], Thickness[0.0075]}, {
              Dashing[{0.0125, 0.0125}], GrayLevel[0.5], Thickness[0.0075]}, {
              Dashing[{0.025, 0.025}], GrayLevel[0.5], Thickness[0.0075]}, {
              Dashing[{0.025, 0.0125}], GrayLevel[0.5], Thickness[0.0075]}, {
              Dashing[{0.025, 0.0125}], GrayLevel[0.5], Thickness[0.0075]}}, 
        PlotRange -> {{-2.0, 2.0}, {0, 1.5}}];

Show[{ATFPFig1, 
        Graphics[
          Text[StyleForm["Observed", FontSize -> 16], {0.7, 0.60}, {
              -1, 0}]], 
        Graphics[
          Text[StyleForm["λ = 0.10", FontSize -> 16], {0.5, 
              1.0}, {-1, 0}]], 
        Graphics[
          Text[StyleForm["λ = 0.3", FontSize -> 16], {0.15, 
              1.2}, {-1, 0}]]}]

Some auxiliary accounting results (Section 9, Table 5)

This function produces some auxiliary calibrated related output, most of which was not used in the actual paper. However the function was used to produce the "return to capital across equilibria" table (Table 5 in Section 9) in the paper.

\(\(ShiftShareAnalysis[rdlist_, nc_, γ_, α_, β_, 
        λ_, cθ_, CR_] := 
      Module[{dlist, \ trap, a, b, s, v, VR, muL, muK}, 
        
dlist = rdlist/100; \ 
		
        trap\  = \ Table[CR[\([i, 1]\)], {i, 1, nc, 1}]; 
		
        VR = \ Table[Exp[CR[\([i, 3]\)] - CR[\([i, 2]\)]], {i, 1, nc, 1}]; 

        		a\  = \ Table[CR[\([i, 5]\)], {i, 1, nc, 1}]; 
		
        b\  = \ Table[CR[\([i, 6]\)], {i, 1, nc, 1}]; 
		
        x\  = \ Table[CR[\([i, 8]\)], {i, 1, nc, 1}]; 
		
        z\  = \ Table[CR[\([i, 9]\)], {i, 1, nc, 1}]; 
		
        s\  = \ Table[
            x[\([i]\)]*
              γ/
                \((x[\([i]\)]*γ + \((1 - x[\([i]\)])\)*α)\), {
              i, 1, nc, 1}]; 				
        v\  = \ Table[
            z[\([i]\)]*
              γ/
                \((z[\([i]\)]*γ + \((1 - z[\([i]\)])\)*α)\), {
              i, 1, nc, 1}]; 	
		
        muL\  = \ 
          Table[\((v[\([i]\)]/s[\([i]\)])\)*\((a[\([i]\)]/b[\([i]\)])\)*
              VR[\([i]\)], {i, 1, nc, 1}]; 					
        muK\  = \ 
          Table[\((v[\([i]\)]/s[\([i]\)])\)*\((x[\([i]\)]/z[\([i]\)])\)*
              VR[\([i]\)], {i, 1, nc, 1}]; 
		
        rL\  = Table[
            γ*
              \((\((1 - s[\([i]\)])\)*
                    Exp[CR[\([i, 2]\)] + rdlist[\([i, 6]\)]])\)/
                \((\((1 - x[\([i]\)])\)*
                    Exp[rdlist[\([i, 4]\)] + rdlist[\([i, 6]\)]])\), {i, 1, 
              nc, 1}]\ ; 
		
        rH\  = Table[
            γ*
              \((\((1 - v[\([i]\)])\)*
                    Exp[CR[\([i, 3]\)] + rdlist[\([i, 6]\)]])\)/
                \((\((1 - z[\([i]\)])\)*
                    Exp[rdlist[\([i, 4]\)] + rdlist[\([i, 6]\)]])\), {i, 1, 
              nc, 1}]\ ; 
		
        AgProdCon\  = \ 
          Table[\((a[\([i]\)] + \((b[\([i]\)] - a[\([i]\)])\)/2)\)*
              \((muL[\([i]\)] - 1)\)*\((s[\([i]\)]/a[\([i]\)])\), {i, 1, nc, 
              1}]; 
		
        NonAgProdCon\  = \ 
          Table[\((1 - \((a[\([i]\)] + \((b[\([i]\)] - a[\([i]\)])\)/2)\))\)*
              \((muL[\([i]\)] - 1)\)*
              \((\((1 - s[\([i]\)])\)/\((1 - a[\([i]\)])\))\), {i, 1, nc, 
              1}]; 
		
        RA1 = Table[
            \((\((b[\([i]\)] - a[\([i]\)])\) + 
                  \((\((b[\([i]\)] - a[\([i]\)])\) \((muL[\([i]\)] - 1)\))\)/
                    2)\)*\((
                s[\([i]\)]/a[\([i]\)] - 
                  \((1 - s[\([i]\)])\)/\((1 - a[\([i]\)])\))\), {i, 1, nc, 
              1}]; 
		
        SSlist = 
          Table[If[j == 2, muK[\([i]\)], 
              If[j == 3, rL[\([i]\)], 
                If[j == 4, rH[\([i]\)], 
                  If[j == 5, VR[\([i]\)] - 1, 
                    If[j == 6, AgProdCon[\([i]\)], 
                      If[j == 7, NonAgProdCon[\([i]\)], 
                        If[j == 8, RA1[\([i]\)], trap[\([i]\)]]]]]]]], {i, 1, 
              nc, 1}, {j, 1, 8, 1}]]; \)\)

\(CR1 = FullCalibration[rdlist, nc, 0.35, 0.4, 0.2, 0.05, 0]; 

    CR2 = FullCalibration[rdlist, nc, 0.35, 0.4, 0.2, 0.1, 0]; 

    CR3 = FullCalibration[rdlist, nc, 0.35, 0.4, 0.2, 0.2, 0]; 

    CR4 = FullCalibration[rdlist, nc, 0.35, 0.4, 0.2, 0.3, 0]; 

    CR5 = FullCalibration[rdlist, nc, 0.35, 0.4, 0.2, 0.5, 0]; \)

\(SSR1 = ShiftShareAnalysis[rdlist, nc, 0.35, 0.4, 0.2, 0.05, 0, CR1]; 

    SSR2 = ShiftShareAnalysis[rdlist, nc, 0.35, 0.4, 0.2, 0.1, 0, CR2]; 

    SSR3 = ShiftShareAnalysis[rdlist, nc, 0.35, 0.4, 0.2, 0.2, 0, CR3]; 

    SSR4 = ShiftShareAnalysis[rdlist, nc, 0.35, 0.4, 0.2, 0.3, 0, CR4]; 

    SSR5 = ShiftShareAnalysis[rdlist, nc, 0.35, 0.4, 0.2, 0.5, 0, CR5]; \)

\(rAct\  = \ 
      Log[Table[
          SSR1[\([i, 1]\)]*SSR1[\([i, 3]\)] + 
            \((1 - SSR1[\([i, 1]\)])\)*SSR1[\([i, 4]\)], {i, 1, nc, 1}]]; 

    rH1 = \ Log[Table[SSR1[\([i, 4]\)], {i, 1, nc, 1}]]; 

    rH2 = \ Log[Table[SSR2[\([i, 4]\)], {i, 1, nc, 1}]]; 

    rH3 = \ Log[Table[SSR3[\([i, 4]\)], {i, 1, nc, 1}]]; 

    rH4 = \ Log[Table[SSR4[\([i, 4]\)], {i, 1, nc, 1}]]; 

    rH5 = \ Log[Table[SSR5[\([i, 4]\)], {i, 1, nc, 1}]]; \)

This figure shows how the distribution of the return to capital is affected by the prevalence of poverty traps. This figure was not reported in the paper.

\(DrAct = CummDenF[Sort[rAct], nc]; 
DrH1 = CummDenF[Sort[rH1], nc]; 

    DrH2 = CummDenF[Sort[rH2], nc]; 
DrH3 = CummDenF[Sort[rH3], nc]; 

    DrH4 = CummDenF[Sort[rH4], nc]; 
DrH5 = CummDenF[Sort[rH5], nc]; \)

\(l1 = 
      Table[If[j == 1, DrAct[\([i, 1]\)], DrAct[\([i, 2]\)]], {i, 1, nc, 1}, {
          j, 1, 2, 1}]; 

    l2 = Table[
        If[j == 1, DrH1[\([i, 1]\)], DrH1[\([i, 2]\)]], {i, 1, nc, 1}, {j, 1, 
          2, 1}]; 

    l3 = Table[
        If[j == 1, DrH2[\([i, 1]\)], DrH2[\([i, 2]\)]], {i, 1, nc, 1}, {j, 1, 
          2, 1}]; 

    l4 = Table[
        If[j == 1, DrH3[\([i, 1]\)], DrH3[\([i, 2]\)]], {i, 1, nc, 1}, {j, 1, 
          2, 1}]; \)

rFig1 = 
      MultipleListPlot[l1 , l2, l4, PlotJoined -> True, 
        Frame -> True, 
        FrameLabel 
          -> {
            StyleForm["Log of real interest rate", FontSize -> 16], 
            StyleForm["CDF", FontSize -> 16], "", ""}, 
        Axes -> False, SymbolShape -> {None, None}, 
        PlotStyle 
          -> {{GrayLevel[0], Thickness[0.0075]}, {
              Dashing[{0.0125, 0.0125}], GrayLevel[0.5], Thickness[0.0075]}, {
              Dashing[{0.0125, 0.0125}], GrayLevel[0.5], Thickness[0.0075]}}, 
        PlotRange -> {{-3, 1}, {0, 1}}];

Show[{rFig1, 
        Graphics[
          Text[StyleForm["Observed", FontSize -> 16], {-1.7, 0.60}, {
              1, 0}]], 
        Graphics[
          Text[StyleForm["λ = 0.2", FontSize -> 16], {-1, 
              0.70}, {-1, 0}]], 
        Graphics[
          Text[StyleForm["λ = 0.05", FontSize -> 16], {-0.1, 
              0.90}, {-1, 0}]]}]

These commands produce the numbers reported in Table 5 of the paper.

NumTrap1 = Sum[SSR1[[i, 1]], {i, 1, nc, 1}]Sum[SSR1[[i, 2]]*SSR1[[i, 1]], {i, 1, nc, 1}]/NumTrap1Sum[SSR1[[i, 2]]*(1 - SSR1[[i, 1]]), {i, 1, nc, 1}]/
      (nc - NumTrap1)\(rC = Sort[SSR1]; 

    Median[Table[rC[\([i, 2]\)], {i, nc - NumTrap1 + 1, nc, 1}]]\)\(Median[Table[rC[\([i, 2]\)], {i, 1, nc - NumTrap1, 1}]]
\)

NumTrap2 = Sum[SSR2[[i, 1]], {i, 1, nc, 1}]Sum[SSR2[[i, 2]]*SSR2[[i, 1]], {i, 1, nc, 1}]/NumTrap2Sum[SSR2[[i, 2]]*(1 - SSR2[[i, 1]]), {i, 1, nc, 1}]/
      (nc - NumTrap2)\(rC = Sort[SSR2]; 

    Median[Table[rC[\([i, 2]\)], {i, nc - NumTrap2 + 1, nc, 1}]]\)Median[Table[rC[[i, 2]], {i, 1, nc - NumTrap2, 1}]]

NumTrap3 = Sum[SSR3[[i, 1]], {i, 1, nc, 1}]Sum[SSR3[[i, 2]]*SSR3[[i, 1]], {i, 1, nc, 1}]/NumTrap3Sum[SSR3[[i, 2]]*(1 - SSR3[[i, 1]]), {i, 1, nc, 1}]/
      (nc - NumTrap3)\(rC = Sort[SSR3]; 

    Median[Table[rC[\([i, 2]\)], {i, nc - NumTrap3 + 1, nc, 1}]]\)Median[Table[rC[[i, 2]], {i, 1, nc - NumTrap3, 1}]]

NumTrap4 = Sum[SSR4[[i, 1]], {i, 1, nc, 1}]Sum[SSR4[[i, 2]]*SSR4[[i, 1]], {i, 1, nc, 1}]/NumTrap4Sum[SSR4[[i, 2]]*(1 - SSR4[[i, 1]]), {i, 1, nc, 1}]/
      (nc - NumTrap4)\(rC = Sort[SSR4]; 

    Median[Table[rC[\([i, 2]\)], {i, nc - NumTrap4 + 1, nc, 1}]]\)Median[Table[rC[[i, 2]], {i, 1, nc - NumTrap4, 1}]]

NumTrap5 = Sum[SSR5[[i, 1]], {i, 1, nc, 1}]Sum[SSR5[[i, 2]]*SSR5[[i, 1]], {i, 1, nc, 1}]/NumTrap5Sum[SSR5[[i, 2]]*(1 - SSR5[[i, 1]]), {i, 1, nc, 1}]/
      (nc - NumTrap5)\(rC = Sort[SSR5]; 

    Median[Table[rC[\([i, 2]\)], {i, nc - NumTrap5 + 1, nc, 1}]]\)Median[Table[rC[[i, 2]], {i, 1, nc - NumTrap5, 1}]]

Robustness of results to parameter choices for α, β, γ and λ (Section 8, Figure 5)

Two key parameters in our model are our choices for α -- the capital share in the agricultural sector -- and our choice for λ. This section evaluates the sensitivity of our results to these and our other parameter choices. This material forms the basis of Figure 5 in Section 8 of the paper.

$TextStyle  = 
       {FontFamily  ->  "Times",  FontSize  ->  16};

<< Graphics`Graphics3D`

Robustness over λ and α

\(NS = 10; 
ls = 0.0; 
as = 0.0; 
lstep = 0.05; 
astep = 0.05; 

    r = Table[0, {h, 1, NS, 1}, {i, 1, NS, 1}, {j, 1, 3, 1}]; \)

\(\(Do[
      \(Do[
		
          wlist = GoodOrBad[rdlist, nc, 0.35, as + q*astep, 0.2, \ 
              ls + j*lstep, 0]; 
		
          PercentTrap = 
            N[Sum[\(Table[wlist[\([i, 1]\)], {i, 1, nc, 1}]\)[\([t]\)], {t, 
                    1, nc, 1}]/nc, 3]; 
		r[\([q, j, 1]\)] = ls + j*lstep; 
          
		r[\([q, j, 2]\)] = as + q*astep; 
		
          \(r[\([q, j, 3]\)] = PercentTrap;\), {j, 1, NS, 1}];\), 
      \ {q, 1, NS, 1}]; \)\)

RFigλα = 
      ListSurfacePlot3D[r, Axes -> True, 
        AxesLabel -> {"λ", "α", ""}, 
        Boxed -> True, Mesh -> True, Lighting -> True,  
        BoxRatios  ->  {5, 4,  3}, 
        ViewPoint  ->  {-2,  -2,  1}]

Robustness over λ and β

\(NS = 10; 
ls = 0.0; 
bs = 0.0; 
lstep = 0.05; 
bstep = 0.05; 

    r = Table[0, {h, 1, NS, 1}, {i, 1, NS, 1}, {j, 1, 3, 1}]; \)

\(\(Do[
      \(Do[
wlist = 
          GoodOrBad[rdlist, nc, 0.35, 0.4, bs + q*bstep, \ ls + j*lstep, 0]; 
        PercentTrap = 
          N[Sum[\(Table[wlist[\([i, 1]\)], {i, 1, nc, 1}]\)[\([t]\)], {t, 1, 
                  nc, 1}]/nc, 3]; \ 		r[\([q, j, 1]\)] = ls + j*lstep; 
        
r[\([q, j, 2]\)] = bs + q*bstep; 
        
\(r[\([q, j, 3]\)] = PercentTrap; 
        \), {j, 1, NS, 1}]; \), \ {q, 1, NS, 1}]; \)\)

RFigλβ = 
      ListSurfacePlot3D[r, Axes -> True, 
        AxesLabel -> {"λ", "β", ""}, Boxed -> True, 
        Mesh -> True, Lighting -> True,  BoxRatios  ->  {5, 4,  3}, 
        ViewPoint  ->  {-6,  -2,  1.5}]

Robustness over λ and γ

\(NS = 12; 
ls = 0.0; 
gs = 0.0; 
lstep = 0.05; 
gstep = 0.05; 

    r = Table[0, {h, 1, NS, 1}, {i, 1, NS, 1}, {j, 1, 3, 1}]; \)

\(\(Do[
      \(Do[
wlist = 
          GoodOrBad[rdlist, nc, gs + q*gstep, 0.4, 0.2, \ ls + j*lstep, 0]; 
        PercentTrap = 
          N[Sum[\(Table[wlist[\([i, 1]\)], {i, 1, nc, 1}]\)[\([t]\)], {t, 1, 
                  nc, 1}]/nc, 3]; \ 		r[\([q, j, 1]\)] = ls + j*lstep; 
        
r[\([q, j, 2]\)] = gs + q*gstep; 
        
\(r[\([q, j, 3]\)] = PercentTrap; 
        \), {j, 1, NS, 1}]; \), \ {q, 1, NS, 1}]; \)\)

RFigλγ = 
      ListSurfacePlot3D[r, Axes -> True, 
        AxesLabel -> {"λ", "γ", ""}, 
        Boxed -> True, Mesh -> True, Lighting -> True,  
        BoxRatios  ->  {5, 4,  3}, 
        ViewPoint  ->  {-6,  -1,  1.5}]

Robustness over γ and α

\(NS = 12; 
as = 0.0; 
gs = 0.0; 
astep = 0.05; 
gstep = 0.05; 

    r = Table[0, {h, 1, NS, 1}, {i, 1, NS, 1}, {j, 1, 3, 1}]; \)

\(\(Do[
      \(Do[
wlist = 
          GoodOrBad[rdlist, nc, gs + q*gstep, as + j*astep, 0.2, \ 0.3, 0]; 
        PercentTrap = 
          N[Sum[\(Table[wlist[\([i, 1]\)], {i, 1, nc, 1}]\)[\([t]\)], {t, 1, 
                  nc, 1}]/nc, 3]; \ 		r[\([q, j, 1]\)] = as + j*astep; 
        
r[\([q, j, 2]\)] = gs + q*gstep; 
        
\(r[\([q, j, 3]\)] = PercentTrap; 
        \), {j, 1, NS, 1}]; \), \ {q, 1, NS, 1}]; \)\)

RFigγα = 
      ListSurfacePlot3D[r, Axes -> True, 
        AxesLabel -> {"α", "γ", ""}, Boxed -> True, 
        Mesh -> True, Lighting -> True,  BoxRatios  ->  {5, 4,  3}, 
        ViewPoint  ->  {-6,  -2,  1.5}]

Summary panel of figures (Figure 5 in the paper)

Show[
      GraphicsArray[{{RFigλα, RFigλβ}, {
            RFigλγ, RFigγα}}]]

Robustness to mis-measured agricultural output (Section 8, Table 4)

This section assesses the extent to which our main results are robust to mis-measured agricultural output as reported in Section 8 of the paper (Table 4).

sOrginal  =  Table[rdlist[[i, 2]], {i, 1, nc, 1}];

sAdjustedData[err_] :=  
      Table[If[j == 2, 
          100*rdlist[[i, 2]]/(100*err + 
                rdlist[[i, 2]]*(1 - err)), 
          If[j == 3, 
            Log[(((rdlist[[i, 2]]/100)*
                      Exp[rdlist[[i, 3]] + rdlist[[i, 6]]])/err) + 
                  (((100 - rdlist[[i, 2]])/100))*
                    Exp[rdlist[[i, 3]] + rdlist[[i, 6]]]] - 
              rdlist[[i, 6]], rdlist[[i, j]]]], {i, 1, nc, 1}, {j, 1, 
          Length[rdlist[[1]]], 1}];

\(Data90 = sAdjustedData[0.9]; 
Data80 = sAdjustedData[0.8]; 

    Data70 = sAdjustedData[0.7]; 
Data60 = sAdjustedData[0.6]; 

    Data50 = sAdjustedData[0.5]; \)

\(wlist1 = WelfareCalibration[Data90, nc, 0.35, 0.4, 0.2, 0.10, 0, 22]; 

    wlist2 = WelfareCalibration[Data80, nc, 0.35, 0.4, 0.2, 0.10, 0, 22]; 

    wlist3 = WelfareCalibration[Data70, nc, 0.35, 0.4, 0.2, 0.10, 0, 22]; 

    wlist4 = WelfareCalibration[Data60, nc, 0.35, 0.4, 0.2, 0.10, 0, 22]; 

    wlist5 = WelfareCalibration[Data50, nc, 0.35, 0.4, 0.2, 0.10, 0, 22]; 
    \)

\(GDPact1 = Table[Data90[\([i, 3]\)], {i, 1, nc, 1}]; 

    GDPact2 = Table[Data80[\([i, 3]\)], {i, 1, nc, 1}]; 

    GDPact3 = Table[Data70[\([i, 3]\)], {i, 1, nc, 1}]; 

    GDPact4 = Table[Data60[\([i, 3]\)], {i, 1, nc, 1}]; 

    GDPact5 = Table[Data50[\([i, 3]\)], {i, 1, nc, 1}]; \)

\(GDPhigh1 = 
      Table[If[wlist1[\([i, 2]\)] <= wlist1[\([i, 1]\)], Data90[\([i, 3]\)], 
          wlist1[\([i, 2]\)]], {i, 1, nc, 1}]; 

    GDPlow1 = 
      Table[If[wlist1[\([i, 2]\)] > wlist1[\([i, 1]\)], Data90[\([i, 3]\)], 
          wlist1[\([i, 2]\)]], {i, 1, nc, 1}]; 

    GDPhigh2 = 
      Table[If[wlist2[\([i, 2]\)] <= wlist2[\([i, 1]\)], Data80[\([i, 3]\)], 
          wlist2[\([i, 2]\)]], {i, 1, nc, 1}]; 

    GDPlow2 = 
      Table[If[wlist2[\([i, 2]\)] > wlist2[\([i, 1]\)], Data80[\([i, 3]\)], 
          wlist2[\([i, 2]\)]], {i, 1, nc, 1}]; 

    GDPhigh3 = 
      Table[If[wlist3[\([i, 2]\)] <= wlist3[\([i, 1]\)], Data70[\([i, 3]\)], 
          wlist3[\([i, 2]\)]], {i, 1, nc, 1}]; 

    GDPlow3 = 
      Table[If[wlist3[\([i, 2]\)] > wlist3[\([i, 1]\)], Data70[\([i, 3]\)], 
          wlist3[\([i, 2]\)]], {i, 1, nc, 1}]; 

    GDPhigh4 = 
      Table[If[wlist4[\([i, 2]\)] <= wlist4[\([i, 1]\)], Data60[\([i, 3]\)], 
          wlist4[\([i, 2]\)]], {i, 1, nc, 1}]; 

    GDPlow4 = 
      Table[If[wlist4[\([i, 2]\)] > wlist4[\([i, 1]\)], Data60[\([i, 3]\)], 
          wlist4[\([i, 2]\)]], {i, 1, nc, 1}]; 

    GDPhigh5 = 
      Table[If[wlist5[\([i, 2]\)] <= wlist5[\([i, 1]\)], Data50[\([i, 3]\)], 
          wlist5[\([i, 2]\)]], {i, 1, nc, 1}]; 

    GDPlow5 = 
      Table[If[wlist5[\([i, 2]\)] > wlist5[\([i, 1]\)], Data50[\([i, 3]\)], 
          wlist5[\([i, 2]\)]], {i, 1, nc, 1}]; \)

Calculations reported in Column 1 of Table 4.

1 - Variance[GDPhigh1]/Variance[GDPact1]\(	1 - Variance[GDPhigh2]/Variance[GDPact2]\)\(	1 - Variance[GDPhigh3]/Variance[GDPact3]\)\(	1 - Variance[GDPhigh4]/Variance[GDPact4]\)\(	1 - Variance[GDPhigh5]/Variance[GDPact5]\)

Calculations reported in Columns 2 to 4 of Table 4.

\(PercentGain1 = 
      Table[If[j == 2, Exp[GDPhigh1[\([i]\)] - GDPlow1[\([i]\)]], 
          If[wlist1[\([i, 2]\)] <= wlist1[\([i, 1]\)], 0, 1]], {i, 1, nc, 
          1}, {j, 1, 2, 1}]; 

    NumTrap1 = Sum[PercentGain1[\([i, 1]\)], {i, 1, nc, 1}]\)nc - NumTrap1N[NumTrap1/nc]N[1 - NumTrap1/nc]PercentGainLow1 = 
      Sum[PercentGain1[[i, 2]]*PercentGain1[[i, 1]], {i, 1, nc, 1}]/
        NumTrap1PercentGainHigh1 = 
      Sum[PercentGain1[[i, 2]]*(1 - PercentGain1[[i, 1]]), {i, 1, 
            nc, 1}]/(nc - NumTrap1)\(PSG1 = Sort[PercentGain1]; 

    Median[Table[PSG1[\([i, 2]\)], {i, nc - NumTrap1 + 1, nc, 1}]]\)Median[Table[PSG1[[i, 2]], {i, 1, nc - NumTrap1 + 1, 1}]]

\(PercentGain2 = 
      Table[If[j == 2, Exp[GDPhigh2[\([i]\)] - GDPlow2[\([i]\)]], 
          If[wlist2[\([i, 2]\)] <= wlist2[\([i, 1]\)], 0, 1]], {i, 1, nc, 
          1}, {j, 1, 2, 1}]; 

    NumTrap2 = Sum[PercentGain2[\([i, 1]\)], {i, 1, nc, 1}]\)nc - NumTrap2N[NumTrap2/nc]N[1 - NumTrap2/nc]PercentGainLow2 = 
      Sum[PercentGain2[[i, 2]]*PercentGain2[[i, 1]], {i, 1, nc, 1}]/
        NumTrap2PercentGainHigh2 = 
      Sum[PercentGain2[[i, 2]]*(1 - PercentGain2[[i, 1]]), {i, 1, 
            nc, 1}]/(nc - NumTrap2)\(PSG2 = Sort[PercentGain2]; 

    Median[Table[PSG2[\([i, 2]\)], {i, nc - NumTrap2 + 1, nc, 1}]]\)Median[Table[PSG2[[i, 2]], {i, 1, nc - NumTrap2 + 1, 1}]]

\(PercentGain3 = 
      Table[If[j == 2, Exp[GDPhigh3[\([i]\)] - GDPlow3[\([i]\)]], 
          If[wlist3[\([i, 2]\)] <= wlist3[\([i, 1]\)], 0, 1]], {i, 1, nc, 
          1}, {j, 1, 2, 1}]; 

    NumTrap3 = Sum[PercentGain3[\([i, 1]\)], {i, 1, nc, 1}]\)nc - NumTrap3N[NumTrap3/nc]N[1 - NumTrap3/nc]PercentGainLow3 = 
      Sum[PercentGain3[[i, 2]]*PercentGain3[[i, 1]], {i, 1, nc, 1}]/
        NumTrap3PercentGainHigh3 = 
      Sum[PercentGain3[[i, 2]]*(1 - PercentGain3[[i, 1]]), {i, 1, 
            nc, 1}]/(nc - NumTrap3)\(PSG3 = Sort[PercentGain3]; 

    Median[Table[PSG3[\([i, 2]\)], {i, nc - NumTrap3 + 1, nc, 1}]]\)Median[Table[PSG3[[i, 2]], {i, 1, nc - NumTrap3 + 1, 1}]]

\(PercentGain4 = 
      Table[If[j == 2, Exp[GDPhigh4[\([i]\)] - GDPlow4[\([i]\)]], 
          If[wlist4[\([i, 2]\)] <= wlist4[\([i, 1]\)], 0, 1]], {i, 1, nc, 
          1}, {j, 1, 2, 1}]; 

    NumTrap4 = Sum[PercentGain4[\([i, 1]\)], {i, 1, nc, 1}]\)nc - NumTrap4N[NumTrap4/nc]N[1 - NumTrap4/nc]PercentGainLow4 = 
      Sum[PercentGain4[[i, 2]]*PercentGain4[[i, 1]], {i, 1, nc, 1}]/
        NumTrap4PercentGainHigh4 = 
      Sum[PercentGain4[[i, 2]]*(1 - PercentGain4[[i, 1]]), {i, 1, 
            nc, 1}]/(nc - NumTrap4)\(PSG4 = Sort[PercentGain4]; 

    Median[Table[PSG4[\([i, 2]\)], {i, nc - NumTrap4 + 1, nc, 1}]]\)Median[Table[PSG4[[i, 2]], {i, 1, nc - NumTrap4 + 1, 1}]]

\(PercentGain5 = 
      Table[If[j == 2, Exp[GDPhigh5[\([i]\)] - GDPlow5[\([i]\)]], 
          If[wlist5[\([i, 2]\)] <= wlist5[\([i, 1]\)], 0, 1]], {i, 1, nc, 
          1}, {j, 1, 2, 1}]; 

    NumTrap5 = Sum[PercentGain5[\([i, 1]\)], {i, 1, nc, 1}]\)nc - NumTrap5N[NumTrap5/nc]N[1 - NumTrap5/nc]PercentGainLow5 = 
      Sum[PercentGain5[[i, 2]]*PercentGain5[[i, 1]], {i, 1, nc, 1}]/
        NumTrap5PercentGainHigh5 = 
      Sum[PercentGain5[[i, 2]]*(1 - PercentGain5[[i, 1]]), {i, 1, 
            nc, 1}]/(nc - NumTrap5)\(PSG5 = Sort[PercentGain5]; 

    Median[Table[PSG5[\([i, 2]\)], {i, nc - NumTrap5 + 1, nc, 1}]]\)Median[Table[PSG5[[i, 2]], {i, 1, nc - NumTrap5 + 1, 1}]]

Multiple equilibria and the world distribution of income (Section 6, Tables 1 & 2, Figures 2 & 3)

This series of commands solves for output levels across different equilibria under different assumptions for λ. This allows us to explore the effects of equilibrium switching on the world distribution of income. The calculations done below formed the basis of Section 6 in the paper.

\(wlist1 = WelfareCalibration[rdlist, nc, 0.35, 0.4, 0.2, 0.05, 0, 22]; 

    wlist2 = WelfareCalibration[rdlist, nc, 0.35, 0.4, 0.2, 0.10, 0, 22]; 

    wlist3 = WelfareCalibration[rdlist, nc, 0.35, 0.4, 0.2, 0.20, 0, 22]; 

    wlist4 = WelfareCalibration[rdlist, nc, 0.35, 0.4, 0.2, 0.30, 0, 22]; 

    wlist5 = WelfareCalibration[rdlist, nc, 0.35, 0.4, 0.2, 0.50, 0, 22]; 
    \)

GDPact = Table[rdlist[[i, 3]], {i, 1, nc, 1}];

\(GDPhigh1 = 
      Table[If[wlist1[\([i, 2]\)] <= wlist1[\([i, 1]\)], rdlist[\([i, 3]\)], 
          wlist1[\([i, 2]\)]], {i, 1, nc, 1}]; 

    GDPlow1 = 
      Table[If[wlist1[\([i, 2]\)] > wlist1[\([i, 1]\)], rdlist[\([i, 3]\)], 
          wlist1[\([i, 2]\)]], {i, 1, nc, 1}]; 

    GDPhigh2 = 
      Table[If[wlist2[\([i, 2]\)] <= wlist2[\([i, 1]\)], rdlist[\([i, 3]\)], 
          wlist2[\([i, 2]\)]], {i, 1, nc, 1}]; 

    GDPlow2 = 
      Table[If[wlist2[\([i, 2]\)] > wlist2[\([i, 1]\)], rdlist[\([i, 3]\)], 
          wlist2[\([i, 2]\)]], {i, 1, nc, 1}]; 

    GDPhigh3 = 
      Table[If[wlist3[\([i, 2]\)] <= wlist3[\([i, 1]\)], rdlist[\([i, 3]\)], 
          wlist3[\([i, 2]\)]], {i, 1, nc, 1}]; 

    GDPlow3 = 
      Table[If[wlist3[\([i, 2]\)] > wlist3[\([i, 1]\)], rdlist[\([i, 3]\)], 
          wlist3[\([i, 2]\)]], {i, 1, nc, 1}]; 

    GDPhigh4 = 
      Table[If[wlist4[\([i, 2]\)] <= wlist4[\([i, 1]\)], rdlist[\([i, 3]\)], 
          wlist4[\([i, 2]\)]], {i, 1, nc, 1}]; 

    GDPlow4 = 
      Table[If[wlist4[\([i, 2]\)] > wlist4[\([i, 1]\)], rdlist[\([i, 3]\)], 
          wlist4[\([i, 2]\)]], {i, 1, nc, 1}]; 

    GDPhigh5 = 
      Table[If[wlist5[\([i, 2]\)] <= wlist5[\([i, 1]\)], rdlist[\([i, 3]\)], 
          wlist5[\([i, 2]\)]], {i, 1, nc, 1}]; 

    GDPlow5 = 
      Table[If[wlist5[\([i, 2]\)] > wlist5[\([i, 1]\)], rdlist[\([i, 3]\)], 
          wlist5[\([i, 2]\)]], {i, 1, nc, 1}]; \)

These numbers of reported in Tables 1 and 2 of the paper. The include the inequality decompositions and calculations of the mean and median output ratios across low and high equilibria.

\(PercentGain1 = 
      Table[If[j == 2, Exp[GDPhigh1[\([i]\)] - GDPlow1[\([i]\)]], 
          If[wlist1[\([i, 2]\)] <= wlist1[\([i, 1]\)], 0, 1]], {i, 1, nc, 
          1}, {j, 1, 2, 1}]; 

    NumTrap1 = Sum[PercentGain1[\([i, 1]\)], {i, 1, nc, 1}]\)nc - NumTrap1N[NumTrap1/nc]N[1 - NumTrap1/nc]PercentGainLow1 = 
      Sum[PercentGain1[[i, 2]]*PercentGain1[[i, 1]], {i, 1, nc, 1}]/
        NumTrap1PercentGainHigh1 = 
      Sum[PercentGain1[[i, 2]]*(1 - PercentGain1[[i, 1]]), {i, 1, 
            nc, 1}]/(nc - NumTrap1)\(PSG1 = Sort[PercentGain1]; 

    Median[Table[PSG1[\([i, 2]\)], {i, nc - NumTrap1 + 1, nc, 1}]]\)Median[Table[PSG1[[i, 2]], {i, 1, nc - NumTrap1 + 1, 1}]]

TRAP1 =  
          Table[If[wlist1[[i,  1]]  <=  wlist1[[i,  2]], 1,  0], 
             {i,  1, nc,  1}];
DCEntropy[Exp[GDPact], nc, TRAP1, 2]DCMeanLogDev[Exp[GDPact], nc, TRAP1]DCTheil[Exp[GDPact], nc, TRAP1]

\(PercentGain2 = 
      Table[If[j == 2, Exp[GDPhigh2[\([i]\)] - GDPlow2[\([i]\)]], 
          If[wlist2[\([i, 2]\)] <= wlist2[\([i, 1]\)], 0, 1]], {i, 1, nc, 
          1}, {j, 1, 2, 1}]; 

    NumTrap2 = Sum[PercentGain2[\([i, 1]\)], {i, 1, nc, 1}]\)nc - NumTrap2N[NumTrap2/nc]N[1 - NumTrap2/nc]PercentGainLow2 = 
      Sum[PercentGain2[[i, 2]]*PercentGain2[[i, 1]], {i, 1, nc, 1}]/
        NumTrap2PercentGainHigh2 = 
      Sum[PercentGain2[[i, 2]]*(1 - PercentGain2[[i, 1]]), {i, 1, 
            nc, 1}]/(nc - NumTrap2)\(PSG2 = Sort[PercentGain2]; 

    Median[Table[PSG2[\([i, 2]\)], {i, nc - NumTrap2 + 1, nc, 1}]]\)Median[Table[PSG2[[i, 2]], {i, 1, nc - NumTrap2 + 1, 1}]]

TRAP2 =  
          Table[If[wlist2[[i,  1]]  <=  wlist2[[i,  2]], 1,  0], 
             {i,  1, nc,  1}];
DCEntropy[Exp[GDPact], nc, TRAP2, 2]DCMeanLogDev[Exp[GDPact], nc, TRAP2]DCTheil[Exp[GDPact], nc, TRAP2]

\(PercentGain3 = 
      Table[If[j == 2, Exp[GDPhigh3[\([i]\)] - GDPlow3[\([i]\)]], 
          If[wlist3[\([i, 2]\)] <= wlist3[\([i, 1]\)], 0, 1]], {i, 1, nc, 
          1}, {j, 1, 2, 1}]; 

    NumTrap3 = Sum[PercentGain3[\([i, 1]\)], {i, 1, nc, 1}]\)nc - NumTrap3N[NumTrap3/nc]N[1 - NumTrap3/nc]PercentGainLow3 = 
      Sum[PercentGain3[[i, 2]]*PercentGain3[[i, 1]], {i, 1, nc, 1}]/
        NumTrap3PercentGainHigh3 = 
      Sum[PercentGain3[[i, 2]]*(1 - PercentGain3[[i, 1]]), {i, 1, 
            nc, 1}]/(nc - NumTrap3)\(PSG3 = Sort[PercentGain3]; 

    Median[Table[PSG3[\([i, 2]\)], {i, nc - NumTrap3 + 1, nc, 1}]]\)Median[Table[PSG3[[i, 2]], {i, 1, nc - NumTrap3 + 1, 1}]]

TRAP3 =  
          Table[If[wlist3[[i,  1]]  <=  wlist3[[i,  2]], 1,  0], 
             {i,  1, nc,  1}];
DCEntropy[Exp[GDPact], nc, TRAP3, 2]DCMeanLogDev[Exp[GDPact], nc, TRAP3]DCTheil[Exp[GDPact], nc, TRAP3]

\(PercentGain4 = 
      Table[If[j == 2, Exp[GDPhigh4[\([i]\)] - GDPlow4[\([i]\)]], 
          If[wlist4[\([i, 2]\)] <= wlist4[\([i, 1]\)], 0, 1]], {i, 1, nc, 
          1}, {j, 1, 2, 1}]; 

    NumTrap4 = Sum[PercentGain4[\([i, 1]\)], {i, 1, nc, 1}]\)nc - NumTrap4N[NumTrap4/nc]N[1 - NumTrap4/nc]PercentGainLow4 = 
      Sum[PercentGain4[[i, 2]]*PercentGain4[[i, 1]], {i, 1, nc, 1}]/
        NumTrap4PercentGainHigh4 = 
      Sum[PercentGain4[[i, 2]]*(1 - PercentGain4[[i, 1]]), {i, 1, 
            nc, 1}]/(nc - NumTrap4)\(PSG4 = Sort[PercentGain4]; 

    Median[Table[PSG4[\([i, 2]\)], {i, nc - NumTrap4 + 1, nc, 1}]]\)Median[Table[PSG4[[i, 2]], {i, 1, nc - NumTrap4 + 1, 1}]]

TRAP4 =  
          Table[If[wlist4[[i,  1]]  <=  wlist4[[i,  2]], 1,  0], 
             {i,  1, nc,  1}];
DCEntropy[Exp[GDPact], nc, TRAP4, 2]DCMeanLogDev[Exp[GDPact], nc, TRAP4]DCTheil[Exp[GDPact], nc, TRAP4]

\(PercentGain5 = 
      Table[If[j == 2, Exp[GDPhigh5[\([i]\)] - GDPlow5[\([i]\)]], 
          If[wlist5[\([i, 2]\)] <= wlist5[\([i, 1]\)], 0, 1]], {i, 1, nc, 
          1}, {j, 1, 2, 1}]; 

    NumTrap5 = Sum[PercentGain5[\([i, 1]\)], {i, 1, nc, 1}]\)nc - NumTrap5N[NumTrap5/nc]N[1 - NumTrap5/nc]PercentGainLow5 = 
      Sum[PercentGain5[[i, 2]]*PercentGain5[[i, 1]], {i, 1, nc, 1}]/
        NumTrap5PercentGainHigh5 = 
      Sum[PercentGain5[[i, 2]]*(1 - PercentGain5[[i, 1]]), {i, 1, 
            nc, 1}]/(nc - NumTrap5)\(PSG5 = Sort[PercentGain5]; 

    Median[Table[PSG5[\([i, 2]\)], {i, nc - NumTrap5 + 1, nc, 1}]]\)Median[Table[PSG5[[i, 2]], {i, 1, nc - NumTrap5 + 1, 1}]]

TRAP5 =  
          Table[If[wlist5[[i,  1]]  <=  wlist5[[i,  2]], 1,  0], 
             {i,  1, nc,  1}];
DCEntropy[Exp[GDPact], nc, TRAP5, 2]DCMeanLogDev[Exp[GDPact], nc, TRAP5]DCTheil[Exp[GDPact], nc, TRAP5]

Graphs of how the CDF of the world distribution of income is affected by eliminating low output equilibria (not reported in the paper), and graphs of how the PDF of the world distribution is affected (reported as Figure 3 in the paper).

\(CDFGDPact = CummDenF[Sort[GDPact], nc]; 

    CDFGDPhigh1 = CummDenF[Sort[GDPhigh1], nc]; 

    CDFGDPhigh2 = CummDenF[Sort[GDPhigh2], nc]; 

    CDFGDPhigh3 = CummDenF[Sort[GDPhigh3], nc]; 

    CDFGDPhigh4 = CummDenF[Sort[GDPhigh4], nc]; 

    CDFGDPhigh5 = CummDenF[Sort[GDPhigh5], nc]; \)

GDPFig1 = 
      MultipleListPlot[CDFGDPact, CDFGDPhigh2, CDFGDPhigh3, 
        PlotJoined -> True, Frame -> True, 
        FrameLabel 
          -> {"Log of GDP per worker", 
            "Cummulative Density Function", "", ""}, 
        Axes -> False, SymbolShape -> {None, None, None}, 
        PlotStyle 
          -> {{GrayLevel[0], Thickness[0.0075]}, {
              Dashing[{0.025, 0.0125}], GrayLevel[0.5], Thickness[0.0075]}, {
              Dashing[{0.025, 0.0125}], GrayLevel[0.5], Thickness[0.0075]}, {
              Dashing[{0.025, 0.0125}], GrayLevel[1], Thickness[0.0075]}}, 
        PlotRange -> {{6.5, 11}, {0, 1}}]

Show[{GDPFig1, 
        Graphics[Text["λ = 0.0", {7.45, 0.20}, {1, 0}]], 
        Graphics[Text["λ = 0.1", {7.575, 0.10}, {1, 0}]], 
        Graphics[Text["λ = 0.2", {8.25, 0.20}, {-1, 0}]]}]

LocationReport[GDPact]DispersionReport[GDPact]

LocationReport[GDPhigh1]DispersionReport[GDPhigh1]

LocationReport[GDPhigh2]DispersionReport[GDPhigh2]

LocationReport[GDPhigh3]DispersionReport[GDPhigh3]

LocationReport[GDPhigh4]DispersionReport[GDPhigh4]

LocationReport[GDPhigh5]DispersionReport[GDPhigh5]

These commands calculate the percentage of the variation in the worldwide distribution of income that can be explained by multiple equilibria. The compute how much, say the standard deviation of the log of output per worker, would fall in percentage terms if all countries were in their "high" equilibrium. These numbers are reported in Column 3 of Table 1.

PercentExplained1 = 
      1 - (Variance[GDPhigh1])/(Variance[GDPact])PercentExplained2 = 
      1 - (Variance[GDPhigh2])/(Variance[GDPact])PercentExplained3 = 
      1 - (Variance[GDPhigh3])/(Variance[GDPact])PercentExplained4 = 
      1 - (Variance[GDPhigh4])/(Variance[GDPact])PercentExplained5 = 
      1 - (Variance[GDPhigh5])/(Variance[GDPact])

hAct = KDenBandPISheather[GDPact, nc, 0, 0, 15, 15]hOpt1 = KDenBandPISheather[GDPhigh1, nc, 0, 0, 15, 15]hOpt2 = KDenBandPISheather[GDPhigh2, nc, 0, 0, 15, 15]hOpt3 = KDenBandPISheather[GDPhigh3, nc, 0, 0, 15, 15]hOpt4 = KDenBandPISheather[GDPhigh4, nc, 0, 0, 15, 15]hOpt5 = KDenBandPISheather[GDPhigh5, nc, 0, 0, 15, 15]

\(Case1 = 
      KDensity[GDPact, nc, 0.33*\((hAct + hOpt2 + hOpt3)\), 0, 6.5, 11, 250, 
        0.05]; 
Case2 = 
      KDensity[GDPhigh2, nc, 0.33*\((hAct + hOpt2 + hOpt3)\), 0, 6.5, 11, 
        250, 0.05]; 

    Case3 = KDensity[GDPhigh3, nc, 0.33*\((hAct + hOpt2 + hOpt3)\), 0, 6.5, 
        11, 250, 0.05]; 

    Case4 = KDensity[GDPhigh4, nc, hOpt4, 0, 6.5, 11, 250, 0.05]; 

    c1 = Table[
        If[j == 1, Case1[\([i, 1]\)], Case1[\([i, 2]\)]], {i, 1, 250, 
          1}, {j, 1, 2, 1}]; 

    c2 = Table[
        If[j == 1, Case2[\([i, 1]\)], Case2[\([i, 2]\)]], {i, 1, 250, 
          1}, {j, 1, 2, 1}]; 

    c3 = Table[
        If[j == 1, Case3[\([i, 1]\)], Case3[\([i, 2]\)]], {i, 1, 250, 
          1}, {j, 1, 2, 1}]; 

    c3 = Table[
        If[j == 1, Case4[\([i, 1]\)], Case4[\([i, 2]\)]], {i, 1, 250, 
          1}, {j, 1, 2, 1}]; \)

\(FigureActualDistribution := 
      ListPlot[c1, PlotJoined -> True, Frame -> True, 
        Axes -> False, 
        FrameLabel 
          -> {"Log of GDP per worker", 
            "Kernel density estimate", "", ""}, 
        PlotRange -> {{6.5, 11}, {0, 0.7}}, 
        PlotStyle -> {Thickness[0.0075]}]; 

    FigureOptimisticDistribution2 := 
      ListPlot[c2, PlotJoined -> True, Frame -> True, 
        PlotStyle 
          -> {Dashing[{0.0125, 0.0125}], GrayLevel[0.3], 
            Thickness[0.0075]}, Axes -> False, 
        FrameLabel 
          -> {"Log of GDP per worker", 
            "Kernel density estimate", "", ""}, 
        PlotRange -> {{6.5, 11}, {0, 0.7}}]; 

    FigureOptimisticDistribution3 := 
      ListPlot[c3, PlotJoined -> True, 
        PlotStyle 
          -> {Dashing[{0.0125, 0.0125}], GrayLevel[0.6], 
            Thickness[0.0075]}, Axes -> False, Frame -> True, 
        FrameLabel 
          -> {"Log of GDP per worker", 
            "Kernel density estimate", "", ""}, 
        PlotRange -> {{6.5, 11}, {0, 0.7}}]; 

    FigureOptimisticDistribution4 := 
      ListPlot[c4, PlotJoined -> True, 
        PlotStyle 
          -> {Dashing[{0.0125, 0.0125}], GrayLevel[0.1], 
            Thickness[0.0075]}, Axes -> False, Frame -> True, 
        FrameLabel 
          -> {"Log of GDP per worker", 
            "Kernel density estimate", "", ""}, 
        PlotRange -> {{6.5, 11}, {0, 0.75}}]; \)

Show[{FigureActualDistribution, FigureOptimisticDistribution2, 
        FigureOptimisticDistribution3, 
        Graphics[Text["λ = 0.0", {7.0, 0.29}, {-1, 0}]], 
        Graphics[Text["λ = 0.1", {7.5, 0.34}, {-1, 0}]], 
        Graphics[Text["λ = 0.2", {8.5, 0.475}, {1, 0}]]}]

Counterfactual labor allocations

These calculations form the basis of Figure 2 in the paper.

\(FCRes1 = FullCalibration[rdlist, nc, 0.35, 0.4, 0.2, 0.0, 0]; 

    FCRes2 = FullCalibration[rdlist, nc, 0.35, 0.4, 0.2, 0.1, 0]; 

    FCRes3 = FullCalibration[rdlist, nc, 0.35, 0.4, 0.2, 0.2, 0]; \)

\(naAct = 1 - Table[FCRes1[\([i, 4]\)], {i, 1, nc, 1}]; 

    naOpt2 = 1 - Table[FCRes2[\([i, 6]\)], {i, 1, nc, 1}]; 

    naOpt3 = 1 - Table[FCRes3[\([i, 6]\)], {i, 1, nc, 1}]; \)

hnaAct = KDenBandPISheather[naAct, nc, 0, 0, 15, 15]hnaOpt2 = KDenBandPISheather[naOpt2, nc, 0, 0, 15, 15]hnaOpt3 = KDenBandPISheather[naOpt3, nc, 0, 0, 15, 15]

\(CaseAct = 
      KDensity[naAct, nc, 0.33*\((hnaAct + hnaOpt2 + hnaOpt3)\), 0, 0, 1, 
        250, 0.05]; 

    CaseOpt2 = 
      KDensity[naOpt2, nc, 0.33*\((hnaAct + hnaOpt2 + hnaOpt3)\), 0, 0, 1, 
        250, 0.05]; 

    CaseOpt3 = 
      KDensity[naOpt3, nc, 0.33*\((hnaAct + hnaOpt2 + hnaOpt3)\), 0, 0, 1, 
        250, 0.05]; \)

\(c1 = 
      Table[If[j == 1, CaseAct[\([i, 1]\)], CaseAct[\([i, 2]\)]], {i, 
          1, 250, 1}, {j, 1, 2, 1}]; 

    c2 = Table[
        If[j == 1, CaseOpt2[\([i, 1]\)], CaseOpt2[\([i, 2]\)]], {i, 1, 
          250, 1}, {j, 1, 2, 1}]; 

    c3 = Table[
        If[j == 1, CaseOpt3[\([i, 1]\)], CaseOpt3[\([i, 2]\)]], {i, 1, 
          250, 1}, {j, 1, 2, 1}]; \)

\(FigureActualAllocation := 
      ListPlot[c1, PlotJoined -> True, Frame -> True, 
        Axes -> False, 
        FrameLabel 
          -> {"Non-agricultural labor share", 
            "Kernel density estimate", "", ""}, 
        PlotRange -> {{0, 1}, {0, 3}}, 
        PlotStyle -> {Thickness[0.0075]}]; 

    FigureOptimisticAllocation2 := 
      ListPlot[c2, PlotJoined -> True, Frame -> True, 
        PlotStyle 
          -> {Dashing[{0.0125, 0.0125}], GrayLevel[0.3], 
            Thickness[0.0075]}, Axes -> False, 
        FrameLabel 
          -> {"Non-agricultural labor share", 
            "Kernel density estimate", "", ""}, 
        PlotRange -> {{0, 1}, {0, 3}}]; 

    FigureOptimisticAllocation3 := 
      ListPlot[c3, PlotJoined -> True, 
        PlotStyle 
          -> {Dashing[{0.0125, 0.0125}], GrayLevel[0.6], 
            Thickness[0.0075]}, Axes -> False, Frame -> True, 
        FrameLabel 
          -> {"Non-agricultural labor share", 
            "Kernel density estimate", "", ""}, 
        PlotRange -> {{0, 1}, {0, 3}}]; \)

Show[{FigureActualAllocation, FigureOptimisticAllocation2, 
        FigureOptimisticAllocation3, 
        Graphics[Text["λ = 0.0", {0.3, 1.25}, {1, 0}]], 
        Graphics[Text["λ = 0.1", {0.6, 2}, {1, 0}]], 
        Graphics[Text["λ = 0.2", {0.9, 2.5}, {1, 0}]]}]

The Evolution of the World Distribution of Income, 1965 to 1988 (Transition Matrix Analysis in CEPR Discussion Paper No. 3046)

This section uses the model to asses how the world distribution of income might evolve over time. By calibrating the model in 1965 and 1988 we can see which countries move between equilibria over the period. This allows us to compute a transition matrix and solve for the ergodic distribution of states over time. This material was not included in the final version of the paper but can be found in the CEPR discussion paper version (#3046).

\(c1965 = 
      ReadList["1965.txt", 
        \ {Number, \ Number, Number, Number, Number, Number, Number, 
          String}]; 

    c1988 = ReadList["1988.txt", 
        \ {Number, \ Number, Number, Number, Number, Number, Number, 
          String}]; 
nc = 121; \)

\(TableForm[c1965]; 
TableForm[c1988]; \)

step1965  =  
            WelfareCalibration[c1965,  nc,  0.35,  0.4,  0.2,  0.3,  0, 
               8]; 
step1988  =  
            WelfareCalibration[c1988,  nc,  0.35,  0.4,  0.2,  0.3,  0, 
               8];

\(LowHigh\  := \ 
\ \ \ \ \ \ 
          Table[If[step1965[\([i, \ 1]\)]\  < \ step1965[\([i, \ 2]\)], \ 

              \ \ \ \ \ \ \ \ \ \ 
              If[step1988[\([i, \ 1]\)]\  > \ step1988[\([i, \ 2]\)], \ 1, \ 
                0], \ 0], \ {i, \ 1, \ 
\ \ \ \ \ \ \ \ \ \ nc, \ 1}]\);
 \(LowLow\  := \ 
\ \ \ \ \ \ 
        Table[If[step1965[\([i, \ 1]\)]\  < \ step1965[\([i, \ 2]\)], \ 

            \ \ \ \ \ \ \ \ \ \ 
            If[step1988[\([i, \ 1]\)]\  < \ step1988[\([i, \ 2]\)], \ 1, \ 
              0], \ 0], \ {i, \ 1, \ 
\ \ \ \ \ \ \ \ \ \ nc, \ 1}]\);
\(HighLow\  := \ 
\ \ \ \ \ \ 
        Table[If[step1965[\([i, \ 1]\)]\  > \ step1965[\([i, \ 2]\)], \ 

            \ \ \ \ \ \ \ \ \ \ 
            If[step1988[\([i, \ 1]\)]\  < \ step1988[\([i, \ 2]\)], \ 1, \ 
              0], \ 0], \ {i, \ 1, \ 
\ \ \ \ \ \ \ \ \ \ nc, \ 1}]\);
\(HighHigh\  := \ 
\ \ \ \ \ \ 
        Table[If[step1965[\([i, \ 1]\)]\  > \ step1965[\([i, \ 2]\)], \ 

            \ \ \ \ \ \ \ \ \ \ 
            If[step1988[\([i, \ 1]\)]\  > \ step1988[\([i, \ 2]\)], \ 1, \ 
              0], \ 0], \ {i, \ 1, \ 
\ \ \ \ \ \ \ \ \ \ nc, \ 1}]\);

TableForm[
      Table[If[j == 1, step1988[[i, 3]], 
          If[j == 2, LowHigh[[i]], 
            If[j == 3, LowLow[[i]], HighHigh[[i]]]]], {i, 1, nc, 1}, {
          j, 1, 4, 1}]];

\(NL65 = 
      Sum[LowLow[\([i]\)], {i, 1, nc, 1}] + 
        Sum[LowHigh[\([i]\)], {i, 1, nc, 1}]; 

    NH65 = nc - 
        \((Sum[LowLow[\([i]\)], {i, 1, nc, 1}] + 
            Sum[LowHigh[\([i]\)], {i, 1, nc, 1}])\); 

    NL88 = Sum[LowLow[\([i]\)], {i, 1, nc, 1}] + 
        Sum[HighLow[\([i]\)], {i, 1, nc, 1}]; 

    NH88 = nc - 
        \((Sum[LowLow[\([i]\)], {i, 1, nc, 1}] + 
            Sum[HighLow[\([i]\)], {i, 1, nc, 1}])\); \)

\(x65 = {N[NL65/nc, 3], N[NH65/nc, 3]}; 
MatrixForm[x65]\)

\(x88 = {N[NL88/nc, 3], N[NH88/nc, 3]}; 
MatrixForm[x88]\)

\(pLL = N[Sum[LowLow[\([i]\)], {i, 1, nc, 1}]/NL65, 3]; 

    pLH = N[Sum[LowHigh[\([i]\)], {i, 1, nc, 1}]/NL65, 3]; 

    pHL = N[Sum[HighLow[\([i]\)], {i, 1, nc, 1}]/NH65, 3]; 

    pHH = N[Sum[HighHigh[\([i]\)], {i, 1, nc, 1}]/NH65, 3]; \)

TM = {{pLL, pHL}, {pLH, pHH}};

MatrixForm[TM]

g[t_] := MatrixForm[MatrixPower[TM, t] . x65];

g[0]g[1]g[2]g[3]g[4]g[100000]

Gathering all the baseline results and putting them in a text file (Appendix Table)

These next few commands gather together all our key results thus far and generate a "master array" which is then exported as a text file so that it can be read in a spreadsheet program. The export datafiles are organized as follows:  S1 - trap, y_low, y_high, a, a_low, a_high, S2 - x, x_low, x_high, wwedge, logTFPm, logTFPa, S3 - lambda, wbcode, region, bgcode, SocInf, EngFrac, S4 - EurFrac, Latitude, EcOrg, YrsOpen, GADP, logFrankRomer, S5 - Mining, a88, s88adj, s88, logKL, logHL, S6 - logD, logL

SetDirectory[ParentDirectory[]]

CR = FullCalibration[rdlist, nc, 0.35, 0.4, 0.2, 0.1, 0];

\(S1 = 
      Table[If[j == 1, CR[\([i, 1]\)], 
          If[j == 2, CR[\([i, 2]\)], 
            If[j == 3, CR[\([i, 3]\)], 
              If[j == 4, CR[\([i, 4]\)], 
                If[j == 5, CR[\([i, 5]\)], CR[\([i, 6]\)]]]]]], {i, 1, nc, 
          1}, {j, 1, 6, 1}]; 

    S2 = Table[
        If[j == 1, CR[\([i, 7]\)], 
          If[j == 2, CR[\([i, 8]\)], 
            If[j == 3, CR[\([i, 9]\)], 
              If[j == 4, CR[\([i, 10]\)], 
                If[j == 5, CR[\([i, 11]\)], CR[\([i, 12]\)]]]]]], {i, 1, nc, 
          1}, {j, 1, 6, 1}]; 

    S3 = Table[
        If[j == 1, CR[\([i, 13]\)], 
          If[j == 2, CR[\([i, 14]\)], 
            If[j == 3, CR[\([i, 15]\)], 
              If[j == 4, rdlist[\([i, 9]\)], 
                If[j == 5, rdlist[\([i, 11]\)], rdlist[\([i, 12]\)]]]]]], {i, 
          1, nc, 1}, {j, 1, 6, 1}]; 

    S4 = Table[
        If[j == 1, rdlist[\([i, 13]\)], 
          If[j == 2, rdlist[\([i, 14]\)], 
            If[j == 3, rdlist[\([i, 15]\)], 
              If[j == 4, rdlist[\([i, 16]\)], 
                If[j == 5, rdlist[\([i, 17]\)], rdlist[\([i, 18]\)]]]]]], {i, 
          1, nc, 1}, {j, 1, 6, 1}]; 

    S5 = Table[
        If[j == 1, rdlist[\([i, 19]\)], 
          If[j == 2, rdlist[\([i, 1]\)], 
            If[j == 3, rdlist[\([i, 2]\)], 
              If[j == 4, rdlist[\([i, 20]\)], 
                If[j == 5, rdlist[\([i, 4]\)], rdlist[\([i, 5]\)]]]]]], {i, 
          1, nc, 1}, {j, 1, 6, 1}]; 

    S6 = Table[
        If[j == 1, rdlist[\([i, 6]\)], rdlist[\([i, 7]\)]], {i, 1, nc, 1}, {
          j, 1, 2, 1}]; \)

OutputForm[TableForm[N[S1, 6], TableAlignments  ->  {Right},\(TableHeadings\  -> 

              \ \ \ \ \ \ \ \ \ \ \ {
                None, {"trap", "y_low", "y_high", "a", 
                  "a_low", "a_high"}}\)]]>>s1.datOutputForm[TableForm[N[S2, 6], TableAlignments  ->  {Right},\(TableHeadings\  -> 

              \ \ \ \ \ \ \ \ \ \ \ {
                None, {"x", "x_low", "x_high", "wwedge", 
                  "logTFPm", "logTFPa"}}\)]]>>s2.datOutputForm[TableForm[N[S3, 6], TableAlignments  ->  {Right},\(TableHeadings\  -> 

              \ \ \ \ \ \ \ \ \ \ \ {
                None, {"lambda", "wbcode", "region", 
                  "bgcode", "SocInf", "EngFrac"}}\)]]>>s3.datOutputForm[TableForm[N[S4, 6], TableAlignments  ->  {Right},\(TableHeadings\  -> 

              \ \ \ \ \ \ \ \ \ \ \ {
                None, {"EurFrac", "Latitude", "EcOrg", 
                  "YrsOpen", "GADP", "FrankRomer"}}\)]]>>s4.datOutputForm[TableForm[N[S5, 6], TableAlignments  ->  {Right},\(TableHeadings\  -> 

              \ \ \ \ \ \ \ \ \ \ \ {
                None, {"Mining", "a88", "s88adj", "s88", 
                  "logKL", "logHL"}}\)]]>>s5.datOutputForm[TableForm[N[S6, 6], TableAlignments  ->  {Right},\(TableHeadings\  -> 

              \ \ \ \ \ \ \ \ \ \ \ {None, {"logL", "LogD"}}\)]]>>s6.dat

TableForm[S1]