# Model Project 2021: Solow model with economical damage from climate change

The aim of the project is to solve the basic Solow model and then extend the model to account for economical damage from climate change. 

**The project is organized as follows:**

1. Model project description
2. Model description
3. Classic Solow model
4. Extended Solow model (with economic damage from climate change)
5. Conclusion

**Group:** SSTL

**Members:** Louise Otte Arildsen (srb330), Signe Kolind (hmk792), Stine Fürst (qvj635) & Thomas Nielsen (qzw392)

**Imports and set magics:**


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


```python
import modelproject as mp
```

# Model description

$$
\begin{align}
Y_{t}&=[1-D] K_{t}^{\alpha}\left(A_{t} L_{t}\right)^{1-\alpha}, 0<\alpha<1 \\
D&=1-\frac{1}{1+\sigma(\Delta T)^{2}}, \sigma>0 \\
K_{t+1}&=s Y_{t}+(1-\delta) K_{t}, 0<s<1,0<\delta<1, K_{0}>0  \\
A_{t+1}&=(1+g) A_{t}, A_{0}>0  \\
L_{t+1}&=(1+n) L_{t}, L_{0}>0
\end{align}
$$

$$
\tilde{y}_{t} \equiv \frac{Y_{t}}{A_{t} L_{t}}=\frac{y_{t}}{A_{t}} ; \quad \tilde{k}_{t} \equiv \frac{K_{t}}{A_{t} L_{t}}=\frac{k_{t}}{A_{t}}
$$

Equation (1) is a Cobb-Douglas production function, which describes production ($Y_t$) as a function of capital ($K_t$), population size ($L_t$), and the level of technology ($A_t$). The model assumes that climate changes reduce production with a share D, meaning that a share of the production ($DY_t$) is damaged and can not be used. The share (1-D) of the production can be used for consumption and investment. Equation 2 describes the damage caused by climate change, where \DeltaT measures the temperature change. Equation (3) describes the development of capital, where s is the saving rate and delta is the depreciation rate. Equation (4) and (5) describes how the level of technology ($A_t$) and the population size develops over time.

# Classic Solow model (D=0)

When assuming $D=0$ the model reduces to the classic Solow model, which we will solve in this section.

## Analytical solution

The transition equation of the model for capital per effective worker is given by: 

$$
\tilde{k}_{t+1}=\frac{1}{(1+n)(1+g)}\left(s \tilde{k}_{t}^{\alpha}+(1-\delta) \tilde{k}_{t}\right)
$$

Note that convergence requires that $n+g+\delta+ng >0$.

The intuition behind the convergence result is decreasing marginal product on capital and constant depletion/attrition. In steady-state this implies that the savings are just as large as they cover for depreciation and thinning of resources.

First we define all the **symbols** using sympy:


```python
y = sm.symbols('y')
k = sm.symbols('k')
K = sm.symbols('K')
alpha = sm.symbols('alpha')
sigma = sm.symbols('sigma')
deltaT = sm.symbols('delta T')
s = sm.symbols('s')
delta = sm.symbols('delta')
g = sm.symbols('g')
n = sm.symbols('n')
```

Then we define the **steady state equation**


```python
ss = sm.Eq(k,(s*k**alpha+(1-delta)*k)/((1+n)*(1+g)))
```

and **solve** it


```python
kss = sm.solve(ss,k)[0]
kss
```




$\displaystyle \left(\frac{\delta + g n + g + n}{s}\right)^{\frac{1}{\alpha - 1}}$



For later use, we turn the solution into a **Python funciton**


```python
#transforming the function to a python function
ss_func = sm.lambdify((s,g,n,alpha,delta),kss)
```


```python
#Defining parameters
s = 0.3
g = 0.02
n = 0.01
alpha = 1/3
delta = 0.05
```


```python
print(f'analytical solution is: {ss_func(s,g,n,alpha,delta):.3f}')
```

    analytical solution is: 7.235


In steady state the capital per effective worker equals 7.235.

## Numerical solution

We can re-write the equation for the steady state capital per capita as

$$
\tilde{k}_{t+1}-\frac{1}{(1+n)(1+g)}\left(s \tilde{k}_{t}^{\alpha}+(1-\delta) \tilde{k}_{t}\right)=0
$$

whereby it clearly becomes a **root-finding problem**. Such a problem can be solved by various methods. 

In the section below, the problem will be solved from the **brent-** and **bisection method** in order to compare the robustness of our results.


### Bisection Method


```python
def solve_for_ss(s,g,n,alpha,delta):
    """ solve for the steady state level of capital

    Args:
        s (float): saving rate
        g (float): technological growth rate
        n (float): population growth rate
        alpha (float): cobb-douglas parameter
        delta (float): capital depreciation rate 

    Returns:
        result (RootResults): the solution represented as a RootResults object

    """ 
    
    # a. define objective function
    obj_kss = lambda kss: kss - (s*kss**alpha + (1-delta)*kss)/((1+g)*(1+n))

    #. b. call root finder
    result = optimize.root_scalar(obj_kss,bracket=[0.1,100],method='bisect')
    
    return result


#Solve numerically for the steady state:
bisect_solution = solve_for_ss(s,g,n,alpha,delta)
print(f' The numerical solution, calculated from the bisect method is: {bisect_solution.root:}')
```

     The numerical solution, calculated from the bisect method is: 7.234696712204852



```python
def solve_for_ss(s,g,n,alpha,delta,method='bisect'):
    """ solve for the steady state level of capital

    Args:
        s (float): saving rate
        g (float): technological growth rate
        n (float): population growth rate
        alpha (float): cobb-douglas parameter
        delta (float): capital depreciation rate 

    Returns:
        result (RootResults): the solution represented as a RootResults object

    """ 
    
    # a. define objective function
    obj_kss = lambda kss: kss - (s*kss**alpha + (1-delta)*kss)/((1+g)*(1+n))

    #. b. call root finder
    result = optimize.root_scalar(obj_kss,bracket=[0.1,100],method=method)
    
    return result


#Solve numerically for the steady state:
bisect_solution = solve_for_ss(s,g,n,alpha,delta)
print(f' The numerical solution, calculated from the bisect method is: {bisect_solution.root:}')
```

     The numerical solution, calculated from the bisect method is: 7.234696712204852



```python
#Solve numerically for the steady state:
brent_solution = solve_for_ss(s,g,n,alpha,delta,method='brentq')
print(f'The numerical solution, calculated from the brent method, is: {brent_solution.root:}')
```

    The numerical solution, calculated from the brent method, is: 7.234696712205987


### Brent Method


```python
def solve_for_brent(s,g,n,alpha,delta):
    """ solve for the steady state level of capital

    Args:
        s (float): saving rate
        g (float): technological growth rate
        n (float): population growth rate
        alpha (float): cobb-douglas parameter
        delta (float): capital depreciation rate 

    Returns:
        result (RootResults): the solution represented as a RootResults object

    """ 
    
    # a. define objective function
    obj_kss = lambda kss: kss - (s*kss**alpha + (1-delta)*kss)/((1+g)*(1+n))

    #. b. call root finder
    result = optimize.root_scalar(obj_kss,bracket=[0.1,100],method='brentq')
    
    return result

#Solve numerically for the steady state:
brent_solution = solve_for_brent(s,g,n,alpha,delta)
print(f'The numerical solution, calculated from the brent method, is: {brent_solution.root:}')
```

    The numerical solution, calculated from the brent method, is: 7.234696712205987




When comparing our numerical and analytical results, it is clear that our results are quite robust, only differing by few decimals:


```python
print(f'Analytical solution : {ss_func(s,g,n,alpha,delta):}')
print(f'Bisect solution     : {bisect_solution.root:}')
print(f'Brent solution      : {brent_solution.root:}')

```

    Analytical solution : 7.234696712206231
    Bisect solution     : 7.234696712204852
    Brent solution      : 7.234696712205987


Because our results seem robust, we choose to rely on a single optimization method, when climate change is introduced from section 4 in the extended model. 

## GDP per effective worker:

From equation 1 we divide both sides with $A_{t}L_{t}$ and obtain the following expression for GDP per effective worker:

$$
\begin{equation*}
\tilde{y}_{t}=\tilde{k}_{t}^\alpha
\end{equation*}
$$

We insert the expression for $\tilde{k}^*=\left(\frac{\delta + gn + g+n}{s}\right)^{\frac{1}{\alpha - 1}}$ found in section 3.1 in $\tilde{k}_{t}$'s  place. 
We then get this expression for GDP per effective worker: 

$$
\tilde{y}^{*}=\left(\frac{s}{n+g+\delta+n g}\right)^{\frac{\alpha}{1-\alpha}}
$$

We calculate $\tilde{y}^{*}$ given our parameter values:


```python
gdp_pr_worker=(s/(n+g+delta+n*g))**(alpha/(1-alpha))
print(f'The GDP per effective worker is: {gdp_pr_worker:.3f}')
```

    The GDP per effective worker is: 1.934


# The Solow model with economical damage from climate change (D>0)

In this section we consider an extension of the classic Solow model, which includes economical damage from climate change. D represents the reduction in production due to climate change and we now assume that D>0.  

## Analytical solution

The transition equation of the model for capital per effective worker is given by: 

$$
\tilde{k}_{t+1}=\frac{1}{(1+n)(1+g)}\left(s[1-D] \tilde{k}_{t}^{\alpha}+(1-\delta) \tilde{k}_{t}\right)
$$

Then we define the **steady state equation**


```python
# introducing climate effects
y = sm.symbols('y')
k = sm.symbols('k')
K = sm.symbols('K')
alpha = sm.symbols('alpha')
sigma = sm.symbols('sigma')
deltaT = sm.symbols('delta T')
s = sm.symbols('s')
delta = sm.symbols('delta')
g = sm.symbols('g')
n = sm.symbols('n')
D = sm.symbols('D')

# climate steady state
css = sm.Eq(k,(((s*(1-D)*(k**alpha))+((1-delta)*k)))/((1+n)*(1+g)))
```

and **solve** it:


```python
kcss = sm.solve(css,k)[0]
kcss
```




$\displaystyle \left(\frac{- \delta - g n - g - n}{s \left(D - 1\right)}\right)^{\frac{1}{\alpha - 1}}$




```python
#transforming the function to a python function
css_func = sm.lambdify((s,g,n,alpha,delta,D),kcss)
```


```python
#Defining parameters
s = 0.3
g = 0.02
n = 0.01
alpha = 1/3
delta = 0.05
D=0.175
print(f'analytical solution is: {css_func(s,g,n,alpha,delta,D):.3f}')
```

    analytical solution is: 5.421


In steady state the capital per effective worker equals 5.421.

## Numerical solution


```python
def solve_for_css(s,g,n,alpha,delta,D):
    """ solve for the steady state level of capital

    Args:
        s (float): saving rate
        g (float): technological growth rate
        n (float): population growth rate
        alpha (float): cobb-douglas parameter
        delta (float): capital depreciation rate 

    Returns:
        result (RootResults): the solution represented as a RootResults object

    """ 
    
    # a. define objective function
    obj_kcss = lambda kcss: kcss - (s*(1-D)*kcss**alpha + (1-delta)*kcss)/((1+g)*(1+n))

    #. b. call root finder
    result = optimize.root_scalar(obj_kcss,bracket=[0.1,100],method='bisect')
    
    return result


```


```python
solution_2 = solve_for_css(s,g,n,alpha,delta,D)

print(f'analytical solution is: {css_func(s,g,n,alpha,delta,D):.3f}')
print(f' numerical solution is: {solution_2.root:.3f}')
```

    analytical solution is: 5.421
     numerical solution is: 5.421


## GDP per effective worker with climate change (D = 0.175)

There are different scenarios for the future temperature changes. A scenario with "little CO2 emission" will imply that the global average temperature increases by approx. 1 degree over the next 100 years. The worst case scenario, with "a lot of CO2 emissions", will imply that the global temperature rise at about 4 degrees over the next 100 years. A meta-study by Nordhaus og Moffat (2017) shows that the most pessimistic estimates in the literature suggests that temperature rises of 4 degress are associated with approx. 17.5% loss of production i.e. D = 0.175.

We want to compare the GDP per effective worker in the classic Solow model (D=0) with the GDP per effective in the extended Solow model when the worst possible climate scenario unfolds (D=0.175). We start by calculating the GDP per effecitive worker with D=0.175:

From equation 1 we divide both sides with $A_{t}L_{t}$ and obtain the following expression for GDP per effective worker:

$$
\begin{equation*}
\tilde{y}_{t}=(1-D)\tilde{k}_{t}^\alpha
\end{equation*}
$$

We insert the expression for $\tilde{k}^*=\left(\frac{\delta + gn + g+n}{s(1-D)}\right)^{\frac{1}{\alpha - 1}}$ found in section 4.1 in $\tilde{k}_{t}$'s  place. 
We then get this expression for GDP per effective worker: 

$$
\tilde{y}^{*}=(1-D)^{\frac{1}{1-\alpha}} \left(\frac{s}{n+g+\delta+n g}\right)^{\frac{\alpha}{1-\alpha}}
$$

We calculate $\tilde{y}^{*}$ given our parameter values:


```python
D=0.175
gdp_pr_worker_climate=(1-D)**(1/(1-alpha))*(s/(n+g+delta+n*g))**(alpha/(1-alpha))
print(f'The GDP per effective worker when accounting for climate change is: {gdp_pr_worker_climate:.3f}')
```

    The GDP per effective worker when accounting for climate change is: 1.449


### Compare GDP per effective worker with the model without climate change


```python
comparison= 1-gdp_pr_worker_climate/gdp_pr_worker
print(f'The reduction in GDP per effective worker when accounting for climate change: {comparison:.3f}')
```

    The reduction in GDP per effective worker when accounting for climate change: 0.251


This imply that GDP per efficient worker decreases by approx. 25 pct. This effect is greater than the initial decline of 17.5 pct. The reason for this difference is due to negative effects via capital accumulation.


```python
def Solaris(D,alpha,s,n,g,delta):
    return 1-((1-D)**(1/(1-alpha))*(s/(n+g+delta+n*g))**(alpha/(1-alpha)))/(s/(n+g+delta+n*g))**(alpha/(1-alpha))
```

We create an interactive widget where climate changes reduces the production with a share D and thus also reduces the GDP per effective worker:


```python
widgets.interact(Solaris, 
                 alpha= widgets.fixed(1/3),
                 s=widgets.fixed(0.3),
                 g=widgets.fixed(0.02),
                 n=widgets.fixed(0.01),
                 delta=widgets.fixed(0.05),
                 D= widgets.FloatSlider(description ='$D$',        
                                        min = 0.0 , max = 1 , step = 0.005, value = 0.17));
```


    interactive(children=(FloatSlider(value=0.17, description='$D$', max=1.0, step=0.005), Output()), _dom_classes…


We see that for any given level of D the reduction in GDP per capita is greater than the reduction in production due to the negative affect froms capital accumulation.

## Visulization

The Solow model attemps to explain how long-run economic growth is determined from capital accumulation, population(labor) growth and increase in productivity. The key takeway, is that the economy will grow towards a steady state, as long as per capita investments sf(k) is larger than what is needed to maintain the current level of capital. In order to maintain the population at the given level of k, both labor and capital must be “employed” and replaced by new units, as of regards to depreciation and population growth k(n+d+g) and in our case *climate change, D*. As shown in the transition equation, D, will essentially decrease the ammount of savings in the economy, which implies that the economy will converge towards a lower steady state, since less capital can be accumulated. 

The process is illustrated in the below visualization of the Solow model. 
Note that the default setting of widget paratemeter "D" is value=0. 

In order to view the effects of climate changes and compare with our prior steady state result, please change the value of D to value = 0.175):


```python
def simulate_basic(k,alpha,delta,s,n,g,K,T,D):
    """
    Args:
    k      (int)  : initial value for physical capital per worker
    alpha  (float): return on physical capital
    delta  (float): depreciation rate on capital
    s      (float): savings/investments
    n      (float): increase in labour
    g             : growth
    K      (int)  : totalfactorproductivity
    T      (int)  : periods to simulate through
    D             : climate change
    
    Returns:
    A plot showing how capital per worker accumulates over time"""
```


```python
widgets.interact(mp.simulate_basic, 
                k     = widgets.fixed(0), 
                alpha = widgets.FloatSlider(description = '$\u03B1$' , min = 0 ,    max = 0.9 , step = 0.05 , value = 1/3),
                delta = widgets.FloatSlider(description = '$\u03B4$' , min = 0 ,    max = 0.5 , step = 0.01 , value = 0.05),
                s     = widgets.FloatSlider(description = 's' ,        min = 0.01 , max = 0.3 , step = 0.005 , value = 0.3),
                n     = widgets.FloatSlider(description ='$n$',        min = 0.01 , max = 0.1 , step = 0.005, value = 0.01),
                g     = widgets.FloatSlider(description ='$g$',        min = 0.01 , max = 0.1 , step = 0.005, value = 0.02),
                D     = widgets.FloatSlider(description ='$D$',        min = 0.0 , max = 1 , step = 0.005, value = 0),
                K     = widgets.fixed(1),
                T     = widgets.IntSlider(description='$T$' ,          min = 1,     max = 1000, step = 10,    value = 10))
```


    interactive(children=(FloatSlider(value=0.3333333333333333, description='$α$', max=0.9, step=0.05), FloatSlide…





    <function modelproject.simulate_basic(k, alpha, delta, s, n, g, K, T, D)>



## Simulation of the classic and extended Solow model 

We want to simulate the classic and the extended Solow model, to show how the two models converges to their steady state, $\tilde{k}^*$.

We choose the initial capital and create a function that simulates the Solow model for $T=100$ periods.


```python
k0 =1
T = 100

def solow_equation(k,alpha,delta,s,n,g, D):
    """ calculate capital in the next-period
        
    Args:
    
        k (float): capital in this period
        alpha (float): cobb-douglas parameter
        delta (float): depreciation rate
        s (float): saving rate
        n : population growth
        g : technology growth
    
    Returns:
    
        k_plus (float): capital in next period
        
    """
    
    saving = s*k**alpha
    depreciation = delta*k
    k_plus = (k + saving*(1-D) - depreciation)/((1+n)*(1+g))    
    
    return k_plus
 

def simulate_solow_model(k0,alpha,delta,s,n,g,T,D):
    """ simulate the solow model
        
    Args:
    
        k0 (float): initial level of kapital
        alpha (float): cobb-douglas parameter
        delta (float): depreciation rate
        s (float): saving rate
        T (int): number of periods to simulate
    
    Returns:
    
        k_path (list): capital path (length T)
        
    """
    
    # a. initialize
    k_path = [k0]  
    
    # b. simulate forward
    for t in range(1,T):
        
        # i. apply solow equation
        k_plus = solow_equation(k_path[t-1],alpha,delta,s,n,g,D)    
        
        # ii. save result
        k_path.append(k_plus)
        
    return k_path
```

This function represents the steady state value:


```python
def find_ssk_k(delta,g,n,alpha,s,D):
    return ((delta+g*n+g+n)/(s*(1-D)))**(1/(alpha-1))
```

**The classic Solow model (D=0)**

We create the Solow path and the steady state value:


```python
k_path = simulate_solow_model(k0,alpha,delta,s,n,g,T,D=0)
k_star = find_ssk_k(delta,g,n,alpha,s,D=0)
k_stars =[k_star for i in range(T)]
```


```python
fig = plt.figure()
ax = fig.add_subplot(1,1,1)

ax.plot(k_path[:T], lw=2, label='Solow path')
ax.plot(k_stars,lw=2, label="Steady state")
ax.set_title('Simulated convergence of capital to steady state')
ax.set_xlabel('t')
ax.set_ylabel('Capital per capita')
ax.set_xlim([0,100])
ax.set_ylim([0,10])
ax.legend(loc="lower right")
ax.grid(True);
```

**The extended Solow model (D=0.175)**

We create the Solow path and the steady state value:


```python
k_path_ext = simulate_solow_model(k0,alpha,delta,s,n,g,T,D=0.175)
k_star_ext = find_ssk_k(delta,g,n,alpha,s,D=0.175)
k_stars_ext =[k_star_ext for i in range(T)]
```


```python
fig = plt.figure()
ax = fig.add_subplot(1,1,1)

ax.plot(k_path_ext[:T], lw=2, label='Solow path')
ax.plot(k_stars_ext,lw=2, label="Steady state")
ax.set_title('Simulated convergence of capital to steady state')
ax.set_xlabel('t')
ax.set_ylabel('Capital per capita')
ax.set_xlim([0,100])
ax.set_ylim([0,10])
ax.legend(loc="lower right")
ax.grid(True);
```

We find that the convergence is little bit faster in the extended Solow model with economic damage from climate change. 

# Conclusion

We have solved the basic Solow model and an extended Solow model, which account for economic damage from climate change. We have analytically solved the model using Pythons Sympy and numerically using two different methods.

For the classical Solow model, we find that in steady-state, the capital per effective worker equals 7.235. The result was constant despite the different approaches. Furthermore, we find that The GDP per effective worker is 1.934.

For the extended model, we conclude that the capital per effective worker equals 5.421 in the analytical and numerical solution, and the GDP per effective worker is 1.449.

Comparing the two models, we find that for any given level of D the reduction in GDP per capita is greater than the reduction in production due to the negative effect from capital accumulation.

578eb426a64cfeae36aacbf799080f45bec4c26e
ipynb
Jupyter Notebook
modelproject/modelproject.ipynb
AskerNC/projects-2021-sstl
f64ab9fb8bd398c245b296e156cf5f5acaf53fa7
MIT
modelproject/modelproject.ipynb
AskerNC/projects-2021-sstl
f64ab9fb8bd398c245b296e156cf5f5acaf53fa7
MIT
modelproject/modelproject.ipynb
AskerNC/projects-2021-sstl
f64ab9fb8bd398c245b296e156cf5f5acaf53fa7
MIT
Qwen/Qwen-72B
1. YES
2. YES