This note illustrates how to code a R&D growth model in Python. The purpose of the note is to walk through Python applications, not to offer a detailed discussion of the Solow Model or to show best coding practices. The note also assumes familiarity with growth theory and a beginner experience with Python.
For a more complete and detailed discussion of Python applications see the material in Quant Econ.
A well-known result of the Solow Model is that growth is exogenous, that is, it goes unexplained by the model. Growth (in the long-run) does not depend on the savings rate, but on the growth rate of technology ($A$). The growth rate of $A$ is an assumption given to the model. Yet, Solow Model provides a number of important insights.
Endogenous growth models do not rely on "exogenous assumptions" to explain growth. One way to endogeneize growth is by adding a "production function of technology" to the model. In this sense, the R&D growth model can be intepreted as an extension of Solow model.
The model assumes two sectors: (1) production of goods and (2) production of new ideas (technology). Factors of production, capital $(K)$ and labor $(L)$, are split between these two sectors. The model is built in continuous time. Labor grows at constnat rate $n$.
Assume a Cobb-Douglas function with constant returns to scale and labor-augmenting technology.
\begin{equation} Y(t) = \left[(1-a_K)K \right]^{\alpha} \left[A(t) (1-a_L) L(t) \right]^{1-\alpha}; \; \alpha \in (0, 1) \end{equation}where $a_K$ and $a_L$ represent the share of capital and labor allocated to the production of new ideas. Therefore, $(1-a_K)$ and $(1-a_L)$ represent the share of capital and labor allocated to the production of new goods.
Similar to Solow model, the accumulation of new capital is an exogenous fixed savings rate $(s)$ of output. For simplicity, assume there is no depreciation.
\begin{equation} \dot{K(t)} = s\cdot Y(t) \end{equation}The production of new ideas follows a Cobb-Douglas function that may or may not have constant returns to scale. Whether the returns to scale in the production of new ideas is decreasing, constant, or increasing depends on the assumption of how factors of production interact in the R&D industry. Because this industry prodices ideas, the argument of "replication" (if you double inputs you double output) does not hold.
\begin{equation} \dot{A(t)} = B \left[a_K K(t) \right]^{\beta} \left[a_l L(t) \right]^{\gamma} \cdot A(t)^{\theta}; \; B>0, \beta\geq 0, \gamma\geq 0 \end{equation}The dot on top of $A$ denotes (instantaneous) change of technology and $B$ is a shift parameters. The exponent $\theta$ captures how the accumulated knowledge affects the discovery of new ideas.
If the more knowledge already discovered makes it more difficulty to discover new ideas, then $\theta$ is negative. If the accumulated knowledge makes it easier to discover new ideas, then $\theta$ is larger than one.
Solving the model implies finding how the growth rates look in the steady-state. Are growth rates constant? Or do they change with time?
Let $g_K(t)$ and $g_A(t)$ denote the growth rates of capital and technology. Then, we want to study the conditions in which $\dot{g_K(t)}=\dot{g_A(t)}=0$. To do this, we need to calculate the growth rates, and then set them constant (set the change in the growth rate equal to zero).
We can procede in three steps to find the steady-state of capital. First, input the production function (of goods) into the motion function of capital. Second, calculate the growth rate of capital. Third, set the change in the growth equal to zero and "solve for" growth of capital in terms of growth of technology. The last step will be useful to plot thedynamics of the model.
The growth rates of capital and technology are $g_K(t) \equiv \frac{\dot{K(t)}}{K(t)}$ and $g_A(t) \equiv \frac{\dot{A(t)}}{A(t)}$.
\begin{align} \dot{K(t)} &= s (1-a_K)^{\alpha} (1-a_L)^{1-\alpha} \cdot K(t)^{\alpha} A(t)^{1-\alpha} L(t)^{1-\alpha} \\[5 pt] g_K(t) &= s (1-a_K)^{\alpha} (1-a_L)^{1-\alpha} \left[\frac{A(t)L(t)}{K(t)} \right] \\[5 pt] \frac{\dot{g_K(t)}}{g_K(t)} &= (1-\alpha) \left[g_A(t) + n - g_K(t) \right] \end{align}$K$ will grow if $\left[g_A(t) + n - g_K(t) \right] > 0$. Now set $g_K(t) = 0$ and "solve for" $g_K(t)$ in terms of $g_A(t)$.
\begin{equation} g_K(t) = n + g_A(t) \end{equation}The last line is a simple linear quation with intercept $n$ and slope 1. Also, note that $g_K(t)$ does not depend on the savings rate or the other parameters of the production function (of final goods). Similar to Solow model, the growth rate of $K$ equals (in the steady state) the growth rate of population plus the growth rate of technology.
We can procede in a similar way than we did with capital to find the steady-state of technology. First, find the growth rate of technology. Then, set the change in the growth rate of technology equal to zero and solve for $g_K(t)$ in terms of $g_A(t)$.
\begin{align} \dot{A(t)} &= B \left[a_K K(t) \right]^{\beta} \left[a_l L(t) \right]^{\gamma} \cdot A(t)^{\theta} \\[5 pt] g_A(t) &= B \left[a_K K(t) \right]^{\beta} \left[a_L L(t) \right]^{\gamma} \cdot A(t)^{\theta - 1} \\[5pt] \frac{\dot{g_A(t)}}{g_A(t)} &= \beta g_K(t) + \gamma n + (\theta - 1) g_A(t) \end{align}$A$ will grow if $\left[\beta g_K(t) + \gamma n + (\theta - 1) g_A(t) \right] > 0$. Now, assuming $g_A(t) = 0$:
\begin{equation} g_K(t) = -\frac{\gamma n}{\beta} + \frac{(1-\theta)}{\beta} g_A(t) \end{equation}The last line also is a simple linear quation, with intercent $-\frac{\gamma n}{\beta}$ and slope $\frac{(1-\theta)}{\beta}$. Slope is positive if $1-\theta>0$, and the model has en equilibrium if $\frac{(1-\theta)}{\beta}>1$ (see the figure below).
Section 3 produces two equations with two unknowns.
\begin{align} g_K(t) &= n + g_A(t) \\[5pt] g_K(t) &= -\frac{\gamma n}{\beta} + \frac{(1-\theta)}{\beta} g_A(t) \end{align}Sovling this system of equation yields the (constant) growth rates in the steady-state.
\begin{align} g_K(t)^* &= n \cdot \frac{\beta + \gamma}{1-(\theta + \gamma)} + n \\[5 pt] g_A(t)^* &= n \cdot \frac{\beta + \gamma}{1-(\theta + \gamma)} \end{align}Consider the following:
To be more precise, this model can be considered a semi-endogenous growth model. Long-run growth rates depends on endogenous variables of the model, but also on the exogenous growth rate of population.
We can now plot both saddle-paths $(\dot{g_A(t)}=0, \dot{g_K(t)}=0)$ using Python. These are two linear equations.
"1|IMPORT PACKAGES"
import numpy as np # Package for scientific computing with Python
import matplotlib.pyplot as plt # Matplotlib is a 2D plotting library
"2|DEFINE PARAMETERS AND ARRAYS"
# Parameters
n = 0.02 # Growth rate of population
beta = 0.50 # Exponent of capital in the production of ideas
gamma = 1-beta # Exponent of labor in the priduction of ideas
theta = 0.25 # Exponent of technology in the production of ideas
# Arrays
gA = np.arange(0, 1, 0.01) # Build an array between 0 and 1 with step size=0.01
gK = n + gA
"3|DEFINE AND POPULATE THE SADDLE-PATH FUNCTIONS"
ss_K = n + gA # Saddle path for K
ss_A = -(gamma * n)/beta + (1-theta)/beta * gA # Saddle path for A
ss_int = -(gamma * n)/beta # Intercept of K saddle-path
"4|EQUILIBRIUM VALUES"
A_star = n * (beta + gamma)/(1 - (theta + gamma))
K_star = n + A_star
"5|PLOT THE SADDLE-PATH"
v = [0, gA.max()*0.2, -0.05, gK.max()*0.2] # Axis range
# Build figure
plt.figure(figsize=(10, 7))
plt.title("DYNAMICS OF GROWTH RATES OF CAPITAL AND TECHNOLOGY")
plt.text(-0.01 , v[3]-0.01, r'$g_K(t)$') # Manually add x-axis label
plt.text(v[1]-0.01, -0.01 , r'$g_A(t)$') # Manually add y-axis label
plt.axis(v)
plt.box(False)
plt.xticks([], []) # Hide x-axis ticks
plt.yticks([], []) # Hide y-axis ticks
plt.axvline(0, color="k") # Add vertical axis
plt.axhline(0, color="k") # Add horizontal axis
# Add both lines
plt.plot(gA, ss_A, "k-", label=r'$g_A(t)$')
plt.plot(gA, ss_K, "k-", label=r'$g_K(t)$')
# Add equilibrium mark
plt.axvline(A_star, -v[2]/(0.25), (0.05+K_star)/0.25, color="k", ls=":")
plt.axhline(K_star, 0.0 , A_star/v[1] , color="k", ls=":")
# Add reference values
plt.text(A_star, -0.01 , r'$g_A(t)^*$', horizontalalignment="center")
plt.text(-0.001, K_star, r'$g_K(t)^*$', horizontalalignment="right")
plt.text(-0.001, n , r'$n$' , horizontalalignment="right")
plt.text(-0.001, ss_int, r'$-\frac{\gamma n}{\beta}$',
horizontalalignment="right")
# Add function labels
plt.text(v[1]-0.03, v[3]-0.02, r'$\dot{g_A(t)^*}=0$')
plt.text(v[1]-0.09, v[3]-0.02, r'$\dot{g_K(t)^*}=0$')
# Build plot
plt.show()
This is a stable model. Any starting point outside the steady-state will return to equilibrium = $\left(g_A(t)^*. g_K(t)^* \right)$. Form sections 3.1 and 3.2, we can follow the model dynamics for $A$ and $K$ the following way:
\begin{align} g_A(t+1) &= g_K(t) + \dot{g_K(t)} \\[5pt] g_K(t+1) &= g_A(t) + \dot{g_A(t)}\\[15pt] g_A(t+1) &= g_K(t) + (1-\alpha) \left[g_A(t) + n - g_K(t) \right] \\[5pt] g_K(t+1) &= g_A(t) + \beta g_K(t) + \gamma n + (\theta - 1)g_A(t) \end{align}Now set four arbitrary starting points: A (blue), B (green), C (red), D (cyan). These strating points will be located in the four quadrants of the model dynamics plot. From these four stating points, the code performs 75 iterations. The plot shows that the dynamics of strating points A and C remain between the two saddle paths, while starting poitns B and D cross over a saddle path to then approach the steady-state.
"1|IMPORT PACKAGES"
import numpy as np # Package for scientific computing with Python
import matplotlib.pyplot as plt # Matplotlib is a 2D plotting library
"2|DEFINE PARAMETERS AND ARRAYS"
# Parameters
n = 0.02 # Growth rate of population
alpha = 0.60 # Exponent of capital in the production of goods
beta = 0.50 # Exponent of capital in the production of ideas
gamma = 1-beta # Exponent of labor in the priduction of ideas
theta = 0.25 # Exponent of technology in the production of ideas
# Arrays
gA = np.arange(0, 1, 0.01) # Build an array between 0 and 1 with step size=0.01
gK = n + gA
size = np.size(gA)
"3|DEFINE AND POPULATE THE SADDLE-PATH FUNCTIONS"
ss_K = n + gA # Saddle path for K
ss_A = -(gamma * n)/beta + (1-theta)/beta * gA # Saddle path for A
"4|EQUILIBRIUM VALUES"
A_star = n * (beta + gamma)/(1 - (theta + gamma))
K_star = n + A_star
"5|STABILITY DYNAMICS"
iterations = 75
" | Starting Poin A (blue)"
# Create arrays to store model dynamics
A_gA = np.zeros(iterations)
A_gK = np.zeros(iterations)
# Set arbitrary initial values
A_gA[0] = 0.015
A_gK[0] = 0.020
for j in range(1, iterations):
A_gA[j] = A_gA[j-1] + beta*A_gK[j-1] + gamma*n + (theta-1)*A_gA[j-1]
A_gK[j] = A_gK[j-1] + (1-alpha)*(A_gA[j-1] + n - A_gK[j-1])
" | Starting Poin B (green)"
# Create arrays to store model dynamics
B_gA = np.zeros(iterations)
B_gK = np.zeros(iterations)
# Set arbitrary initial values
B_gA[0] = 0.015
B_gK[0] = 0.180
for j in range(1, iterations):
B_gA[j] = B_gA[j-1] + beta*B_gK[j-1] + gamma*n + (theta-1)*B_gA[j-1]
B_gK[j] = B_gK[j-1] + (1-alpha)*(B_gA[j-1] + n - B_gK[j-1])
" | Starting Poin C (red)"
# Create arrays to store model dynamics
C_gA = np.zeros(iterations)
C_gK = np.zeros(iterations)
# Set arbitrary initial values
C_gA[0] = 0.150
C_gK[0] = 0.180
for j in range(1, iterations):
C_gA[j] = C_gA[j-1] + beta*C_gK[j-1] + gamma*n + (theta-1)*C_gA[j-1]
C_gK[j] = C_gK[j-1] + (1-alpha)*(C_gA[j-1] + n - C_gK[j-1])
" | Starting Poin D (cyan)"
# Create arrays to store model dynamics
D_gA = np.zeros(iterations)
D_gK = np.zeros(iterations)
# Set arbitrary initial values
D_gA[0] = 0.150
D_gK[0] = 0.020
for j in range(1, iterations):
D_gA[j] = D_gA[j-1] + beta*D_gK[j-1] + gamma*n + (theta-1)*D_gA[j-1]
D_gK[j] = D_gK[j-1] + (1-alpha)*(D_gA[j-1] + n - D_gK[j-1])
"5|PLOT THE SADDLE-PATH"
v = [0, gA.max()*0.2, -0.05, gK.max()*0.2] # Axis range
# Build figure
plt.figure(figsize=(10, 7))
plt.title("MODEL STABILITY")
plt.text(-0.01 , v[3]-0.01, r'$g_K(t)$') # Manually add x-axis label
plt.text(v[1]-0.01, -0.01 , r'$g_A(t)$') # Manually add y-axis label
plt.axis(v)
plt.box(False)
plt.xticks([], []) # Hide x-axis ticks
plt.yticks([], []) # Hide y-axis ticks
plt.axvline(0, color="k") # Add vertical axis
plt.axhline(0, color="k") # Add horizontal axis
# Add both lines
plt.plot(gA, ss_A, "k-", label=r'$g_A(t)$')
plt.plot(gA, ss_K, "k-", label=r'$g_K(t)$')
# Add equilibrium mark
plt.axvline(A_star, -v[2]/(0.25), (0.05+K_star)/0.25, color="k", ls=":")
plt.axhline(K_star, 0.0 , A_star/v[1] , color="k", ls=":")
# Add reference values
plt.text(A_star, -0.01 , r'$g_A(t)^*$', horizontalalignment="center")
plt.text(-0.001, K_star, r'$g_K(t)^*$', horizontalalignment="right")
# Add function labels
plt.text(v[1]-0.03, v[3]-0.02, r'$\dot{g_A(t)^*}=0$')
plt.text(v[1]-0.09, v[3]-0.02, r'$\dot{g_K(t)^*}=0$')
# Model dynamics
plt.plot(A_gA[0], A_gK[0], "bo")
plt.plot(A_gA , A_gK , "b:")
plt.plot(B_gA[0], B_gK[0], "go")
plt.plot(B_gA , B_gK , "g:")
plt.plot(C_gA[0], C_gK[0], "ro")
plt.plot(C_gA , C_gK , "r:")
plt.plot(D_gA[0], D_gK[0], "co")
plt.plot(D_gA , D_gK , "c:")
plt.text(A_gA[0]-0.005, A_gK[0] , "A", color="b")
plt.text(B_gA[0]-0.005, B_gK[0] , "B", color="g")
plt.text(C_gA[0] , C_gK[0]+0.005, "C", color="r")
plt.text(D_gA[0]+0.005, D_gK[0] , "D", color="c")
# Arrows
plt.arrow(A_gA[0], A_gK[0], 0.005, 0 , color="b")
plt.arrow(A_gA[0], A_gK[0], 0 , 0.005, color="b")
plt.arrow(B_gA[0], B_gK[0], 0.005, 0 , color="g")
plt.arrow(B_gA[0], B_gK[0], 0 ,-0.005, color="g")
plt.arrow(C_gA[0], C_gK[0], -0.005, 0 , color="r")
plt.arrow(C_gA[0], C_gK[0], 0 ,-0.005, color="r")
plt.arrow(D_gA[0], D_gK[0], -0.005, 0 , color="c")
plt.arrow(D_gA[0], D_gK[0], 0 , 0.005, color="c")
# Build plot
plt.show()