#!/usr/bin/env python
# coding: utf-8

# # <font color="880000">Ancient Economies: A Malthusian Model</font>

# ## <font color="880000"> The Basic Setup </font>
# 
# The basic setup model has production (or output) $ Y $ as a function of the capital stock $ K $, the labor force $ L $, the efficiency of labor $ E $, and the decreasing-returns parameter $ \alpha $:
# 
# >(1) $ Y = K^{\alpha}(EL)^{1-\alpha} $ :: production
# 
# >(2) $ \ln(Y) = \alpha\ln(K) + (1-\alpha)\left(\ln(E)+\ln(L)\right) $ :: log production
# 
# Growth rates of the capital stock, the labor force, the efficiency of labor, and of production as functions of the savings-investment rate $ s $, the population and labor force growth rate $ n $, the efficiency-of-labor growth rate $ g $, and the depreciation rate $ \delta $:
# 
# >(3) $ \frac{dK/dt}{K} = \frac{d\ln(K)}{dt} = g_k = \frac{sY}{K} - \delta $ :: the proportional rate of growth of the capital stock
# 
# >(4) $ \frac{dL/dt}{L} = \frac{d\ln(L)}{dt} = n $ :: the proportional rate of growth of the labor force and population
# 
# >(5) $ \frac{dE/dt}{E} = \frac{d\ln(E)}{dt}  = g $ :: the proportional rate of growth of the efficiency of labor
# 
# >(6) $ \frac{d\ln(Y)}{dt} = g_{n+y} = \alpha g_k + (1-\alpha)n + (1-\alpha)g  $ :: the proportional rate of growth of total production
# 
# And define the capital-output ratio:
# 
# >(7) $ \kappa = \frac{K}{Y} $
# 
# >(8) $ \ln(\kappa) = \ln(K) - \ln(Y) $
# 
# >(9) $ \frac{d\ln(\kappa)}{dt} = g_\kappa = g_k - g_y $ :: proportional growth rate of the capital-output ratio
# 
# &nbsp;

# ### <font color="000088"> Determining the Equilibrium Capital-Output Ratio $ \kappa^* $ </font>
# 
# Now fix the parameters $ \alpha, \beta, s, h $, and look for a situation in which $ g_\kappa = 0 $: in which the capital stock $ K $ and production $ Y $ are growing at the same rate so that the capital-output ratio $ g_\kappa $ is constant.
# 
# Starting from (9) and substituting:
# 
# >(10) $ g_\kappa = g_k - \left( \alpha g_k + (1-\alpha)n + (1-\alpha)g \right)   $
# 
# >(11) $ g_\kappa - g_{y+n} = (1-\alpha)\left( g_k - n - g \right) $
# 
# >(12) $ g_\kappa = (1-\alpha)\left( \frac{sY}{K} - \delta - n - g \right) $
# 
# >(13) $ n + g + \delta = \frac{sY}{K} $ whenever $ g_\kappa = 0 $
# 
# >(14) $ \kappa = \frac{K}{Y} = \frac{s}{n + g + \delta} $ whenever $ g_\kappa = 0 $
# 
# So we define the _steady-state growth equilibrium_ capital-output ratio:
# 
# >(15) $ \kappa^*_{[s,n,g,\delta]} = \left( \frac{K}{Y} \right)^* = \frac{s}{n + g + \delta} $
# 
# &nbsp;

# ### <font color="000088"> Determining Steady-State Growth-Path Production per Worker </font>
# 
# Recall 2:
# 
# >(2) $ \ln(Y) = \alpha\ln(K) + (1-\alpha)\left(\ln(E)+\ln(L)\right) $
# 
# from that derive:
# 
# >(16) $ \ln(Y) = \alpha\left(\ln(\kappa) + \ln(Y) \right) + (1-\alpha)\left(\ln(E)+\ln(L)\right) $
# 
# >(17) $ (1-\alpha)\ln(Y) = \alpha\ln(\kappa) + (1-\alpha)\left(\ln(E)+\ln(L)\right) $
# 
# >(18) $ \ln(Y) = \left( \frac{\alpha}{1-\alpha} \right)\ln(\kappa) + \ln(L) + \ln(E) $
# 
# >(19) $ \ln \left( \frac{Y}{L} \right) = \left( \frac{\alpha}{1-\alpha} \right)\ln(\kappa) + \ln(E) $
# 
# These $ \alpha/(1-\alpha) $ terms are getting annoying:
# 
# >(20) $ \theta = \frac{\alpha}{1-\alpha} $
# 
# >(21) $ \ln \left( \frac{Y}{L} \right) = \theta\ln(\kappa) + \ln(E) $
# 
# And so we define steady-state growth-path production-per-worker as:
# 
# >(22) $ \ln \left( \frac{Y}{L} \right)^* = \theta \ln(\kappa^*) =+ \ln(E)  $ 
# 
# >(23) $ \left(\frac{Y}{L}\right)^* = \left(\kappa^*\right)^\theta E $
# 
# >(24) $ \left(\frac{Y}{L}\right)^* = \left( \frac{s}{n+g+\delta} \right)^\theta E $
# 
# >(25) $ \frac{d}{dt} \left(\frac{Y}{L}\right)^* = g $
# 
# Along the steady-state growth path, production per worker $ Y/L $, capital per worker $ K/L $, and the efficiency of labor $ E $ both grow at the proportional rate $ g $; the population and labor force $ L $ grows at the proportional rate $ n $; total production $ Y $ and the capital stock $ K $ grow at the proportional rate $ n + g $; and the capital-output ratio $ K/L $ is constant.
# 
# &nbsp;

# ## <font color="880000"> Population, Resource Scarcity, and the Efficiency of Labor </font>
# 
# Now let's complicate the determinants of the efficiency of labor. Let's make its growth a function of the rate $ h $ at which economically useful ideas are generated, but also of the rate of population and labor force growth $ n $ because a higher population makes resources per capita scarce, as determined by an effect-of-resource scarcity parameter $ \gamma $:
# 
# >(26) $ \frac{dE/dt}{E} = \frac{d\ln(E)}{dt} = g = h - \frac{n}{\gamma} $
# 
# Thus:
# 
# >(27) $ \frac{d}{dt} \left(\frac{Y}{L}\right)^* = 0 $ whenever $ h - \frac{n}{\gamma} = 0 $
# 
# >(28) $ n^{*mal} = \gamma h $ is the population growth rate at which $ \frac{d}{dt} \left(\frac{Y}{L}\right)^* = 0 $
# 
# When population is growing at the rate $ n^{*mal} $, the efficiency of labor—and thus the steady-state growth-path level of production per worker $ Y/L $—is constant.
# 
# &nbsp;

# ### <font color="000088"> Determinants of Population and Labor Force Growth </font>
# 
# Now let's make the rate of growth of the population and labor force depend on the level of prosperity $ y = Y/L $; on the "subsistence" standard of living for necessities $y^{sub} $; and also on the fraction $ 1/\phi $ of production that is devoted to necessities, not conveniences and luxuries, and thus enters into reproductive and survival fitness:
# 
# >(29) $ \frac{dL/dt}{L} = \frac{d\ln(L)}{dt} = n = \beta  \left( \frac{y}{\phi y^{sub}}-1 \right) $
# 
# Then for population to be growing at its Malthusian rate:
# 
# >(30) $ \gamma h = \beta \left(\frac{1}{\phi}\right) \left( \frac{y}{y^{sub}}- \phi \right)  $
# 
# >(31) $ \frac{\phi \gamma h}{\beta} = \left( \frac{y}{y^{sub}}- \phi \right) $
# 
# >(32) $ \frac{y}{y^{sub}} = \phi \left( 1 + \frac{\gamma h}{\beta} \right) $
# 
# >(33) $ y^{*mal} = \phi y^{sub} \left( 1 + \frac{ n^{*mal}}{\beta}\right) = \phi y^{sub} \left( 1 + \frac{ \gamma h}{\beta}\right) $
# 
# &nbsp;

# ### <font color="000088"> The Full Equilibrium </font>
# 
# Thus we have our equilibrium for the pre-industrial Malthusian economy:
# 
# * Start with the rate $ h $ at which new economically-useful ideas are being generated and with the responsiveness $ \beta $ of population growth to increased prosperity.
# 
# * From those derive the Malthusian rate of population growth: $ n^{*mal} = \gamma h $
# 
# * Then the Malthusian standard of living is: $ y^{*mal} = \phi y^{sub} \left( 1 + \frac{ \gamma h}{\beta}\right) $
# 
# From 26, we can determine the log level $ E $ of the efficiency of labor:
# 
# >(32) $ \ln(E) = \ln(H) - \frac{\ln(L)}{\gamma} $
# 
# Recall (24):
# 
# >(24) $ y^* = \left( \frac{s}{n+g+\delta} \right)^\theta E $
# 
# Then:
# 
# >(33) $ y^{*mal} = \left( \frac{s}{\gamma h +\delta} \right)^\theta E $
# 
# >(34) $ \ln(\phi) + \ln\left( y^{sub} \right) + \ln\left(1 + \frac{\gamma h}{\beta} \right) = \theta \ln(s) - \theta \ln(\gamma h +\delta) + \ln(E) $
# 
# >(35) $ \ln(y^{sub}) = \theta \ln(s) - \theta \ln(\gamma h +\delta) + \ln(H_t) - \frac{\ln(L_t)}{\gamma} - \ln(\phi) - \ln\left(1 + \frac{\gamma h}{\beta} \right)  $
# 
# >(36) $ \frac{\ln(L_t)}{\gamma} = \theta \ln(s) - \theta \ln(\gamma h +\delta) + \ln(H_t) - \ln(\phi) - \ln( y^{sub}) -ln\left(1 + \frac{\gamma h}{\beta} \right) $
# 
# Thus the population and labor force in the full Malthusian equilibrium will be:
# 
# >(37) $ \ln(L_t^{*mal}) = \gamma \left[ \theta \ln(s) - \theta \ln(\gamma h +\delta) + \ln(H_t) - \ln(\phi) -  \ln( y^{sub}) -ln\left(1 + \frac{\gamma h}{\beta} \right) \right] $
# 
# It might be worthwhile decomposing this into terms depending on: (1) the level and rate of growth of the stock of _ideas_; (2) the rule of law and thrift that drive investment; (3) luxuries (including urbanization and elites); and demography in form of sociological (and law-and-order) determinants of "subsistence":
# 
# >(38) $ \ln(L_t^{*mal}) = \gamma \left[ \ln(H_t)  - \theta \ln(\gamma h +\delta) -ln\left(1 + \frac{\gamma h}{\beta} \right) \right] +  \gamma \left[ \theta \ln(s) \right] -  \gamma \left[ \ln(\phi) \right]  - \gamma \left[ \ln( y^{sub}) \right] $
# 
# And recall (31) for the full Malthusian equilibrium standard of living:
# 
# >(31) $ y^{*mal} = \phi y^{sub} \left( 1 + \frac{ \gamma h}{\beta}\right) $
# 
# Plus for the rate of population growth:
# 
# >(28) $ n^{*mal} = \gamma h $

# Production per worker and thus prosperity is determined by (a) true subsistence, (b) the wedge between prosperity and reproductive fitness produced by spending on conveniences and luxuries that do not impact reproductive success, and (c) the wedge above subsistence needed to generate population growth consonant with the advance of knowledge and population pressure's generation of resource scarcity.
# 
# The level of population is:
# 
# * raised by increases in the level of economically-useful knowledge $ H_t $
# * raised by increases in the importance $ \theta $ of capital accumulation for production
# * raised by increases in the share of production saved and invested $ s $
# * lowered by increases in the rate of depreciation $ \delta $, in the rate of increase of knowledge $ h $, and in the parameter $ \gamma $ capturing how much resource scarcity reduces the efficiency of labor $ E $
# * lowered by an increase in the wedge $ \phi $ between prosperity and subsistence created by spending on luxuries and conveniences (middle-class luxuries that do not affect reproduction, but also the "luxury" of having an upper class, and the additional luxuries of living in cities and having trade networks that can spread plagues).
# * lowered by an increase in basic biological subsistence requirements $ y^{sub} $ (produced by, say, Western European marriage patterns or lineage-family control of reproduction by clan heads).
# * lowered by an increase in the wedge $ \gamma h / \beta $ between basic subsistence and spending on necessities needed to generate population growth consonant with the advance of knowledge and population pressure's generation of resource scarcity
# 
# 
# 
# &nbsp;

# In[1]:


# DEFINING CLASS MALTHUSIAN
#
# kept in delong_classes

import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')

import numpy as np

class malthusian:
    
    """
    Implements the Malthusian Model with:
    
    1. population growth 
        n =  β*(y/(ϕ ysub)-1)
        
    2. growth of efficiency-of-labor 
        g = h-n/γ
    """
    def __init__(self,
                 L = 1,               # initial labor force
                 E = 1/3,             # initial efficiency of labor
                 K = 3.0,             # initial capital stock
                 
                               # determinants of n (population growth):
                 β = 0.025,           # responsiveness of population growth to increased prosperity.
                 ϕ = 1,               # luxuries parameter
                 ysub = 1,            # subsistence level
                 
                               # determinants of g(efficiency-of-labor growth):
                 h = 0,               # rate at which useful ideas are generated
                 γ = 2.0,             # effect-of-resource scarcity parameter  
                 
                 s = 0.15,            # savings-investment rate
                 α = 0.5,             # orientation-of-growth-toward-capital parameter
                 δ = 0.05,            # deprecation rate on capital parameter
                ):
        self.L, self.E, self.K, self.h, self.γ, self.s, self.α, self.δ = L, E, K, h, γ, s, α, δ 
        self.β, self.ϕ, self.ysub = β, ϕ, ysub
        
        # production (or output)
        self.Y = self.K**self.α*(self.E*self.L)**(1-self.α) 
        self.y = self.Y/self.L
        
        # capital-output ratio 
        self.κ = self.K/self.Y    
        
        # population growth 
        self.n = self.β*((self.y/(self.ϕ*self.ysub)) - 1)
        
        # growth rate of efficiency-of-labor
        self.g = self.h-self.n/self.γ
        
        # store initial data
        self.initdata = vars(self).copy()
    
    def update(self):
        # unpack parameters
        K, s, Y, δ, L, n, E, g, α =self.K, self.s, self.Y, self.δ, self.L, self.n, self.E, self.g, self.α
        β, ϕ, ysub, h, γ = self.β, self.ϕ, self.ysub, self.h, self.γ
        
        #update variables        
        K = s*Y + (1-δ)*K
        L = L*np.exp(n)
        E = E*np.exp(g)
        Y = K**α*(E*L)**(1-α)
        y = Y/L
        κ = K/Y
        n = β*(y/(ϕ*ysub)-1)
        g = h-n/γ
                
        #store variables
        self.K, self.s, self.Y, self.δ, self.L, self.n, self.E, self.g, self.α = K, s, Y, δ, L, n, E, g, α
        self.κ, self.y = κ, y
        
    def gen_seq(self, t, var = 'κ', init = True, log = False):
        "Generate and return time series of selected variable. Variable is κ by default."
        
        path = []
        
        # initialize data 
        if init == True:
            for para in self.initdata:
                 setattr(self, para, self.initdata[para])

        for i in range(t):
            path.append(vars(self)[var])
            self.update()
        
        if log == False:
            return path
        else:
            return np.log(np.asarray(path))

    def steady_state(self, disp = True):
        "Calculate variable values in the steady state"
        #unpack parameters
        s, γ, h, δ, ϕ, ysub, β, α= self.s, self.γ, self.h, self.δ, self.ϕ, self.ysub, self.β, self.α
        
        self.mal_κ = s/(γ*h+δ)
        # malthusian rate of population growth
        self.mal_n = γ*h
        # malthusian standard of living
        self.mal_y = ϕ*(ysub+γ*h/β)
        self.mal_E = self.mal_y*((γ*h+δ)/s)**(α/(1-α))
        
        if display == True:
            return(f'steady-state capital-output ratio κ: {self.mal_κ:.2f}')
            return(f'Malthusian rate of population growth n: {self.mal_n: .2f}')
            return(f'Malthusian standard of living y: {self.mal_y:.2f}')
            return(f'steady-state efficiency-of-labor E: {self.mal_E:.2f}') 
        else: 
            return(self.mal_κ,self.mal_n,self.mal_y,self.mal_E)


# In[2]:


# TESTING THE MALTHUSIAN CLASS
# 
# No idea stock growth (h = 0)

import numpy as np
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')

import delong_classes

T = 1000

m_base = delong_classes.malthusian(K=3.0)
m_base.scenario = "base scenario"
m_alt = delong_classes.malthusian(K=3.0)
m_alt.scenario = "alt scenario"

figcontents = {
        (0,0):('κ','Capital-Output Ratio', False),
        (0,1):('E','Efficiency of Labor', False),
        (1,0):('L','Log Labor Force', True),
        (1,1):('K','Log Capital Stock', True),
        (2,0):('Y','Log Output', True),
        (2,1):('y','Output-per-worker', False)
       }

num_rows, num_cols = 3,2
fig, axes = plt.subplots(num_rows, num_cols, figsize=(12, 12))
for i in range(num_rows):
    for j in range(num_cols):
        for m in m_base, m_alt:
            lb = f'{m.scenario}: initial κ = {m.initdata["κ"]:.2f}'
            axes[i,j].plot(m.gen_seq(T, var = figcontents[i,j][0], log = figcontents[i,j][2]),'o-', lw=2, alpha=0.5, label=lb)
            axes[i,j].set(title=figcontents[i,j][1])

#   global legend
axes[(0,0)].legend(loc='upper center', bbox_to_anchor=(1.1,1.3))
plt.suptitle('Malthusian Model: Simulation Run', size = 20)
plt.show()    


# In[3]:


# TESTING THE MALTHUSIAN CLASS
#
# Adding idea stock growth fast enough that
# (with γ=2) population will double every
# 700 years
# 
# Starting the economy at the h=0 Malthusian
# steady-state

import numpy as np
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')

import delong_classes

T = 1000

m_base = delong_classes.malthusian(h=.0005)
m_base.scenario = "base scenario"
m_alt = delong_classes.malthusian(h=.0005)
m_alt.scenario = "alt scenario"

figcontents = {
        (0,0):('κ','Capital-Output Ratio', False),
        (0,1):('E','Efficiency of Labor', False),
        (1,0):('L','Log Labor Force', True),
        (1,1):('K','Log Capital Stock', True),
        (2,0):('Y','Log Output', True),
        (2,1):('y','Output-per-worker', False)
       }

num_rows, num_cols = 3,2
fig, axes = plt.subplots(num_rows, num_cols, figsize=(12, 12))
for i in range(num_rows):
    for j in range(num_cols):
        for m in m_base, m_alt:
            lb = f'{m.scenario}: initial κ = {m.initdata["κ"]:.2f}'
            axes[i,j].plot(m.gen_seq(T, var = figcontents[i,j][0], log = figcontents[i,j][2]),'o-', lw=2, alpha=0.5, label=lb)
            axes[i,j].set(title=figcontents[i,j][1])

#   global legend
axes[(0,0)].legend(loc='upper center', bbox_to_anchor=(1.1,1.3))
plt.suptitle('Malthusian Model: Simulation Run', size = 20)
plt.show()    


# Notice that the steady state is not affected by where we start.

# ### What if a certain parameter changes value in the course of development?
# To start with, here we consider the time period 0AD - 500AD.
# 
# ### Assume useful ideas stopped developing since 250AD:

# In[4]:


import numpy as np
import matplotlib.pyplot as plt
get_ipython().run_line_magic('matplotlib', 'inline')

import delong_classes

m  = delong_classes.malthusian(h=0.0005, E=0.354, K=3.06)

# generate and store sequences before the change:
T1 = 250 # time before change
T2 = 250 # time after change

figcontents = {
        (0,0):('κ','Capital-Output Ratio', False),
        (0,1):('E','Efficiency of Labor', False),
        (1,0):('L','Log Labor Force', True),
        (1,1):('K','Log Capital Stock', True),
        (2,0):('Y','Log Output', True),
        (2,1):('y','Output-per-worker', False)
       }

num_rows, num_cols = 3,2
fig, axes = plt.subplots(num_rows, num_cols, figsize=(12, 12))
for i in range(num_rows):
    for j in range(num_cols):
        for scenario in {'base', 'with shock'}:
            seq = m.gen_seq(T1, var = figcontents[i,j][0], log = figcontents[i,j][2])
            lb = f'{scenario}'
            if scenario == 'with shock':
                m.h = 0
            seq = np.append(seq, m.gen_seq(T2, var = figcontents[i,j][0], log = figcontents[i,j][2], init = False))
            axes[i,j].plot(seq,'o-', lw=2, alpha=0.5, label=lb)
            axes[i,j].set(title=figcontents[i,j][1])

axes[(0,0)].legend(loc='upper center', bbox_to_anchor=(1.1,1.3))
plt.suptitle('Malthusian Model: Simulation Run with Idea Shock', size = 20)
plt.show()   


# In[5]:


# THE COMING OF AN IMPERIAL PEACE

# φ = 1.25, s = 0.25

import numpy as np
import matplotlib.pyplot as plt

import delong_classes

m  = delong_classes.malthusian(h=0.0005, E=0.354, K=3.06)

# generate and store sequences before the change:
T1 = 250 # time before change
T2 = 250 # time after change

figcontents = {
        (0,0):('κ','Capital-Output Ratio', False),
        (0,1):('E','Efficiency of Labor', False),
        (1,0):('L','Log Labor Force', True),
        (1,1):('K','Log Capital Stock', True),
        (2,0):('Y','Log Output', True),
        (2,1):('y','Output-per-worker', False)
       }

num_rows, num_cols = 3,2
fig, axes = plt.subplots(num_rows, num_cols, figsize=(12, 12))
for i in range(num_rows):
    for j in range(num_cols):
        for scenario in {'base', 'with shock'}:
            seq = m.gen_seq(T1, var = figcontents[i,j][0], log = figcontents[i,j][2])
            lb = f'{scenario}'
            if scenario == 'with shock':
                m.φ = 1.25
                m.s = 0.25
            seq = np.append(seq, m.gen_seq(T2, var = figcontents[i,j][0], log = figcontents[i,j][2], init = False))
            axes[i,j].plot(seq,'o-', lw=2, alpha=0.5, label=lb)
            axes[i,j].set(title=figcontents[i,j][1])

axes[(0,0)].legend(loc='upper center', bbox_to_anchor=(1.1,1.3))
plt.suptitle('Malthusian Model: Simulation Run with Idea Shock', size = 20)
plt.show()   


# In[6]:


# THE COLLAPSE OF AN EMPIRE

# φ = 1.25, s = 0.25

import numpy as np
import matplotlib.pyplot as plt

import delong_classes

m  = delong_classes.malthusian(h=0.0005, E=0.27, K=6.5, φ = 1.25, s = 0.25)

# generate and store sequences before the change:
T1 = 250 # time before change
T2 = 250 # time after change

figcontents = {
        (0,0):('κ','Capital-Output Ratio', False),
        (0,1):('E','Efficiency of Labor', False),
        (1,0):('L','Log Labor Force', True),
        (1,1):('K','Log Capital Stock', True),
        (2,0):('Y','Log Output', True),
        (2,1):('y','Output-per-worker', False)
       }

num_rows, num_cols = 3,2
fig, axes = plt.subplots(num_rows, num_cols, figsize=(12, 12))
for i in range(num_rows):
    for j in range(num_cols):
        for scenario in {'base', 'with shock'}:
            seq = m.gen_seq(T1, var = figcontents[i,j][0], log = figcontents[i,j][2])
            lb = f'{scenario}'
            if scenario == 'with shock':
                m.φ = 1.0
                m.s = 0.15
            seq = np.append(seq, m.gen_seq(T2, var = figcontents[i,j][0], log = figcontents[i,j][2], init = False))
            axes[i,j].plot(seq,'o-', lw=2, alpha=0.5, label=lb)
            axes[i,j].set(title=figcontents[i,j][1])

axes[(0,0)].legend(loc='upper center', bbox_to_anchor=(1.1,1.3))
plt.suptitle('Malthusian Model: Simulation Run with Idea Shock', size = 20)
plt.show()   


# ### Consider a plague happened on 250AD:
# The plague decreased the total population (and thus labor force) to 1/10. 

# In[7]:


import numpy as np
import matplotlib.pyplot as plt

import delong_classes

m  = delong_classes.malthusian(h=0.0005, E=0.354, K=3.06)

# generate and store sequences before the change:
T1 = 250 # time before change
T2 = 250 # time after change

figcontents = {
        (0,0):('κ','Capital-Output Ratio', False),
        (0,1):('E','Efficiency of Labor', False),
        (1,0):('L','Log Labor Force', True),
        (1,1):('K','Log Capital Stock', True),
        (2,0):('Y','Log Output', True),
        (2,1):('y','Output-per-worker', False)
       }

num_rows, num_cols = 3,2
fig, axes = plt.subplots(num_rows, num_cols, figsize=(12, 12))
for i in range(num_rows):
    for j in range(num_cols):
        for scenario in {'base', 'with plague'}:
            seq = m.gen_seq(T1, var = figcontents[i,j][0], log = figcontents[i,j][2])
            lb = f'{scenario}'
            if scenario == 'with plague':
                m.L = 2/3*m.L
                m.E = m.E*(2/3)**(-1/m.γ)
            seq = np.append(seq, m.gen_seq(T2, var = figcontents[i,j][0], log = figcontents[i,j][2], init = False))
            axes[i,j].plot(seq,'o-', lw=2, alpha=0.5, label=lb)
            axes[i,j].set(title=figcontents[i,j][1])

axes[(0,0)].legend(loc='upper center', bbox_to_anchor=(1.1,1.3))
plt.suptitle('Malthusian Model: Simulation Run with Plague', size = 20)
plt.show()   


# To summarize, the shock on idea development will change the steady state and the development path, while a plague would not change the steady state.

# ### <font color="000088"> Shocks to the System... </font>
# 
# * A plague, pushing $ L_t $ down substantially below $ L_t^{*mal} $...
# * A breakdown (or buildup) of law-and-order, raising incentives to save and invest, and so raising $ s $...
# * Invaders and other raiders, raising $ \delta $...
# * A speed-up or slowdown in innovation, raising or lowering $ h $...
# * Large-scale destruction of the societal web, lowering $ H $...
# * An increased rate of death from disease or violence, raising or lowering $ Y^{sub} $...
# 
# &nbsp;
# 
# ----

# ### <font color="000088"> Local Convergence... </font>
# 
# >(41) $ \ln(L_t^{*mal}) = \gamma \left[ \ln(H_t)  - \theta \ln(\gamma h +\delta) -ln\left(1 + \frac{\gamma h}{\beta} \right) \right] +  \gamma \left[ \theta \ln(s) \right] -  \gamma \left[ \ln(\phi) \right]  - \gamma \left[ \ln( y^{sub}) \right] $
# 
# And recall (31) for the full Malthusian equilibrium standard of living:
# 
# >(42) $ y^{*mal} = \phi y^{sub} \left( 1 + \frac{ \gamma h}{\beta}\right) $
# 
# Plus for the rate of population growth:
# 
# >(43) $ n^{*mal} = \gamma h $
# 
# Convergence:
# 
# >(44) $ \frac{d\kappa_t}{dt} = - (1-\alpha)((1-1/\gamma)n(y_t) + h + \delta)(\kappa_t - \kappa^*) $
# 
# >(45) $ \kappa^* = s/((1-1/\gamma)n(y_t) + h + \delta) $
# 
# >(46) $ n(y_t) = \beta  \left( \frac{y_t}{\phi y^{sub}}-1 \right) $
# 
# >(47) $ y_t = \left(\kappa_t\right)^\theta E_t   $
# 
# >(48) $ \frac{1}{E_t}\frac{dE_t}{dt} = h - n(y_t)/\gamma $
# 
# ----
# 
# >(49) $ \frac{d\kappa_t}{dt} = s - (1-\alpha)\left(\left(1-\frac{1}{\gamma}\right)\beta  \left( \frac{\left(\kappa_t\right)^\theta E_t }{\phi y^{sub}}-1 \right) + h + \delta \right)\kappa_t $
# 
# >(50) $ \frac{dE_t}{dt} = \left( h - \frac{\beta  \left( \frac{\left(\kappa_t\right)^\theta E_t }{\phi y^{sub}}-1 \right)}{\gamma} \right) E_t $
# 
# &nbsp;
# 
# ----

# \kappa_t & E_t ➔ y_t 
# 
# y_t ➔ n_t
# 
# n_t ➔ dE_t/dt
# 
# n_t ➔ d\kappa_t/dt
# 
# &nbsp;
# 
# ----

# ## <font color="880000"> Ancient Economies: A Malthusian Model </font>
# 
# <img src="https://tinyurl.com/20190119a-delong" width="300" style="float:right" />
# 
# ### <font color="000088">Catch Our Breath—Further Notes:</font>
# 
# <br clear="all" />
# 
# ----
# 
# * Weblog Support <https://github.com/braddelong/LS2019/blob/master/2019-08-17-Ancient_Economies.ipynb>
# * nbViewer <https://nbviewer.jupyter.org/github/braddelong/LS2019/blob/master/2019-08-17-Ancient_Economies.ipynb>
# 
# &nbsp;
# 
#