# Numerical Methods for Economists: Lab Assignment The Solow (1956) Model¶

In this lab assignment you will analyze a version of the Solow model with a constant elasticity of substituion (CES) production function and attempt to "explain" observed patterns in capital's share in the U.S. between 1950-2011.

A CES production function looks as follows...

$$Y(t) = \bigg[\alpha K(t)^{\rho} + (1-\alpha) (A(t)L(t))^{\rho}\bigg]^{\frac{1}{\rho}} \tag{1.2}$$

where $0 < \alpha < 1$ and $-\infty < \rho < 1$. The parameter $\rho = \frac{\sigma - 1}{\sigma}$ where $\sigma$ is the elasticity of substitution between factors of production. The CES production technology is popular because it nests several interesing special cases. In particular, if factors of production are perfect substitutes (i.e., $\sigma = +\infty \implies \rho = 1$), then output is just a linear combination of the inputs.

$$\lim_{\rho \rightarrow 1} Y(t) = \alpha K(t) + (1-\alpha)A(t)L(t) \tag{1.3}$$

On the other hand, if factors of production are perfect complements (i.e., $\sigma = 0 \implies \rho = -\infty$), then we recover the Leontief production function.

$$\lim_{\rho \rightarrow -\infty} Y(t) = \min\left\{\alpha K(t), (1-\alpha) A(t)L(t)\right\} \tag{1.4}$$

Finally, if the elasticity of substitution is unitary (i.e., $\sigma=1 \implies \rho=0$), then output is Cobb-Douglas.

$$\lim_{\rho \rightarrow 0} Y(t) = K(t)^{\alpha}(A(t)L(t))^{1-\alpha} \tag{1.5}$$
In [ ]:
import numpy as np
import pandas as pd
from scipy import integrate, linalg, optimize
import matplotlib as mpl
import matplotlib.pyplot as plt

# for the first few labs we will be working with models of growth
import growth
import pwt


### (5 points) Part a)¶

Show that the CES production function as defined by equation 1.2 exhibits constant returns to scale.

In [ ]:



### (10 points) Part b)¶

Derive the intensive forms for both the general CES production function as defined by equation 1.2 and for its Cobb-Douglas special case defined by equation 1.5. Show that both of these intensive production functions are concave.

In [ ]:



### (5 points) Part c)¶

Using your results from part b, complete the Python function below which defines output (per person/effective person) in terms of capital (per person/effective person) and model parameters.

In [ ]:
def ces_output(t, k, params):
"""
Intensive form for a CES production
function.

Arguments:

t:      (array-like) Time.
k:      (array-like) Capital (per person/effective person).
params: (dict) Dictionary of parameter values.

Returns:

y: (array-like) Output (per person/ effective person)

"""
# extract params
alpha = params['alpha']
sigma = params['sigma']
rho   = # INSERT CODE HERE DEFINING RHO IN TERMS OF SIGMA!

# nest Cobb-Douglas technology as special case
if rho == 0:
y = # INSERT CODE HERE!
else:
y = # INSERT CODE HERE!

return y


### (5 points) Part d)¶

Using your results from part b, complete the Python function below which defines the marginal product of capital (per person/effective person) in terms of capital (per person/effective person) and model parameters.

In [ ]:
def marginal_product_capital(t, k, params):
"""
Marginal product of capital with CES production function.

Arguments:

t:      (array-like) Time.
k:      (array-like) Capital (per person/effective person).
params: (dict) Dictionary of parameter values.

Returns:

mpk: (array-like) Derivative of output with respect to capital, k.

"""
# extract params
alpha = params['alpha']
sigma = params['sigma']
rho   = # INSERT CODE HERE DEFINING RHO IN TERMS OF SIGMA!

# nest Cobb-Douglas technology as special case
if rho == 0:
mpk = # INSERT CODE HERE!
else:
mpk = # INSERT CODE HERE!

return mpk


### (5 points) Part e)¶

Derive the equation of motion for capital per effective worker when the general CES production function defined by equation 1.2 and for the Cobb-Douglas special case defined by equation 1.5.

In [ ]:



### (5 points) Part f)¶

Using your results from parts d) and e) and the Python function defining the equation of motion for capital (per person/effective person) complete the Python function which defines the Jacobian for the Solow model.

In [ ]:
def equation_of_motion_capital(t, k, params):
"""
Equation of motion for capital (per person/effective person).

Arguments:

t:      (array-like) Time.
k:      (array-like) Capital (per person/effective person).
params: (dict) Dictionary of parameter values.

Returns:

k_dot: (array-like) Time-derivative of capital (per person/effective
person).

"""
# extract params
s     = params['s']
n     = params['n']
g     = params['g']
delta = params['delta']

y = ces_output(t, k, params)

k_dot = s * y - (n + g + delta) * k

return k_dot

In [ ]:
def solow_jacobian(t, k, params):
"""
The Jacobian of the Solow model.

Arguments:

t:      (array-like) Time.
k:      (array-like) Capital (per person/effective person).
params: (dict) Dictionary of parameter values.

Returns:

jac: (array-like) Value of the derivative of the equation of
motion for capital (per worker/effective worker) with
respect to k.

"""
# extract params
s     = params['s']
n     = params['n']
g     = params['g']
delta = params['delta']

mpk = #INSERT YOUR CODE HERE

k_dot = s * mpk - (n + g + delta)

return k_dot


### (5 points) Part g)¶

In the cell below, create a Python dictionary called default_params using the following values: $\alpha$ = 0.33, $\delta$ = 0.04, $\sigma$ = 1.0, $g$ = 0.02, $n$ = 0.01, $s$ = 0.15, $L(0)$=1.0, $A(0)$=1.0.

In [ ]:
# INSERT CODE HERE!


### (5 points) Part h)¶

In the cell below, create instance of the SolowModel class called model using the results from parts c), d), f) and g).

In [ ]:
# INSERT CODE HERE!


### (10 points) Part i)¶

Derive an analytic expression for the steady state value of capital (per person/effective person) using your result from part e) and use your result to complete the Python function in the cell below defining the analytic steady state value for $k$ as a function of model parameters.

In [ ]:
def analytic_k_star(params):
"""Steady-state level of capital (per person/effective person)."""
# extract params
s     = params['s']
n     = params['n']
g     = params['g']
alpha = params['alpha']
delta = params['delta']
sigma = params['sigma']
rho   = # INSERT CODE HERE DEFINING RHO IN TERMS OF SIGMA!

# nest Cobb-Douglas technology as special case
if rho == 0:
k_star = # INSERT CODE HERE!
else:
k_star = # INSERT CODE HERE!

return k_star


### (5 points) Part j)¶

Using the code from the lab session as a guide, create a Python dictionary called solow_steady_state_funcs containing your result from part i) and add it to the model.

In [ ]:
# INSERT YOUR CODE HERE!!


Execute the code in the cell below to calibrate the model for the U.S. using data from the Penn World Tables.

In [ ]:
# calibrate the model!
growth.calibrate_ces(model, 'USA', x0=[0.5, 0.5], bounds=None)

# display the parameters...
model.params


### (5 points) Part k)¶

Using the code examples from the lab as a guide, generate a plot of the Solow diagram demonstrating how the output, actual investment, and breakeven investment curves shift as a result of a 50% increase in $\sigma$ (i.e., the elasticity of substitution between capital and effective labor). Explain why, and in which direction, each curve shifts. Discuss the impact on the steady state levels of capital, output, and consumption per effective person.

In [ ]:
# INSERT YOUR CODE HERE!


### (5 points) Part l)¶

Using code from the lab as a guide, generate impulse response functions (IRFs) in for the Solow model following a 50% increase in the elasticity of substitution between capital and effective labor, $\sigma$. Discuss the resulting IRFS.

In [ ]:
# INSERT YOUR CODE HERE!


### (5 points) Part m:¶

Using your results from part b), derive an expression for the elasticity of output per effective person with respect to capital per effective person, $\alpha_k$, as a function of capital per effective person, $k$, and model parameters.

In [ ]:



### (5 points) Part n)¶

Differentiate your result from part m) with respect to capital per effective person, $k$, and discuss under what conditions the derivative will be positive, negative, or zero.

In [ ]:



### (5 points) Part o)¶

Use the fact that $\frac{\partial k}{\partial s} > 0$ and the result from part n) to relate the sign of $\frac{\partial \alpha_k(k)}{\partial s}$ to the elasticity of substitution between capital and effective labor.

In [ ]:



### (15 points) Part p)¶

The code in the cell below regresses $\alpha_k(k)$ on the investment's share of GDP for the U.S., $s$ during two different time periods: 1950-1970 and 1970-2011. Can you "explain" these pattern using your result from part o)?

In [ ]:
fig = plt.figure(figsize=(8,6))

# capital's share for US (pre-1970)
y = 1 - model.data.labsh.loc[:1970]

# savings rate for US (pre-1970)
x = model.data.csh_i.loc[:1970]

# regress capital's share on savings rate
res1   = pd.ols(y=y, x=x)

ax1.scatter(x, y)
ax1.set_xlabel('s', family='serif', fontsize=15)
ax1.set_ylabel(r'$\alpha_k(k)$', fontsize=15, family='serif', rotation='horizontal')
ax1.set_title('U.S. pre-1970', family='serif', fontsize=15)

# plot regression line pre-1970
grid = np.linspace(0.15, 0.25, 100)
ax1.plot(grid, res1.beta[1] + res1.beta[0] * grid, 'r')

# capital's share for US (post-1970)
y = 1 - model.data.labsh.loc[1970:]

# savings rate for US (post-1970)
x = model.data.csh_i.loc[1970:]

# regress capital's share on savings rate post-1970
res2   = pd.ols(y=y, x=x)

ax2.scatter(x, y)
ax2.set_xlabel('s', family='serif', fontsize=15)
ax2.set_title('U.S. post-1970', family='serif', fontsize=15)

# plot regression line post-1970
grid = np.linspace(0.15, 0.25, 100)
ax2.plot(grid, res2.beta[1] + res2.beta[0] * grid, 'r')

fig.suptitle("Patterns in capital's share for the U.S?", x=0.5, y=1.05,
fontsize=20, family='serif')
fig.tight_layout()
plt.show()

In [ ]:
# summary of the regression results for 1950-1970
print res1

In [ ]:
# summary of the regression results for 1970-2011
print res2

In [ ]: