# The basic solow model and extension with human capital

Imports for use in project:


```python
import numpy as np
from scipy import optimize
import sympy as sm
from types import SimpleNamespace
import matplotlib.pyplot as plt
import ipywidgets as widgets
from types import SimpleNamespace
import warnings

# Autoreload modules when code is run
%load_ext autoreload
%autoreload 2

# Import our python file
import modelproject as mp
```

# Basic Solow model

**Write out the model in equations here.** 

Make sure you explain well the purpose of the model and comment so that other students who may not have seen it before can follow.  

## Analytical solution

If your model allows for an analytical solution, you should provide here.

You may use Sympy for this. Then you can characterize the solution as a function of a parameter of the model.

To characterize the solution, first derive a steady state equation as a function of a parameter using Sympy.solve and then turn it into a python function by Sympy.lambdify. See the lecture notes for details. 

## Numerical solution

To find a numerical solution for the steady state of capital we use the Solow-equation:
$$ k_{t+1} - k_t = \frac{1}{1+n} \left[sBk_t^\alpha - (n+\delta) k_t \right] $$
And the fact that in steady state we have: $$ k_{t+1}=k_t=k^* $$

This gives us that in steady state it must follow that:
$$ 0 = \frac{1}{1+n} \left[sBk^{*\alpha} - (n+\delta) k^* \right]$$ 
$$\Leftrightarrow 0 = sBk^{*\alpha} - (n+\delta) k^* \tag{1} $$ 

As plausible values for our parameters we look at table A pages 349-351 in the book: **Introducing Advanced Macroeconomics - Growth and business cycles** by *Peter Birch Sørensen and Hans Jørgen Whitta-Jacobsen*, second edition. We choose values for Denmark.  


```python
# Plausible parameter values as found in table A
par = SimpleNamespace()
par.alpha = 1/3
par.s = 0.333
par.B = 1 
par.n = 0.008
par.delta = 0.01
```

We construct an algorithm to find the steady state of capital per worker and use it with the plausible values.


```python
# Algorithm to find the root of a function on the interval [a,b] with the bisect method
def bisection(f,a,b,tol=1e-8):
    """  
    Args:
    f      (function): Function you wish to find root of
    a      (float)   : Lower bound
    b      (float)   : Upper bound
    tol    (float)   : Tolerance on solution
        
    Returns:
    The x value that solves equation f(x) = 0 for a <= x <= b.
    
    """
    
    # Test inputs
    if f(a)*f(b) >= 0: 
        print("Bisection method fails.") # If f(a) and f(b) both have same sign it will not have a root on given interval
        return None
    
    # Step 1: initialize
    a_n = a
    b_n = b
    
    # Step 2-4:
    while True:
        
        # Step 2: midpoint and associated value
        m_n = (a_n+b_n)/2 # The midpoint of interval is found
        f_m_n = f(m_n) # The function value of the midpoint (m) is found
        
        # Step 3: Check if solution is found and if not determine sub-interval for evaluation in next iteration
        if abs(f_m_n) < tol: # If the function value of the midpoint is close enough to 0 we have found our solution
            return m_n
        elif f(a_n)*f_m_n < 0: # If f(a) and f(m) have different signs then m is the new upper bound
            a_n = a_n
            b_n = m_n
        elif f(b_n)*f_m_n < 0: # If f(b) and f(m) have different signs then m is the new lower bound
            a_n = m_n
            b_n = b_n
        else:
            print("Bisection method fails.")
            return None
    # Step 4: results 
    return (a_n + b_n)/2

# Use the algorithm to find the steady state of capital (k*)
f = lambda k: par.s*par.B*k**par.alpha - (par.n + par.delta)*k # Equation (1) is the one we wish to solve
print(f'The steady state of capital found with algorithm is : k* = {bisection(f,0.1,100,1e-8):.3f}')
```

    The steady state of capital found with algorithm is : k* = 79.572


We can also solve the problem by using a built in optimizer from the scipy package that uses the bisection method.


```python
print(f'The steady state of capital found by using built in optimizer is :k* = {mp.ss_bas(par.n, par.s, par.B, par.alpha, par.delta):.3f}')
```

    The steady state of capital found by using built in optimizer is :k* = 79.572


Unsurprisingly the results are the same.

From equation (1) we can see that in order for the model to converge to a steady state of capital then the stability condition of $$ n + \delta > 0 $$ must be fulfilled. If this condition is met then the model will converge for any starting values.

## Further analysis

We plot the basic solow model with interactive widgets to examine how changes in parameters affect the model. We also find the steady states of: capital per worker (k*), output per worker (y*) and consumption per worker (c*). 


```python
widgets.interact(mp.ss_bas_plot, 
                alpha = widgets.FloatSlider(description = r'$\alpha$', min = 0.01, max = 0.6, step = 0.005, value = 0.3, readout_format='.3f'),
                delta = widgets.FloatSlider(description = r'$\delta$', min = 0, max = 0.1, step = 0.005, value = 0.01, readout_format='.3f'),
                s = widgets.FloatSlider(description = '$s$', min = 0.01, max = 0.7, step = 0.005, value = 0.4, readout_format='.3f'),
                n = widgets.FloatSlider(description ='$n$', min = -0.05, max = 0.05, step = 0.005, value = 0.01, readout_format='.3f'),
                B = widgets.fixed(1), #B is assumed to be the same across countries
                k_max = widgets.IntSlider(description='k_max', min = 1, max = 1000, step = 10, value = 100)) # Changes the size of x-axis to zoom
```


    interactive(children=(FloatSlider(value=0.01, description='$n$', max=0.05, min=-0.05, readout_format='.3f', st…





    <function modelproject.ss_bas_plot(n, s, B, alpha, delta, k_max)>



From figure 1 and the interactive widgets we can see the ceteris paribus effects of changes in the different parameters. We find that: n is positively correlated with k*, y* and c*, s is positively correlated with k* and y* but negatively correlated with c*, $\alpha$ is positively correlated with k*, y* and c* while $\delta$ is negatively correlated with k*, y* and c*.  
  
In order to ensure convergance we must have that $ n + \delta > 0 $

# Solow model with human capital

**Write out the model in equations here.** 

Make sure you explain well the purpose of the model and comment so that other students who may not have seen it before can follow.  


```python
# Plausible parameter values as found in table A
par.s_H = 0.1
par.s_K = 0.2
par.g = 0.02
par.n = 0.01
par.alpha = 1/3
par.phi = 1/3
par.delta = 0.01
```

## Analytical solution

If your model allows for an analytical solution, you should provide here.

You may use Sympy for this. Then you can characterize the solution as a function of a parameter of the model.

To characterize the solution, first derive a steady state equation as a function of a parameter using Sympy.solve and then turn it into a python function by Sympy.lambdify. See the lecture notes for details. 

## Numerical solution

The equations for the steady states of physical capital per capita and human capital per capita can be written as

\\[ 0 = \left(\dfrac{s_k^{1-\varphi}s_h^{\varphi}}{n+g+\delta+ng}\right)^{\frac{1}{1-\alpha-\varphi}} - \tilde{k}^{\ast} \\]

\\[ 0 = \left(\dfrac{s_k^{\alpha}s_h^{1-\alpha}}{n+g+\delta+ng}\right)^{\frac{1}{1-\alpha-\varphi}} - \tilde{h}^{\ast} \\]


```python
warnings.filterwarnings("ignore", category=RuntimeWarning) # Some values of parameters give a runtimewarning, this is ignored 

# We define a function containing our h- and k-functions, and a vector x = [h,k]
obj = lambda x: [mp.h_func(x[1],par.s_H,par.s_K,par.g,par.n,par.alpha,par.phi,par.delta,x[0]),mp.k_func(x[0],par.s_H,par.s_K,par.g,par.n,par.alpha,par.phi,par.delta,x[1])]

# We solve the vector functions, have nonlinear system of equation and therefore use broyden1 method
sol = optimize.root(obj,[1,1],method = 'broyden1')

# The numerical solution is found
num_sol1 = sol.x

print(f'Numerical solution is: k* = {num_sol1[1]:.3f}, h* = {num_sol1[0]:.3f}')
```

    Numerical solution is: k* = 61.572, h* = 30.786



```python
par.incr = 0.01 # Define a variable called incr to examine the effects of increase in human capital

# Solve the functions, first with original value of human capital and then with increased human capital
k_vec, h_vec_DeltaK0, h_vec_DeltaH0 = mp.solve_ss(par.s_H,par.s_K,par.g,par.n,par.alpha,par.phi,par.delta)
k1_vec, h1_vec_DeltaK0, h1_vec_DeltaH0 = mp.solve_ss(par.s_H+par.incr,par.s_K,par.g,par.n,par.alpha,par.phi,par.delta)

# Find steady state values for new s_H

# # We define a function containing our h- and k-functions, and a vector x = [h,k] for increased human capital
obj1 = lambda x: [mp.h_func(x[1],par.s_H+par.incr,par.s_K,par.g,par.n,par.alpha,par.phi,par.delta,x[0]),mp.k_func(x[0],par.s_H+par.incr,par.s_K,par.g,par.n,par.alpha,par.phi,par.delta,x[1])]

# We solve the vector functions, have nonlinear system of equation and therefore use broyden1 method
sol1 = optimize.root(obj1,[1,1],method = 'broyden1')

# The numerical solution is found for increased human capital
num_sol2 = sol1.x

# Create the plot
fig2 = plt.figure(figsize=(7,7))
ax = fig2.add_subplot(1,1,1)
ax.plot(k_vec,h_vec_DeltaK0, label=r'$\Delta \tilde{k}=0$', c='lime')
ax.plot(k_vec,h_vec_DeltaH0, label=r'$\Delta \tilde{h}=0$', c='aqua')
ax.plot(k1_vec,h1_vec_DeltaH0, label=r'$\Delta \tilde{h_1}=0$', c='teal')
ax.set_xlabel(r'$\tilde{k}$')
ax.set_ylabel(r'$\tilde{h}$')
ax.legend() 

# We mark the steady state
plt.scatter(sol.x[1],sol.x[0],color='black',s=80,zorder=3)
plt.scatter(sol1.x[1],sol1.x[0],color='black',s=80,zorder=3)

# Lines are drawn to mark ss-value on the axes
plt.axvline(sol.x[1],ymax=sol.x[0]/60,color='gray',linestyle='--') 
plt.axhline(sol.x[0],xmax=sol.x[1]/80,color='gray',linestyle='--') 
plt.axvline(sol1.x[1],ymax=sol1.x[0]/60,color='black',linestyle='--') 
plt.axhline(sol1.x[0],xmax=sol1.x[1]/80,color='black',linestyle='--') 

# Text is added to the plot
ax.text(0.15,sol.x[0]+1 , r'$\tilde{h}^*$', fontsize=12)
ax.text(sol.x[1]+1,0.15 , r'$\tilde{k}^*$', fontsize=12)
ax.text(0.15,sol1.x[0]+1 , r'$\tilde{h_1}^*$', fontsize=12)
ax.text(sol1.x[1]+1,0.15 , r'$\tilde{k_1}^*$', fontsize=12)

#The axis limits are chosen
ax.set(xlim=(0, 80), ylim=(0, 60))

ax.set_title('Figure 2. Phase diagram');

print(f'Numerical solution for original s_H is: k* = {num_sol1[1]:.3f}, h* = {num_sol1[0]:.3f}')
print(f'Numerical solution for increased s_H is: k* = {num_sol2[1]:.3f}, h* = {num_sol2[0]:.3f}')
```

Figure 2 shows the phase diagram and the steady state is found where the nullclines intersect. We can see that an increase in $s_H$ increases the steady states of both $\tilde{h}$ and $\tilde{k}$.

# Conclusion

Add concise conclusion. 

ef01ea8660a312a14cff06af889f70554d3bbc54
ipynb
Jupyter Notebook
modelproject/Solow model.ipynb
NumEconCopenhagen/projects-2022-msm
619192f54879f2494b3eab121f183a0869a1064d
MIT
modelproject/Solow model.ipynb
NumEconCopenhagen/projects-2022-msm
619192f54879f2494b3eab121f183a0869a1064d
MIT
2022-03-28T15:23:11.000Z
2022-03-31T10:45:29.000Z
modelproject/Solow model.ipynb
NumEconCopenhagen/projects-2022-msm
619192f54879f2494b3eab121f183a0869a1064d
MIT
Qwen/Qwen-72B
1. YES
2. YES