An Exciting New Project (Preview)

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


Solow-Swan Model¶

where:

(1)

Per capita Output Function and Transition Equations¶

where:

(2)

Next, we use sympy to solve analytically for the steady state values of output, capital and consumption per capita:

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))

'Steady state per capita capital'

[(s*z/(d + n))**(-1/(alpha - 1))]
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))

'Steady state per capita output'

[z*((s*z/(d + n))**(-1/(alpha - 1)))**alpha]
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))

'Steady state per capita output'

[z*(-s + 1)*((s*z/(d + n))**(-1/(alpha - 1)))**alpha]
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]:
[1, 2, 3, 4, 5, 6, 7, 8, 9]
In [7]:
np.not_equal([1, 1],[0.,0.])

Out[7]:
array([ True,  True])
In [8]:
list(range(1,10))

Out[8]:
[1, 2, 3, 4, 5, 6, 7, 8, 9]
In [9]:
m = SolowSwan(z=3)

In [10]:
m.transition_paths(T=100)
m.k_path

Out[10]:
array([ 1.        ,  1.51960784,  2.09066104,  2.69750111,  3.32852554,
3.97486797,  4.62961207,  5.28728921,  5.94353882,  6.59486839,
7.23847763,  7.87212591,  8.49403   ,  9.10278357,  9.69729289,
10.27672474, 10.84046374, 11.38807715, 11.91928544, 12.43393774,
12.93199104, 13.41349267, 13.87856531, 14.32739431, 14.76021678,
15.17731223, 15.57899467, 15.96560566, 16.33750845, 16.69508291,
17.03872122, 17.36882408, 17.68579758, 17.99005044, 18.2819917 ,
18.56202875, 18.83056566, 19.08800176, 19.33473048, 19.57113841,
19.79760443, 20.01449915, 20.22218434, 20.42101255, 20.61132681,
20.79346045, 20.96773691, 21.13446972, 21.29396243, 21.44650871,
21.59239236, 21.73188745, 21.86525844, 21.99276035, 22.11463892,
22.23113083, 22.34246391, 22.44885731, 22.55052179, 22.64765992,
22.74046632, 22.8291279 , 22.91382408, 22.99472709, 23.07200214,
23.14580768, 23.21629566, 23.28361172, 23.34789542, 23.40928051,
23.46789507, 23.52386179, 23.57729814, 23.62831658, 23.67702473,
23.72352562, 23.7679178 , 23.81029557, 23.85074915, 23.8893648 ,
23.92622504, 23.96140876, 23.99499139, 24.02704504, 24.05763865,
24.08683812, 24.1147064 , 24.14130368, 24.16668745, 24.19091268,
24.21403185, 24.23609511, 24.25715037, 24.2772434 , 24.29641791,
24.31471565, 24.33217649, 24.34883852, 24.3647381 , 24.37990995])
In [11]:
m.ss_analytic()

Out[11]:
(24.69543792333115, 8.643403273165902, 6.914722618532721)
In [12]:
m.ss_numerical()

Out[12]:
(24.692713458647063, 8.64308858586437, 6.914470868691495)
In [13]:
m.transition_paths(T=10, k0=None)
display(m.k_path, m.y_path, m.c_path)

array([1.        , 1.51960784, 2.09066104, 2.69750111, 3.32852554,
3.97486797, 4.62961207, 5.28728921, 5.94353882, 6.59486839])
array([3.        , 3.44423407, 3.82661572, 4.16234998, 4.46133032,
4.73039872, 4.97451762, 5.19742423, 5.40201939, 5.59061105])
array([2.4       , 2.75538726, 3.06129258, 3.32987998, 3.56906425,
3.78431897, 3.97961409, 4.15793938, 4.32161551, 4.47248884])
In [14]:
m.transition_paths_alt(T=10, k0 = None)

Out[14]:
{'k': array([1.        , 1.51960784, 2.09066104, 2.69750111, 3.32852554,
3.97486797, 4.62961207, 5.28728921, 5.94353882, 6.59486839]),
'y': array([3.        , 3.44423407, 3.82661572, 4.16234998, 4.46133032,
4.73039872, 4.97451762, 5.19742423, 5.40201939, 5.59061105]),
'c': array([2.4       , 2.75538726, 3.06129258, 3.32987998, 3.56906425,
3.78431897, 3.97961409, 4.15793938, 4.32161551, 4.47248884])}
In [15]:
m.paths['k']

Out[15]:
array([1.        , 1.51960784, 2.09066104, 2.69750111, 3.32852554,
3.97486797, 4.62961207, 5.28728921, 5.94353882, 6.59486839])
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 [ ]: