We will discuss the mathematical model and the coding section in detail in the upcoming post. I also need to correct some errors I made in the code section. For example, I made a couple of errors in the ‘while’ loop! There are more, though the code works just fine. So, for now, this is just a preview.
Still, for students of economics, even in their early years, this should look quite familiar (except for the coding part). Exposure to coding and applications of other programming languages and/or statistical software’s in economics depends on the program/school you attend. Nonetheless, you can start learning Python now with QuantEcon. QuanEcon, founded by Thomas J. Sargent of New York University and John Stachurski of the Australian National University, is an excellent source for students/practitioners who wish to make the jump. Visit QuantEcon via the following link https://quantecon.org/index.html
In [1]:
import numpy as np
import matplotlib.pyplot as plt
%matplotlib inline
import sympy as sp
from sympy import pretty_print, init_printing, solve
from textwrap import dedent
from IPython.display import display
from ipywidgets import interactive
![\[Y_t = zK_t^\alpha L_t^{1-\alpha}\]](https://i0.wp.com/quicklatex.com/cache3/13/ql_208eac90ba150708eb4057461e735613_l3.png?resize=105%2C20&ssl=1 "Rendered by QuickLaTeX.com")
![\[y_t = zk_t^\alpha\]](https://i0.wp.com/quicklatex.com/cache3/bd/ql_8dfbb67fc90bb3017cf8b953051efdbd_l3.png?resize=61%2C17&ssl=1 "Rendered by QuickLaTeX.com")
In [2]:
# SOLVING FOR THE STEADY STATE VALUES
## variables:
y, k, c = sp.symbols('y, k, c')
## parameters:
z, s, d, n, alpha = sp.symbols('z, s, d, n, alpha')
## capital per capita
y = z*k**alpha
kss = k - (s*y + (1-d)*k)/(1+n)
init_printing(pretty_print=True)
display("Steady state per capita capital", solve(kss, k))
init_printing(pretty_print=False)
display(solve(kss, k))
![\[\left [ \left(\frac{s z}{d + n}\right)^{- \frac{1}{\alpha - 1}}\right ]\]](https://i0.wp.com/quicklatex.com/cache3/92/ql_61c81bc8977d976e4a383995cd636092_l3.png?resize=115%2C53&ssl=1 "Rendered by QuickLaTeX.com")
In [3]:
## output per capita
k = (s*z/(d + n))**(-1/(alpha - 1))
yss = y - z*k**alpha
init_printing(pretty_print=True)
display("Steady state per capita output" ,solve(yss, y))
init_printing(pretty_print=False)
display(solve(yss, y))
![\[\left [ z \left(\left(\frac{s z}{d + n}\right)^{- \frac{1}{\alpha - 1}}\right)^{\alpha}\right ]\]](https://i0.wp.com/quicklatex.com/cache3/99/ql_80a8d8626ec45fc812b7fa9bc5f21d99_l3.png?resize=162%2C53&ssl=1 "Rendered by QuickLaTeX.com")
In [4]:
## consumption per capita
y = z*((s*z/(d + n))**(-1/(alpha - 1)))**alpha
CSS = c - (1-s)*y
init_printing(pretty_print=True)
display("Steady state per capita output", solve(CSS, c))
init_printing(pretty_print=False)
display(solve(CSS, c))
![\[\left [ z \left(- s + 1\right) \left(\left(\frac{s z}{d + n}\right)^{- \frac{1}{\alpha - 1}}\right)^{\alpha}\right ]\]](https://i0.wp.com/quicklatex.com/cache3/f8/ql_725afcd83cbbef72fb406d3c5465def8_l3.png?resize=227%2C53&ssl=1 "Rendered by QuickLaTeX.com")
In [5]:
class SolowSwan():
""" SolowSwan Model Simulation """
def __init__(self, z = 1.0, alpha = 0.33, L = 10, K = 10,
s = 0.2, d = 0.05, n = 0.02):
""" Model Initialization """
# we admit only non-negative values
positive_array = np.array([z, alpha, L, K, s, d, n])
zeros = np.zeros(len(positive_array))
if np.any(np.less_equal(positive_array, zeros)):
raise ValueError("Only positive values admitted.")
# create class attributes for SolowSwan
self.z = z
self.alpha, self.beta = alpha, 1-alpha
self.L, self.K = L, K
self.s, self.d, self.n = s, d, n
def __repr__(self):
" Basic Information about the model "
return self.__str__()
def __str__(self):
message = """Solow-Swan model with the following features:
- Capital (K) = {K}
- Labor (L) = {L}
- Total factor productivity = {z}
- Income share of Capital (alpha) = {alpha}
- Income share of Labor (beta) = {beta}
- Saving rate (s) = {s}
- Depreciation rate (d) = {d}
- Population growth rate (n) = {n}"""
return dedent(message.format(K = self.K, L = self.L, z =self.z,
alpha = self.alpha, beta = self.beta,
s = self.s, d = self.d, n = self.n))
def TotalOutput(self):
""" Aggregate Output """
z, alpha, beta = self.z, self.alpha, self.beta
return z*(self.K**alpha)*(self.L**self.beta)
def perCapOutput(self, k=None):
""" Per Capit Output """
if k is None:
k = self.K/self.L
return self.z*k**self.alpha
def ss_analytic(self):
""" Steady state values calculated using the derived formulae above """
z, alpha, beta, s, d, n = self.z, self.alpha, self.beta, self.s, self.d, self.n
kss = (s*z/(d + n))**(-1/(alpha - 1))
yss = z*kss**alpha
css = (1-s)*yss
# store values as attributes
self.kss, self.yss, self.css = kss, yss, css
return kss, yss, css
def ss_numerical(self, k0 = None, tol = 1e-4, maxiter=200):
""" Compute the steady state values via an iterative procedure """
z, alpha, beta, s, d, n = self.z, self.alpha, self.beta, self.s, self.d, self.n
if k0 is None:
k0 = self.K/self.L
y0 = self.perCapOutput(k0)
error = 1
i = 0
while i < maxiter and error > tol:
knext = (s*y0 + (1-d)*k0)/(1+n)
error = np.abs(knext-k0)
if error < tol:
break
k0 = knext
y0 = self.perCapOutput(k0)
i += 1
kss = knext
yss = y0
css = (1-s)*yss
return kss, yss, css
def next_value_generator(self, T, k0=None):
""" Create a generator for the transition paths of k, y and c"""
z, alpha, beta, s, d, n = self.z, self.alpha, self.beta, self.s, self.d, self.n
## if k0 is not provided, then adopt the ss value of k
if k0 is None:
k0 = self.K/self.L
## set initial values
k_next = k0
y_next = self.perCapOutput(k_next)
c_next = (1-s)*y_next
for t in range(T):
yield k_next, y_next, c_next
k_next = (s*y_next + (1-d)*k_next)/(1+n)
y_next = self.perCapOutput(k=k_next)
c_next = (1-s)*y_next
def transition_paths(self, **kwargs):
""" Simulate transition paths using the next_value_generator module"""
z, alpha, beta, s, d, n = self.z, self.alpha, self.beta, self.s, self.d, self.n
paths = np.array(list(self.next_value_generator(**kwargs)))
self.k_path = paths[:,0]
self.y_path = paths[:,1]
self.c_path = paths[:,2]
def transition_paths_alt(self, T=100, k0=None):
""" Simulate transition paths using an alternative method """
z, alpha, beta, s, d, n = self.z, self.alpha, self.beta, self.s, self.d, self.n
## if k0 is not provided, then adopt the ss value of k
if k0 is None:
k0 = self.K/self.L
## NOTE: I can also just use paths = np.zeros([T, 3]) or np.empty([T, 3]) to store simulated values
## where the 3 columns are for k, y, c
## and T is just time T
## But, creating a dictionary as below is more intuitive in my opinion
paths = {'k': np.zeros(T), 'y': np.zeros(T), 'c': np.zeros(T)}
# intialize value
paths['k'][0] = k0
paths['y'][0] = self.perCapOutput(k0)
paths['c'][0] = (1-s)*paths['y'][0]
for t in range(1, T):
paths['k'][t] = (s*paths['y'][t-1] + (1-d)*paths['k'][t-1])/(1+n)
paths['y'][t] = self.perCapOutput(paths['k'][t])
paths['c'][t] = (1-s)*paths['y'][t]
self.paths = paths
return paths
In [6]:
list(range(1,10))
Out[6]:
In [7]:
np.not_equal([1, 1],[0.,0.])
Out[7]:
In [8]:
list(range(1,10))
Out[8]:
In [9]:
m = SolowSwan(z=3)
In [10]:
m.transition_paths(T=100)
m.k_path
Out[10]:
In [11]:
m.ss_analytic()
Out[11]:
In [12]:
m.ss_numerical()
Out[12]:
In [13]:
m.transition_paths(T=10, k0=None)
display(m.k_path, m.y_path, m.c_path)
In [14]:
m.transition_paths_alt(T=10, k0 = None)
Out[14]:
In [15]:
m.paths['k']
Out[15]:
In [16]:
## Visualization using m.transition_paths()
m = SolowSwan()
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(15, 20))
ss = m.ss_analytic() # this generates 3 ss values for m instance
kss_path = np.array([m.kss,]*50)
yss_path = np.array([m.yss,]*50)
css_path = np.array([m.css,]*50)
plt.rcParams['axes.xmargin'] = 0.01
plt.rcParams['axes.ymargin'] = 0.01
plt.style.use('seaborn-whitegrid')
time = 100
for z in [0.1, 0.5, 1.0, 1.5, 2]:
m.z = z
m.transition_paths(T=time, k0=m.kss);
## ax.plot(np.array([kss_path, m.k_path]).flatten())
k_series = np.concatenate((kss_path, m.k_path), axis=0)
y_series = np.concatenate((yss_path, m.y_path), axis=0)
c_series = np.concatenate((css_path, m.c_path), axis=0)
ax1.plot(k_series, label=f"
, z = {z}", alpha=3)
ax2.plot(y_series, label=f"
, z = {z}", alpha=3)
ax3.plot(c_series, label=f"
, z = {z}", alpha=3)
ax_list = (ax1, ax2, ax3)
ylabel_list = ('k', 'y', 'c')
for ax, i in zip(ax_list, ylabel_list):
ax.set_xlabel("Time, t")
ax.set_ylabel("f
")
ax.legend(loc="upper left", fontsize=10)
ax.set_ylim(0,20)
ax.set_xticks(np.arange(0,time+50,5))
ax1.set_title("Capital (equilibrium) Transition Path")
ax2.set_title("Output (equilibrium) Transition Path")
ax3.set_title("Consumption (equilibrium) Transition Path")
plt.show()
In [17]:
## Visualization using m.transition_paths_alt
m = SolowSwan()
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(15, 20))
ss = m.ss_analytic() # this generates 3 ss values for m instance
kss_path = np.array([m.kss,]*50)
yss_path = np.array([m.yss,]*50)
css_path = np.array([m.css,]*50)
plt.rcParams['axes.xmargin'] = 0.01
plt.rcParams['axes.ymargin'] = 0.01
plt.style.use('seaborn-whitegrid')
time = 100
for z in [0.1, 0.5, 1.0, 1.5, 2]:
m.z = z
m.transition_paths_alt(T=time, k0=m.kss);
## ax.plot(np.array([kss_path, m.k_path]).flatten())
k_series = np.concatenate((kss_path, m.paths['k']), axis=0)
y_series = np.concatenate((yss_path, m.paths['y']), axis=0)
c_series = np.concatenate((css_path, m.paths['c']), axis=0)
ax1.plot(k_series, label=f", z = {z}", alpha=3)
ax2.plot(y_series, label=f", z = {z}", alpha=3)
ax3.plot(c_series, label=f", z = {z}", alpha=3)
ax_list = (ax1, ax2, ax3)
ylabel_list = ('k', 'y', 'c')
for ax, i in zip(ax_list, ylabel_list):
ax.set_xlabel("Time, t")
ax.set_ylabel("f")
ax.legend(loc="upper left", fontsize=10)
ax.set_ylim(0,20)
ax.set_xticks(np.arange(0,time+50,5))
ax1.set_title("Capital (equilibrium) Transition Path")
ax2.set_title("Output (equilibrium) Transition Path")
ax3.set_title("Consumption (equilibrium) Transition Path")
plt.show()
In [18]:
def plotting(z=1, alpha=0.33, s=0.2, d=0.05, n=0.02):
m = SolowSwan()
fig, (ax1, ax2, ax3) = plt.subplots(3, 1, figsize=(8,8))
ss = m.ss_analytic() # this generates 3 ss values for m instance
kss_path = np.array([m.kss,]*50)
yss_path = np.array([m.yss,]*50)
css_path = np.array([m.css,]*50)
plt.rcParams['axes.xmargin'] = 0.01
plt.rcParams['axes.ymargin'] = 0.01
plt.style.use('seaborn-whitegrid')
# introduce shock
m.z, m.alpha = z, alpha
m.s, m.d, m.n = s, d, n
m.transition_paths_alt(T=200, k0=m.kss);
## ax.plot(np.array([kss_path, m.k_path]).flatten())
k_series = np.concatenate((kss_path, m.paths['k']), axis=0)
y_series = np.concatenate((yss_path, m.paths['y']), axis=0)
c_series = np.concatenate((css_path, m.paths['c']), axis=0)
ax1.plot(k_series, label=f", z = {z}", alpha=3)
ax2.plot(y_series, label=f", z = {z}", alpha=3)
ax3.plot(c_series, label=f", z = {z}", alpha=3)
ax_list = (ax1, ax2, ax3)
ylabel_list = ('k', 'y', 'c')
for ax, i in zip(ax_list, ylabel_list):
ax.set_xlabel("Time, t")
ax.set_ylabel("f")
ax.legend(loc="upper left", fontsize=10)
ax.set_ylim(0,20)
ax1.set_title("Capital (equilibrium) Transition Path")
ax2.set_title("Output (equilibrium) Transition Path")
ax3.set_title("Consumption (equilibrium) Transition Path")
plt.show()
In [19]:
interactive_plot = interactive(plotting, z=(0.1,3,0.1),
alpha=(0.1, 1, 0.1),
L=(5,100,1),
K=(5,500,10),
s=(0.1,0.9,0.1),
d=(0.01,0.9,0.1),
n=(0.01,0.3,0.01))
interactive_plot
What is interesting is that as we keep increasing s holding everything else constant, while output and capital continue to move towards higher values, the consumption side is quite different. This suggests that a society with high saving rate is not necessarily one with high standard of living if we make the very basic assumption that the standard of living depends only on consumption (of whatever).
In [ ]: