Housekeeping...

In [2]:
# Classic libraries
from __future__ import division
import numpy as np
import sympy as sym
# Plotting libraries
import matplotlib.pyplot as plt
import seaborn as sn # library for pretty plots. See https://stanford.edu/~mwaskom/software/seaborn/ 
sn.set_style('whitegrid')
cp = sn.xkcd_rgb["pinkish"]
cb = "#3498db"
cr = "#1fa774"
# Specialist functions
from scipy.optimize import minimize, root, fsolve
from scipy.optimize import leastsq, least_squares
from scipy.interpolate import interp1d
from scipy.special import erf
# Core functions - solve the problem and plot the results
import capital
# Other extensions
%matplotlib inline
%load_ext autoreload
%autoreload 2

Loading data from Adamopoulos and Restuccia (2014)

In [3]:
CalUS = np.loadtxt('model_data.csv', delimiter=',').T
CalP = np.loadtxt('model_data_p.csv', delimiter=',').T
y_grid, y_prob, farm_grid, k_grid, output_m = CalUS
y_grid_p, y_prob_p, farm_grid_p, k_grid_p, output_m_p = CalP

Replicating Adamopoulos and Restuccia (2014) and some more using Eeckhout and Kircher (2016)

(or how heterogeneity in the quality of land cahnges the results or A&R)

Index

  1. Introduction
  2. Calculating $k^*$
  3. Solving original model
  4. Uneven spread of $x$
    1. Generating the distributions
    2. Solving the model
    3. Results
  5. The PAM condition

1. Introduction

The original A&R paper's production function is:

$$ p_aA\kappa(\eta(k)^{\rho} + (1 -\eta)(yl))^{\rho})^{\frac{\gamma}{\rho}} $$

Where:

  • $y \to$ farmer skill
  • $l \to$ farm size
  • $k \to$ capital

The rest are parameters that are calibrated to match the farm size distribution of the US and the capital to land ratio.

E&K introduce two more variables:

  • $x \to$ land quality
  • $r \to$ farmer supervising intentisty, normalized to 1.

Adding heterogenity in land quality ($x$) to the A&R model requires a different solution method.

The optimal choice of $k$ can be calculated as a function of the other variables so that the problem can be solved as an initial value problem (IVP) in the standard E&K way. That is, finding $k^*$ such that $$k^*(x,y,l,r) = \mathrm{argmax}_{k}F(x,y,l,r,k)-rRk$$

There is a number of ways of introducing these $x$ and $r$, of which 5 were considered. Below is the list of all tehm, with their accompaning first order condiction (FOC) for solving for $k^*(x,y,l,r)$.

# name function FOC
0 $\color{red}{x}$ with $l/r$ $ rp_aA\kappa(\eta(k)^{1/4} + (1 -\eta)(y(l/r)^{\color{red}{x}})^{1/4})^{2} - rRk$ $ K^{3} + -\frac{Ap_a\kappa\eta^{2}}{2R} K + \frac{Ap_a\kappa\eta(1 -\eta)(y(\frac{l}{r})^{x})^{1/4}}{2R} = 0$
1 $\color{red}{x}$ with $A$ $ rp_a\color{red}{x}A\kappa [\eta (k)^{1/4} + (1-\eta)(y(l/r))^{1/4} ]^2 - rRk $ $K^3 - \frac{xAp_a\kappa\eta^2}{2R}K - \frac{xAp_a\kappa\eta(1-\eta)(y\frac{l}{r})^{1/4}}{2R}= 0 $
2 $\color{red}{x}$ with $k$ $ rp_aA\kappa [\eta (k\color{red}{x})^{1/4} + (1-\eta)(y(l/r))^{1/4} ]^2 - rRk \color{red}{x} $ $K^3 - \frac{Ap_a\kappa\eta^2 x^{-1/2}}{2R}K - \frac{Ap_a\kappa\eta(1-\eta)(y\frac{l}{r})^{1/4} x^{-3/4}}{2R}= 0 $
3 $\color{red}{x}$ with $k$ (out of the bracket) $ rp_aA\kappa [\eta \color{red}{x} k^{1/4} + (1-\eta)(y(l/r))^{1/4} ]^2 - rRk \color{red}{x} $ $K^3 - \frac{Ap_a\kappa\eta^2 x^2}{2R}K - \frac{Ap_a\kappa\eta(1-\eta)(y\frac{l}{r})^{1/4} x}{2R}= 0 $
4 $k^{\color{red}{x}}$ $ rp_aA\kappa(\eta(k)^{\color{red}{x}/4} + (1 -\eta)(y(l/r))^{1/4})^{2} - rRk$ None

Note: $K = k^{1/4}$

For the case of $k^{\color{red}{x}}$, there is no closed-form solution and $k^*$ has to be solved numerically. Solving the model with this especification renders the shooting solver I use unstable, therefore I do not use it. However, capital.py includes a function to carry this approximation out.

For the rest of the specifications, $k^*$ is found by solving a cubic root. This is done in section 2. capital.py approximates the solution with cubic or quadratic "tensor" polynomials - that is, polynomials with cross-products. This solves a multiple root problem and gives an easily differentiable expression (tht would make the solver work faster).

Once $k^*$ has been found, I proceed to solve for the baseline case where $x=1$ (section 3) and then increase the spread up to $x\in[0.75,1.25]$ (section 4). I use a shooting algorithm with logaritmic steps.

I chose to give $x$ a lognormal distribution (with mean 1) because is the one that makes the results move closest to the data - that is, small farms increasing size, but not much. I have left the code here to try other distributions and spreads.

Finally, in section 5 I provide the symolic derivatives to the production function to show that with A&R parameters, PAM holds for any combination of $x,y,l,r$ and $k$.

2 Calculating $k^*$

In the past, I calculated the optimal choice of $k$ for different functional forms - that's why I creaded dictionaries with the keys being the different functional forms the production function could take. I have keep it that way to bear in mind that I chose function (2), but you can solve for $k$ for any of the four different production functions.

I also have different $k^*$ approximations depending on the maximum spread of $x$ (here [0.90,1.10] and [0.75,1.25]). Shorter spreads lead to more accurate results.

Rich country (US calibration):

In [4]:
k_stars25 = {}
k_stars10 = {}

k_stars25[2] = capital.get_k("R", ftype=2,fchoice=1,plot=False, save_in="Rkpams25_{}".format(2), x_range=[0.75,1.25])
k_stars10[2] = capital.get_k("R", ftype=2,fchoice=0,plot=False, save_in="Rkpams{}".format(2), x_range=[0.9,1.1])

Poor country:

In [5]:
k_starsP25 = {}
k_starsP10 = {}

k_starsP25[2] = capital.get_k("P", ftype=2,fchoice=1,plot=False, save_in="Pkpams25_{}".format(2), x_range=[0.75,1.25])
k_starsP10[2] = capital.get_k("P", ftype=2,fchoice=0,plot=False, save_in="Pkpams{}".format(2), x_range=[0.9,1.1])

These $k$s will be feed to the solver.

3 Solving original model ($x=1$)

capital.py solves for no spread of $x$ first - the default values of the spread of $x$ being $(0.99999,1.00001)$.

In [6]:
sol_R1 = capital.solve("R", k_stars10[2], ftype=2,assort='positive',scaling_x=0.90607341901223637)
sol_P1 = capital.solve("P", k_starsP10[2], ftype=2,assort='positive',scaling_x=0.90607341901223637)

Note: In order for the solver to match the lowest values of $y$ (0.000003) capital.solve uses an equivalent production function

$$lp_aA\kappa(\eta(yk)^{1/4} + (1 -\eta)(x(r/l))^{1/4})^{2} - lRk$$

This does not alter the results, but requires some care when dealing with the results: farm size ($\theta = \frac{l}{r}$) coming from the solver would be the inverse of model $\theta$, and in the code below you'll find xs refering to farmer skill ($y$ in the model).

Histogram comparing with the original paper (Figure 7 in the appendix):

In [7]:
capital.plot_histogram(sol_R1,sol_P1,CalUS,CalP,save_in='Figure7.pdf')

For comparison later on, I solve for increasing the spread evenly (uniform distribution) for different ranges of $x$:

In [7]:
sol_R10 = {}
sol_P10 = {}
x_ranges = [(0.975,1.025),(0.95,1.05),(0.925,1.075),(0.9,1.1)]
for i in range(len(x_ranges)):
    sol_R10[i] = capital.solve("R", k_stars10[2], ftype=2,assort='positive',x_range=x_ranges[i],scaling_x=0.90630885374008696)
    sol_P10[i] = capital.solve("P", k_starsP10[2], ftype=2,assort='positive',x_range=x_ranges[i],scaling_x=0.90613690136320324)
    print "Done with range {}".format(x_ranges[i])
Done with range (0.975, 1.025)
Done with range (0.95, 1.05)
Done with range (0.925, 1.075)
Done with range (0.9, 1.1)

4. Uneven spread of x

Here I introduce different ways of spreading $x$. Although in the end only the lognormal is used, I've left the code written for other distributions as well to illustrate grphically how $x$ is unevenly changed.

4.1 Generating the distributions

Lognormal distribution:

In [8]:
_modules = [{'ImmutableMatrix': np.array, 'erf': erf, 'sqrt': np.sqrt}, 'numpy']

x, loc2, mu2, sigma2 = sym.var('x, loc2, mu2, sigma2')
productivity_cdf = (0.5 + 0.5 * sym.erf((sym.log(x - loc2) - mu2) / sym.sqrt(2 * sigma2**2)))

productivity_params62 = {'loc2': 0.0, 'mu2': 0.00, 'sigma2': 0.20}

pdf62 = productivity_cdf.diff(x).subs(productivity_params62)
pdf_exe62 = sym.lambdify(x, pdf62, 'numpy')
cdf_exe62 = sym.lambdify(x, productivity_cdf.subs(productivity_params62), _modules)

Bimodal

In [9]:
sig1, sig2, x, mu1, mu2= sym.var("sig1, sig2, x, mu1, mu2")
c1 = 0.5 * (1 + sym.erf((x-mu1)/ sym.sqrt(2 * sig1**2)))
c2 = 0.5 * (1 + sym.erf((x-mu2)/ sym.sqrt(2 * sig2**2)))
cdf = c1*0.5+c2*0.5

params72 = {'mu1': 0.93333, 'mu2': 1.06667, 'sig1': 0.033,'sig2': 0.033}

pdf72 = cdf.diff(x).subs(params72)
cdf72 = cdf.subs(params72)

pdf_exe72 = sym.lambdify(x, pdf72,  _modules)
cdf_exe72 = sym.lambdify(x, cdf72,  _modules)

Normal

In [10]:
params92 = {'mu1': 1.0, 'sig1': 0.06666 }

pdf92 = c1.diff(x).subs(params92)
cdf92 = c1.subs(params92)

pdf_exe92 = sym.lambdify(x, pdf92,  _modules)
cdf_exe92 = sym.lambdify(x, cdf92,  _modules)

Normalising constants - so total ammount of land integrates to 1.

In [11]:
norm72 = 1/(cdf_exe72(1.1)-cdf_exe72(0.9))
norm92 = 1/(cdf_exe92(1.1)-cdf_exe92(0.9))
norm62 = 1/(cdf_exe62(1.1)-cdf_exe62(0.9))

c72 = cdf_exe72(1.1)*norm72-1
c92= cdf_exe92(1.1)*norm92-1
c62= cdf_exe92(1.1)*norm62-1

Plots of the different distributions

In [184]:
xss = np.linspace(0.9,1.1,6000)
plt.figure(figsize=(6,6))
sn.set_palette('colorblind', 5)
plt.axvline(1.0,label="x=1", c= 'black')
plt.plot(xss,5*np.ones(len(xss)), label="Uniform")
plt.plot(xss,pdf_exe72(xss)*norm72, label='Bimodal')
plt.plot(xss,pdf_exe62(xss)*norm62, label="Lognormal")
plt.xlim(0.9,1.1)
plt.title("Pdf of land quality", fontsize=16)
plt.plot(xss,pdf_exe92(xss)*norm92, label='Normal')


plt.tight_layout()
plt.legend(loc='best')
plt.show()

4.2 Solving the model

The distribution that used in the paper is the Lognormal (code 26xx). This distribution delivers the results closer to the data provided in Adamopoulos and Restuccia (2014).

In [14]:
# Execute this cell to store results
sol_log = {}
sol_logP = {}

Lognormal

In [15]:
sol_log[2605] = capital.solve("R", k_stars10[2], ftype=2,assort='positive',spread='lognorm',dispams=[0.00,0.2],
                             x_range=[0.95,1.05],scaling_x=norm62*1.7611607186330467, normconst=c62)
sol_logP[2605]= capital.solve("P", k_starsP10[2], ftype=2,assort='positive',spread='lognorm',dispams=[0.00,0.2],
                             x_range=[0.95,1.05],scaling_x=norm62*1.7611607186330467, normconst=c62)
In [16]:
sol_log[2610] = capital.solve("R", k_stars10[2], ftype=2,assort='positive',spread='lognorm',dispams=[0.00,0.2],
                             x_range=[0.9,1.1],scaling_x=norm62*0.9061110680375417, normconst=c62)
sol_logP[2610]= capital.solve("P", k_starsP10[2], ftype=2,assort='positive',spread='lognorm',dispams=[0.00,0.2],
                             x_range=[0.9,1.1],scaling_x=norm62*0.9061110680375417, normconst=c62)
In [17]:
sol_log[2615] = capital.solve("R", k_stars25[2], ftype=2,assort='positive',spread='lognorm',dispams=[0.00,0.2],
                             x_range=[0.85,1.15],scaling_x=1.6494247873450758, normconst=c62)
sol_logP[2615]= capital.solve("P", k_starsP25[2], ftype=2,assort='positive',spread='lognorm',dispams=[0.00,0.2],
                             x_range=[0.85,1.15],scaling_x=1.6494247873450758, normconst=c62)
In [18]:
sol_log[2620] = capital.solve("R", k_stars25[2], ftype=2,assort='positive',spread='lognorm',dispams=[0.00,0.2],
                             x_range=[0.8,1.20],scaling_x=1.3196670116025528, normconst=c62)
sol_logP[2620]= capital.solve("P", k_starsP25[2], ftype=2,assort='positive',spread='lognorm',dispams=[0.00,0.2],
                             x_range=[0.8,1.20],scaling_x=1.3196670116025528, normconst=c62)
In [19]:
sol_log[2625] = capital.solve("R", k_stars25[2], ftype=2,assort='positive',spread='lognorm',dispams=[0.00,0.2],
                             x_range=[0.75,1.25],scaling_x=1.1434977494042775, normconst=c62)
sol_logP[2625]= capital.solve("P", k_starsP25[2], ftype=2,assort='positive',spread='lognorm',dispams=[0.00,0.2],
                             x_range=[0.75,1.25],scaling_x=1.1434977494042775, normconst=c62)

4.3 Results

Changes in farm size for widening spreads of $x$

In [ ]:
solsR = (sol_R1,sol_log[2605],sol_log[2610],sol_log[2615],sol_log[2620],sol_log[2625])
solsP = (sol_P1,sol_logP[2605],sol_logP[2610],sol_logP[2615],sol_logP[2620],sol_logP[2625])
labs = ('$x=0$','$x\in[0.95,1.05]$ ','$x\in[0.90,1.10]$','$x\in[0.85,1.15]$','$x\in[0.80,1.20]$','$x\in[0.75,1.25]$')

In logs:

In [20]:
capital.plot_land_choice(solsR,labs,solsP,labs, color=True, logs=True)

In levels:

In [21]:
capital.plot_land_choice(solsR,labs,solsP,labs, color=True, logs=False)

Focusing on the poor country (Figure 1 in the paper):

In [23]:
capital.plot_land_choice(solsP,labs,None,None, color=True, logs=False,save_in='Figure1.pdf')

$k/l$ ratios

These show the $k/l$ ratio normalized by the largest farm. That is, is $k/l=100$ for the largest farm, the $k/l$ ratio of the rest is divided by 100.

Rinch country on the left, poor country on the right.

In [350]:
capital.plot_kl_ratio(solsR,labs,solsP,labs,save_in=None,color=True)

Distribution plots

Rich country:

In [301]:
capital.plot_dischange(sol_R1,sol_log[2620],extras=[sol_log[2610],(0.9,1.1)],xrange=(0.8,1.2),
                            save_in=None,allplots=False,color=True, uc=cr)

Poor country:

In [281]:
capital.plot_dischange(sol_P1,sol_logP[2620],extras=[sol_logP[2610],(0.9,1.1)],xrange=(0.8,1.2),
                            save_in=None, allplots=False,color=True,uc=cp)

5. PAM condition

Does the PAM condition hold? Here I check for the PAM condition, so if lhs-rhs > 0, then PAM holds.

First without substituiting parameters:

In [478]:
x, y, k = sym.var('x, y, k')
R, l, r, A, kappa, p_a, rho, gamma, eta = sym.var('R, l, r, A, kappa, p_a, rho, gamma, eta')

sym.init_printing()
F = r*p_a*A*kappa*(eta*(x*k)**rho + (1- eta)*((l/r)*y)**rho)**(gamma/rho)
lhs = F.diff(x,y)*F.diff(l,r)*F.diff(k,k) - F.diff(x,y)*F.diff(l,k)*F.diff(r,k) - F.diff(x,k)*F.diff(y,k)*F.diff(l,r)
rhs = F.diff(x,r)*F.diff(y,l)*F.diff(k,k) - F.diff(x,r)*F.diff(y,k)*F.diff(l,k) - F.diff(x,k)*F.diff(y,l)*F.diff(r,k)
sym.simplify(lhs-rhs)
Out[478]:
$$- \frac{A^{3} \eta^{2} \gamma^{3} \kappa^{3} p_{a}^{3} r^{2} \rho \left(k x\right)^{2 \rho} \left(\frac{l y}{r}\right)^{\rho} \left(\eta - 1\right) \left(\eta \left(k x\right)^{\rho} - \eta \left(\frac{l y}{r}\right)^{\rho} + \left(\frac{l y}{r}\right)^{\rho}\right)^{\frac{3 \gamma}{\rho}}}{k^{2} l x y \left(\eta^{3} \left(k x\right)^{3 \rho} - 3 \eta^{3} \left(k x\right)^{2 \rho} \left(\frac{l y}{r}\right)^{\rho} + 3 \eta^{3} \left(k x\right)^{\rho} \left(\frac{l y}{r}\right)^{2 \rho} - \eta^{3} \left(\frac{l y}{r}\right)^{3 \rho} + 3 \eta^{2} \left(k x\right)^{2 \rho} \left(\frac{l y}{r}\right)^{\rho} - 6 \eta^{2} \left(k x\right)^{\rho} \left(\frac{l y}{r}\right)^{2 \rho} + 3 \eta^{2} \left(\frac{l y}{r}\right)^{3 \rho} + 3 \eta \left(k x\right)^{\rho} \left(\frac{l y}{r}\right)^{2 \rho} - 3 \eta \left(\frac{l y}{r}\right)^{3 \rho} + \left(\frac{l y}{r}\right)^{3 \rho}\right)}$$

Then substituiting all the parameters fixed across countries:

In [479]:
eta = 0.89
rho = 0.25
gamma = 0.5
kappa = 1

F = r*p_a*A*kappa*(eta*(x*k)**rho + (1- eta)*((l/r)*y)**rho)**(gamma/rho)
lhs = F.diff(x,y)*F.diff(l,r)*F.diff(k,k) - F.diff(x,y)*F.diff(l,k)*F.diff(r,k) - F.diff(x,k)*F.diff(y,k)*F.diff(l,r)
rhs = F.diff(x,r)*F.diff(y,l)*F.diff(k,k) - F.diff(x,r)*F.diff(y,k)*F.diff(l,k) - F.diff(x,k)*F.diff(y,l)*F.diff(r,k)
sym.simplify(lhs-rhs)
Out[479]:
$$\frac{0.00272284375 A^{3} p_{a}^{3} r^{2} \left(k x\right)^{0.5} \left(\frac{l y}{r}\right)^{0.25}}{k^{2} l x y} \left(0.89 \left(k x\right)^{0.25} + 0.11 \left(\frac{l y}{r}\right)^{0.25}\right)^{3.0}$$