

In [ ]:

    
%matplotlib inline



In [ ]:

    
import matplotlib.pyplot as plt
import numpy as np
import pandas as pd
import sympy as sym

import solowpy

6. Calibration:

In this task we will look at a simple strategy for calibrating a Solow model with Cobb-Douglas production using data from the Penn World Tables (PWT).

Capital depreciation rate, $\delta$:

The PWT provides estimated depreciation rates that vary across both time and countries. As an estimate for the rate of capital depreciation for country $i$ I use the time-average of $\verb|delta_k|_{it}$ as provided by the PWT.

Capital's share for country $i$ in year $t$, $\alpha_{it}$ is computed as $\alpha_{it} = 1 - \verb|labsh|_{it}$, where $\verb|labsh|_{it}$ is the labor share of output/income for country $i$ in year $t$ provided by the PWT. I then use the time-average of $\alpha_{it}$ as the estimate of capital's share for country $i$.

Savings rate, $s$:

As a measure of the savings rate for country $i$, I take the simple time-average of the annual investment share of real GDP, $\verb|i_sh|$, for country $i$ from the PWT.

Labor force growth rate, $n$:

To obtain a measure of the labor force growth rate for country $i$, I regress the logarithm of total employed persons, $\verb|emp|$, in country $i$ from the PWT on a constant and a linear time trend.

$$ \ln\ \verb|emp|_i = \beta_0 + \beta_1 \verb|t| + \epsilon_i \tag{2.1}$$

The estimated coefficient $\beta_1$ is then used as my estimate for the $n$. To estimate the initial number of employed persons, $L_0$, I use $e^{\beta_0}$ as the estimate for $L_0$.

Growth rate of technology, $g$:

To obtain a measure of the growth rate of technology for country $i$, I first adjust the total factor productivity measure reported by the PWT, $\verb|rtfpna|$ (which excludes the human capital contribution to TFP), and then regress the logarithm of this adjusted measure of TFP, $\verb|atfpna|$, for country $i$ on a constant and a linear time trend.

$$ \ln\ \verb|atfpna|_i = \beta_0 + \beta_1 \verb|t| + \epsilon_i \tag{2.2}$$

The estimated coefficient $\beta_1$ is then used as my estimate for the $g$. To estimate the initial level of technology, $A_0$, I use $e^{\beta_0}$ as the estimate for $A_0$.

All of this is being done "behind the scenes". All you need to in order to calibrate the model is pick an ISO 3 country code! Now let's calibrate a Solow model for the UK.



In [ ]:

    
import pypwt



In [ ]:

    
pwt = pypwt.load_pwt_data()



In [ ]:

    
fig, ax = plt.subplots(1, 1, figsize=(8,6))

for ctry in pwt.major_axis:
    tmp_data = pwt.major_xs(ctry)
    tmp_data.labsh.plot(color='gray', alpha=0.5)
    
# plot some specific countries
pwt.major_xs('USA').labsh.plot(color='b', ax=ax, label='USA')
pwt.major_xs('IND').labsh.plot(color='g', ax=ax, label='IND')

# plot global average
avg_labor_share = pwt.labsh.mean(axis=0)
avg_labor_share.plot(color='r', ax=ax)

ax.set_xlabel('Year', family='serif', fontsize=15)
ax.set_ylabel('Labor share of income', family='serif', fontsize=15)
ax.set_ylim(0, 1)

plt.show()



In [ ]:

    
def match_moments(model, data, iso3_code, bounds=None):
    r"""
    Simple calibration scheme for a Solow model with Cobb-Douglas production
    based on data from the Penn World Tables (PWT).

    Parameters
    ----------

    model : solow.CobbDouglasModel
        An instance of the CobbDouglasModel class that you wish to calibrate.
    data : pandas.Panel
        An instance of the pandas.Panel class containing PWT data.
    iso3_code : str
        A valid ISO3 country code. For example, to calibrate the model using
        data for the United States, one would set iso3_code='USA'; to calibrate
        a model using data for Zimbabwe, one would set iso3_code='ZWE'. For a
        complete listing of ISO3 country codes see `wikipedia`_.
    bounds : tuple (default=None)
        Start and end years for the subset of the PWT data to use for
        calibration. Note that start and end years should be specified as
        strings. For example, to calibrate a model using data from 1983 to 2003
        one would set

            bounds=('1983', '2003')

        By default calibration will make use of all available data for the
        specified country.

    .. `wikipedia`: http://en.wikipedia.org/wiki/ISO_3166-1_alpha-3

    """
    # get the PWT data for the iso_code
    tmp_data = data.major_xs(iso3_code)
    required_vars = ['rgdpna', 'rkna', 'emp', 'labsh', 'csh_i', 'delta_k']
    tmp_data = tmp_data[required_vars].dropna()
    assert not tmp_data.empty, "Insufficient data to estimate model parameters"
    
    # set bounds
    if bounds is None:
        start = tmp_data.index[0]
        end = tmp_data.index[-1]
    else:
        start = bounds[0]
        end = bounds[1]
    tmp_data = tmp_data.loc[start:end]

    # define the data used in the calibration
    output = tmp_data['rgdpna']
    capital = tmp_data['rkna']
    labor = tmp_data['emp']
    labor_share = tmp_data['labsh']
    savings_rate = tmp_data['csh_i']
    depreciation_rate = tmp_data['delta_k']

    # define a time trend variable
    N = tmp_data.index.size
    linear_trend = pd.Series(np.linspace(0, N - 1, N), index=tmp_data.index)
    time_trend = linear_trend.loc[start:end]

    # estimate capital's share of income/output
    capital_share = 1 - labor_share
    alpha = capital_share.mean()

    # compute solow residual (note dependence on alpha!)
    solow_residual = model.evaluate_solow_residual(output, capital, labor)
    technology = solow_residual.loc[start:end]

    # estimate the fraction of output saved
    s = savings_rate.mean()

    # regress log employed persons on linear time trend
    res = pd.ols(y=np.log(labor), x=time_trend)
    n = res.beta[0]
    L0 = np.exp(res.beta[1])

    # regress log TFP on linear time trend
    res = pd.ols(y=np.log(technology), x=time_trend)
    g = res.beta[0]
    A0 = np.exp(res.beta[1])

    # estimate the depreciation rate for total capital
    delta = depreciation_rate.mean()

    # create a dictionary of model parameters
    tmp_params = {'s': s, 'alpha': alpha, 'delta': delta, 'n': n, 'L0': L0,
                  'g': g, 'A0': A0}

    # update the model's parameters
    return tmp_params



In [ ]:

    
# define model parameters
cobb_douglas_params = {'A0': 1.0, 'L0': 1.0, 'g': 0.02, 'n': 0.03, 's': 0.15,
                      'delta': 0.05, 'alpha': 0.33}

# create an instance of the solow.Model class
cobb_douglas_model = solow.CobbDouglasModel(params=cobb_douglas_params)



In [ ]:

    
match_moments(cobb_douglas_model, pwt, 'CHN')

Estimating model with CES production

Maximum likelihood



In [ ]:

    
# define model parameters
ces_params = {'K0': 1.0, 'A0': 1.0, 'L0': 1.0, 'g': 0.02, 'n': 0.03, 's': 0.15,
              'delta': 0.05, 'alpha': 0.33, 'sigma': 0.95, 'sigma_eps': 0.01}

# create an instance of the solow.Model class
ces_model = solow.CESModel(params=ces_params)



In [ ]:

    
from scipy import optimize, stats



In [ ]:

    
def _intialize_labor_supply(model, labor):
    """Initial condition for labor supply is not a free parameter."""
    
    
def _cobb_douglas_initial_guess(model, data, iso3_code, bounds, sigma0, sigma_eps0):
    """Estimate parameters assumming Cobb Douglas production."""
    tmp_params = match_moments(model, data, iso3_code, bounds)
    
    # initial guesses for sigma, and sigma_eps
    tmp_params['sigma'] = sigma0
    tmp_params['sigma_eps'] = sigma_eps0
    
    return _params_to_array(tmp_params) 

def _array_to_params(array):
    """Converts array to dictionary of model params."""
    keys = ['g', 'n', 's', 'alpha', 'delta', 'sigma', 'sigma_eps']
    params = dict(zip(keys, array))
    return params

def _params_to_array(params):
    """Converts dictionary of model params to an array."""
    g = params['g']
    n = params['n']
    s = params['s']
    alpha = params['alpha']
    delta = params['delta']
    sigma = params['sigma']
    sigma_eps = params['sigma_eps']
    return np.array([g, n, s, alpha, delta, sigma, sigma_eps])
    
def _likelihood_function(params_array, model, output, capital, labor, labor_share):
    """Objective for maximum likelihood estimation."""
    # create tmp parameter dictionary
    tmp_params = model.params.copy()
    new_params = _array_to_params(params_array)
    tmp_params.update(new_params)
    model.params = tmp_params

    # compute solow residual (note dependence on alpha and sigma!)
    technology = model.evaluate_solow_residual(output, capital, labor)
    
    # compute capital (per unit effective labor)
    k = capital / (technology * labor)

    # model predicted labor's share of income/output
    predicted_labor_share = 1 - model.evaluate_output_elasticity(k)

    # residual difference between model predictions and data
    residual = labor_share - predicted_labor_share

    # assumes residuals are gaussian
    total_log_likelihood = np.sum(np.log(stats.norm.pdf(residual, 0, model.params['sigma_eps'])))

    return -total_log_likelihood



In [ ]:

    
def maximize_likelihood(model, data, iso3_code, bounds=None, method='BFGS', sigma0=1.01, sigma_eps0=0.01):
    r"""
    Maximum likelihood estimation scheme for a Solow model with CES
    production based on data from the Penn World Tables (PWT).

    Parameters
    ----------

    model : solow.CESModel
        An instance of the CESModel class that you wish to estimate.
    data : pandas.Panel
        An instance of the pandas.Panel class containing PWT data.
    iso3_code : str
        A valid ISO3 country code. For example, to calibrate the model using
        data for the United States, one would set iso3_code='USA'; to calibrate
        a model using data for Zimbabwe, one would set iso3_code='ZWE'. For a
        complete listing of ISO3 country codes see `wikipedia`_.
    bounds : tuple (default=None)
        Start and end years for the subset of the PWT data to use for
        calibration. Note that start and end years should be specified as
        strings. For example, to calibrate a model using data from 1983 to 2003
        one would set

            bounds=('1983', '2003')

        By default calibration will make use of all available data for the
        specified country.

    .. `wikipedia`: http://en.wikipedia.org/wiki/ISO_3166-1_alpha-3

    """
    # get the PWT data for the iso_code
    tmp_data = data.major_xs(iso3_code)
    required_vars = ['rgdpna', 'rkna', 'emp', 'labsh']
    tmp_data = tmp_data[required_vars].dropna()
    assert not tmp_data.empty, "Insufficient data to estimate model parameters"
    
    # set bounds
    if bounds is None:
        start = tmp_data.index[0]
        end = tmp_data.index[-1]
    else:
        start = bounds[0]
        end = bounds[1]
    tmp_data = tmp_data.loc[start:end]

    # define the data used in the calibration
    output = tmp_data['rgdpna'].loc[start:end]
    capital = tmp_data['rkna'].loc[start:end]
    labor = tmp_data['emp'].loc[start:end]
    labor_share = tmp_data['labsh'].loc[start:end]
    technology = model.evaluate_solow_residual(output, capital, labor)

    initial_guess = _cobb_douglas_initial_guess(model, data, iso3_code, bounds,
                                                sigma0, sigma_eps0)
    
    eps = 1e-6
    bnds = [(None, None), (None, None), (eps, 1-eps), (eps, 1-eps),
            (eps, None), (eps, None), (eps, None)]
    result = optimize.minimize(_likelihood_function,
                               x0=initial_guess,
                               args=(model, output, capital, labor, labor_share),
                               method=method,
                               bounds=bnds)
    if result.success:
        tmp_params = _array_to_params(result.x)
        tmp_params['A0'] = technology[0]
        tmp_params['L0'] = labor[0]
        tmp_params['K0'] = capital[0]
    
        return result, tmp_params
    
    else:
        return result



In [ ]:



In [ ]:

    
# Add assertion error to evaluate solow_residual method@
result, params = maximize_likelihood(ces_model, pwt, 'USA', bounds=(None, None),
                                     method='Nelder-Mead', sigma0=0.5, sigma_eps0=0.01)



In [ ]:

    
result



In [ ]:

    
ces_model.params = params



In [ ]:

    
result



In [ ]:

    
def _initial_condition(model):
    k0 = model.params['K0'] / (model.params['A0'] * model.params['L0'])
    return k0

# need to specify some initial conditions
t0, k0 = 0.0, _initial_condition(ces_model)
numeric_soln = ces_model.ivp.solve(t0, k0, T=62, integrator='dopri5')



In [ ]:

    
fig, ax = plt.subplots(1, 1, figsize=(8,6))

ax.plot(numeric_soln[:,0], 1 - ces_model.evaluate_output_elasticity(numeric_soln[:,1]),
        'bo', markersize=3.0)

ax.plot(pwt.major_xs('USA').labsh.loc[:].values)

# axes, labels, title, etc
ax.set_xlabel('Time, $t$', fontsize=15, family='serif')
ax.set_ylabel(r'$\alpha_K(k(t))$', rotation='horizontal', fontsize=20, family='serif')
ax.set_title('Labor share of income', fontsize=20, family='serif')
ax.legend(loc=0, frameon=False, bbox_to_anchor=(1.0, 1.0))
ax.grid('on')

plt.show()



In [ ]:

    
fig, ax = plt.subplots(1, 1, figsize=(8,6))

intensive_output = ces_model.evaluate_intensive_output(numeric_soln[:,1])
technology = ces_model.params['A0'] * np.exp(ces_model.params['g'] * numeric_soln[:,0])
output_per_capita = intensive_output * technology

ax.plot(numeric_soln[1:,0], np.diff(np.log(output_per_capita)),
        'bo', markersize=3.0)

Y = pwt.major_xs('USA').rgdpna.loc[:]
E = pwt.major_xs('USA').emp.loc[:]
ax.plot(np.log(( Y / E )).diff().values)

# axes, labels, title, etc
ax.set_xlabel('Time, $t$', fontsize=15, family='serif')
ax.set_ylabel(r'$g$', rotation='horizontal', fontsize=20, family='serif')
ax.set_title('Growth rate of GDP (per unit labor)', fontsize=20, family='serif')
ax.legend(loc=0, frameon=False, bbox_to_anchor=(1.0, 1.0))
ax.grid('on')

plt.show()



In [ ]:

    
np.diff?



In [ ]:

    
(d.rgdpna / d.emp).pct_change()



In [ ]: