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

# ## <font color="880000"> Malthusian Model: Convergence </font>
# 
# Recall our Solow-Malthus model. The rate of growth of the labor force and population depends on luxury-adjusted income per worker $ y/\phi $ divided by subsistence $ y^{sub} $:
# 
# >(1) $ \frac{dL/dt}{L} = \frac{d\ln(L)}{dt} = n = \beta  \left( \frac{y}{\phi y^{sub}}-1 \right) $
# 
# Income per worker depends on the capital-output ratio $ \kappa $, the level of the useful ideas stock $ H $, and the amount of resource scarcity induced by the labor force $ L $:
# 
# >(2) $ \ln(y) = \theta\ln(\kappa) + \ln(H) - \ln(L)/\gamma $
# 
# >(3) $ \frac{d\kappa}{dt} = (1-\alpha)s - (h + (1 - 1/\gamma)n + ֿֿ\delta)\kappa $
# 
# With the parameters $ \alpha $ and $\theta $—the capital share of income and the capital-intensity elasticity of income—related by:
# 
# >(4) $ \theta = \frac{\alpha}{1-\alpha}$ and $ \alpha = \frac{\theta}{1+\theta} $
# 
# Substitute:
# 
# >(5) $ \frac{1}{L}\frac{dL}{dt} = \frac{d\ln(L)}{dt} = n = \beta  \left( \frac{\kappa^\theta H L^{-\gamma}}{\phi y^{sub}}-1 \right) $
# 
# >(6) $ \frac{d\kappa}{dt} = -(1-\alpha)(h + (1-1/\gamma)n +\delta)\kappa + (1-\alpha)s $
# 
# Define population-adjusted-for-ideas $ P $:
# 
# >(7) $ P = L H^{-\gamma} $
# 
# Then we have two state variables—capital-intensity $ \kappa $, the capital-output ratio, and population-adjusted-for-ideas $ P $. We have two dynamic equations: The rate of change of population-adjusted-for-ideas $ P $ is a function of the capital-output ratio and itself: 
# 
# >(8) $ \frac{1}{P}\frac{dP}{dt} = p = \beta  \left( \frac{\kappa^\theta }{P^{1/\gamma} \phi y^{sub}}-1 \right) -h\gamma $
# 
# And the rate of change of capital-intensity $ \kappa $ is a function of itself and of the rate of change of population-adjusted-for-ideas $ P $:
# 
# >(9) $ \frac{d\kappa}{dt} = (1-\alpha)s -(1-\alpha)(\gamma h +\delta)\kappa -(1-\alpha)\left(1-\frac{1}{\gamma}\right)\left(\frac{1}{P}\frac{dP}{dt}\right)\kappa  $

# The steady state is then:
# 
# >(10) $ P^{*mal} = \frac{L}{H^\gamma} =   \left[  \left( \frac{1}{\phi} \right) \left( \frac{\kappa^\theta }{  y^{sub}} \right) \left( \frac{1}{ 1 + h\gamma/\beta} \right) \right]^\gamma $
# 
# >(11) $ \kappa^{*mal} = \frac{s}{\gamma h + \delta} $

# 
# Define:
# 
# >(12) $ P = \Pi P^{*mal}     $
# 
# >(13) $ \kappa = k \kappa^{*mal} = k (s/(\delta + \gamma h)) $

# Then:
# 
# >(14) $ \kappa^{*mal} \frac{dk}{dt} =  \left[ (1-\alpha)s -(1-\alpha)((1-1/\gamma)p + \gamma h +\delta) k (s/(\delta + \gamma h)) \right]  $
# 
# >(15) $ \frac{1}{\Pi}\frac{d\Pi}{dt} = p =  \beta \Pi^{(-1/\gamma)} k^\theta  (1 + \gamma h/\beta) - \beta - h \gamma $
# 
# >(16) $ \Pi^{(-1/\gamma)} = 1 - \frac{\Delta \Pi}{\gamma} $
# 
# >(17) $ k^\theta = 1 + \theta k $
# 
# >(18) $ \frac{1}{\Pi}\frac{d\Pi}{dt} = p =  \beta  (1 + \gamma h/\beta)\left[ 1 - \frac{\Delta \Pi}{\gamma} \right]\left[ 1 + \theta \Delta k \right] - \beta - h \gamma $
# 
# >(19) $ \frac{dk}{dt} =  (1-\alpha)( \gamma h +\delta) - (1-\alpha)(1-1/\gamma)(1 + \Delta k) \frac{d\Pi}{dt} - (1-\alpha)(\delta + h \gamma)(1 + \Delta k)  $
# 
# >(20) $ \frac{dk}{dt} =  - (1-\alpha)(\delta + h \gamma)\Delta k - (1-\alpha)(1-1/\gamma)\frac{d\Pi}{dt} - (1-\alpha)(1-1/\gamma)\Delta k \frac{d\Pi}{dt}   $
# 
# >(21) $ \frac{d\Pi}{dt} = -((\beta + \gamma h)/\gamma)\Delta \Pi + (\beta + \gamma h)\theta \Delta k $
# 
# 

# (20) and (21) are a linearized two-equation two state-variable system for $ k = \frac{\kappa}{\kappa^*} $ and $ \Pi = \frac{P}{P^{*mal}} $, with $ \kappa = K/Y $ the capital intensity of the economy, and $ P = L H^{-\gamma} $ the population/labor force adjusted for the ideas stock.
# 
# &nbsp;

# **Memo**: $ \alpha = 1/2, \gamma = 2, \theta = 1, \delta = 0.05, \beta = 0.025, h=0.0005, h \gamma = 0.001, \delta + h \gamma = .051 $
# 
# >(20') $ \frac{dk}{dt} =  - 0.0255\Delta k - 0.25\frac{d\Pi}{dt} - 0.25\Delta k \frac{d\Pi}{dt}   $
# 
# >(21') $ \frac{d\Pi}{dt} = - 0.013\Delta \Pi + 0.026 \Delta k  $

# >(22') $ \frac{dk}{dt} =  - 0.032\Delta k + 0.00325 \Delta \Pi   $
# 
# &nbsp;
# 
# ----

# In[22]:


import numpy as np
import pandas as pd

P_vector = []
kappa_vector = []

P_vector = P_vector + [1.941371*.75]
kappa_vector = kappa_vector + [1.980198]



# In[23]:


for i in range(0,100):
    P = (1-0.025-0.0005)*P_vector[i] + 0.025*kappa_vector[i]
    kappa = (1-0.02525)*kappa_vector[i]+0.05
    P_vector = P_vector + [P]
    kappa_vector = kappa_vector + [kappa]
    
malthus_converge_df = pd.DataFrame()
malthus_converge_df['P'] = P_vector
malthus_converge_df['kappa'] = kappa_vector

print(malthus_converge_df)


# ## <font color="880000"> Malthusian: Convergence </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/long-form-drafts/blob/master/malthusian_convergence.ipynb>
# * nbViewer <https://nbviewer.jupyter.org/github/braddelong/long-form-drafts/blob/master/malthusian_convergence.ipynb>
# 
# &nbsp;