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

Aggregate Output Function and Transition Equations

    \[Y_t = zK_t^\alpha L_t^{1-\alpha}\]

where:

(1)   \begin{align*} K_{t+1} &= sY_t - (n+d)K_t \\ C_{t} &= (1-s)Y_t \\ N_{t+1} &= (1+n)N_t \end{align*}

Per capita Output Function and Transition Equations

    \[y_t = zk_t^\alpha\]

where:

(2)   \begin{align*} k_{t+1} &= \frac{sy_t + (1-d)k_t}{1+n} \\ c_{t} &= (1-s)y_t \\ \end{align*}

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'

    \[\left [ \left(\frac{s z}{d + n}\right)^{- \frac{1}{\alpha - 1}}\right ]\]

[(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'

    \[\left [ z \left(\left(\frac{s z}{d + n}\right)^{- \frac{1}{\alpha - 1}}\right)^{\alpha}\right ]\]

[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'

    \[\left [ z \left(- s + 1\right) \left(\left(\frac{s z}{d + n}\right)^{- \frac{1}{\alpha - 1}}\right)^{\alpha}\right ]\]

[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"k_t, z = {z}", alpha=3)
    ax2.plot(y_series, label=f"y_t, z = {z}", alpha=3)
    ax3.plot(c_series, label=f"c_t, 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</span><span class="si">{i}</span><span class="s2">_t")
    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"k_t, z = {z}", alpha=3)
    ax2.plot(y_series, label=f"y_t, z = {z}", alpha=3)
    ax3.plot(c_series, label=f"c_t, 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</span><span class="si">{i}</span><span class="s2">_t")
    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"k_t, z = {z}", alpha=3)
    ax2.plot(y_series, label=f"y_t, z = {z}", alpha=3)
    ax3.plot(c_series, label=f"c_t, 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</span><span class="si">{i}</span><span class="s2">_t")
        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 [ ]:
 

Leave a Reply