# Model Project 

We start by importing necessary packages:


```python
import pandas as pd
import numpy as np

import matplotlib.pyplot as plt
import ipywidgets as widgets
import time
from scipy import linalg
from scipy import optimize
import sympy as sm
```

# Solow model with climate change 

Consider the standard Solow-model with an indicator for climate change where:

$K_t$ is capital
    $L_t$ is labor (growing with a constant rate of $n$)
    $A_t$ is technology (growing with a constant rate of $g$)
    $D_t$ is the damage function (growing with a constant rate of $n$)
    $Y_t = F(K_t,A_t,L_t,D_t)$ is GDP

    
The Damage-function depends on temperature T

$$ D_t = 1-\frac{1}{1+\theta_1T_t^{\theta_2}} =< 1 $$

$T_t^{\theta_2}$ is the global mean temperature change at time t. $\theta_1$ is the damage parameter.
    
Cobb-Douglas production-function


$$ F(K_{t},A_{t}L_{t},D_{t})=  (1-D_{t})A_{t}K_{t}^{\alpha}L_{t}^{1-\alpha}$$

It is convenient to express all variables of interest in terms of per effective worker, as this can we used as a measure of welfare of sociely. We use small letters to denote per effective worker, so $y=Y/AL$ and $k=K/AL $, which gives: 

$$ f({k}_{t})=(1-D){k}_{t}^{\alpha} $$


Saving is a constant fraction of GDP

$$ S_t = sY_t,\,s\in(0,1) $$


 **Capital accumulates** according to

\\[ K_{t+1}=S_{t}+(1-\delta)K_{t}=sF(K_{t},A_{t}L_{t})+(1-\delta)K_{t}, \delta \in (0,1) \\]

The economy is growing over time (due to exogenous technological progress and population growth), which makes it useful to focus on the behavior of capital stock per unit of effective labor. The transition in the model can be decribed by


\\[ \tilde{k}_{t+1}= \frac{1}{(1+n)(1+g)}[sf(\tilde{k}_{t})+(1-\delta)\tilde{k}_{t}] \\]




# Level of output for changes in K at different levels of Damage

The output level depends on the damage. Therefore, different levels of damage will lead to different levels of output.

Define parameters and arrays


```python
K_size = 100                # Model domain
alpha = (1/3)                
D = 0.1

k = np.arange(K_size)       
        
```

Output values


```python
def output(k, D):           
    y = (1-D)  * (k)**(alpha)
    return y

y  = output(k, D)
y2 = output(k, D + 0.1)
y3 = output(k, D + 0.2)
y4 = output(k, D + 0.3)
y5 = output(k, D + 0.4)
y6 = output(k, D + 0.5)
```

Plot the production function for different values of damage


```python
fig, ax = plt.subplots(figsize=(12, 12))
ax.set(title="Output", xlabel="Capital", ylabel="Output")
ax.set_ylim([0,20])
ax.grid()
ax.plot(k, y,  "b-", alpha=1.00, label="D= 0.1")
ax.plot(k, y2, "b-", alpha=0.85, label="D= 0.2")
ax.plot(k, y3, "b-", alpha=0.70, label="D= 0.3")
ax.plot(k, y4, "b-", alpha=0.55, label="D= 0.4")
ax.plot(k, y5, "b-", alpha=0.40, label="D= 0.5")
ax.plot(k, y6, "b-", alpha=0.25, label="D= 0.6")
ax.legend() 
ax.set_ylim(0, 6)  
ax.set_xlim(0, 100)
plt.show()
```

With an increase in damage the output decreases. The level of damage seems to have an immense  effect on the output. 

# Graph: The Steady-State (Model with Climate Change)

In the following cells we are going to code the graph for our solow model with climate change. The graph shows the output, investments (=savings) and depreciation in our model. 

Define parameters and arrays


```python
K_size = 101                     
A = 1                           
L = 5                            
alpha = (1/3)                    
s = 0.2                       
delta = 0.10                         
g = 0.02
n = 0.01
D_ = 0.1


K = np.arange(K_size)            # Create empty array of K
```

Define Functions


```python
def output(K):   # Cobb-Douglas Production Function
    Y = (1-D_) * A * (K)**(alpha) * (L)**(1-alpha)    
    return Y
```

Populate arrays


```python
Y = output(K)
D = (delta)*K
S = s*Y
```

Steady State Values

In a closed economy with no government, consumption per capita is the difference between income and investment: $$c=y−i$$

There is a capital depreciation: $$ d_t = {\delta}*k_t $$




```python
Kstar = ((delta+(g*n)+n+g)/s*(1-D_))**(1/(alpha-1))
Ystar = (1-D_)*Kstar**alpha
Sstar = s*Ystar
Cstar = Ystar - Sstar
Dstar = (delta)*Kstar
```

Plot the graph


```python

fig, ax = plt.subplots(figsize=(10, 8))
ax.plot(K, Y, "k", ls = '-', label="Output")
ax.plot(K, S, "b", ls = '-', label="Investment")
ax.plot(K, D, "r", ls = '-', label="Depreciation")
ax.set(title="Solow Model", xlabel="Capital Stock")
plt.text(77, 12,  r'$Y = A \cdot K^{\alpha} L^{1-\alpha}$')
plt.text(90, 10,  r'$D = {\delta}K$')
plt.text(90, 3, r'$S = sY$')
plt.legend(loc=2)
ax.set_xlim(0, 100)
ax.set_ylim(0, 15)
plt.show()
```

# Behavior of capital stock per unit of effective labor over time: Comparing Model with and without Climate Change

## Model including climate change 


```python
k = sm.symbols('k')
alpha = sm.symbols('alpha')
delta = sm.symbols('delta')
s = sm.symbols('s')
A = sm.symbols('A')
D = sm.symbols('D')
n = sm.symbols('n')
g = sm.symbols('g')

sm.init_printing(use_unicode=True)
```

Equation of capital stock per unit of effective labour


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

Create a vector for time, t 


```python
x_vec = np.zeros((100,1))
```


```python
x_vec[0] = 1
```

Then we select some parameter values, Let
$$s=0.2 $$
$$ g=0.02 $$ 
$$ n=0.01 $$
$$ delta=0.1$$
$$ alpha=1/3 $$
$$ D=0,1 $$


```python
s = 0.2
g = 0.02
n = 0.01
delta = 0.1
alpha = (1/3)
D = 0.1
```

Now, plug vector x_vec into the capital transition equation: 


```python
for i in range(1, 100):
    x_vec[i] = (s*(1-D)*x_vec[i-1]**alpha+(1-delta)*x_vec[i-1])/((1+g)*(1+n))
```

Now we can check whether the x_vec [steady state capital] converges to a steady state.


```python
x_vec
```




    array([[1.        ],
           [1.04834013],
           [1.09334207],
           [1.13516079],
           [1.17396056],
           [1.20991008],
           [1.24317873],
           [1.27393383],
           [1.30233866],
           [1.328551  ],
           [1.35272218],
           [1.3749964 ],
           [1.39551038],
           [1.4143932 ],
           [1.43176626],
           [1.44774341],
           [1.46243112],
           [1.47592874],
           [1.48832876],
           [1.49971719],
           [1.51017383],
           [1.51977267],
           [1.5285822 ],
           [1.53666575],
           [1.54408184],
           [1.5508845 ],
           [1.55712355],
           [1.56284493],
           [1.56809096],
           [1.57290059],
           [1.5773097 ],
           [1.58135125],
           [1.58505559],
           [1.58845058],
           [1.59156185],
           [1.59441292],
           [1.5970254 ],
           [1.59941913],
           [1.60161232],
           [1.60362169],
           [1.60546255],
           [1.60714899],
           [1.60869391],
           [1.61010912],
           [1.6114055 ],
           [1.61259299],
           [1.6136807 ],
           [1.61467699],
           [1.61558955],
           [1.61642538],
           [1.61719092],
           [1.61789208],
           [1.61853426],
           [1.61912241],
           [1.61966108],
           [1.62015441],
           [1.62060623],
           [1.62102003],
           [1.62139899],
           [1.62174605],
           [1.62206389],
           [1.62235497],
           [1.62262154],
           [1.62286567],
           [1.62308924],
           [1.62329398],
           [1.62348148],
           [1.62365319],
           [1.62381044],
           [1.62395444],
           [1.62408632],
           [1.62420708],
           [1.62431768],
           [1.62441896],
           [1.6245117 ],
           [1.62459664],
           [1.62467442],
           [1.62474564],
           [1.62481087],
           [1.6248706 ],
           [1.6249253 ],
           [1.6249754 ],
           [1.62502127],
           [1.62506328],
           [1.62510174],
           [1.62513697],
           [1.62516923],
           [1.62519877],
           [1.62522583],
           [1.6252506 ],
           [1.62527329],
           [1.62529406],
           [1.62531309],
           [1.62533051],
           [1.62534646],
           [1.62536107],
           [1.62537445],
           [1.6253867 ],
           [1.62539792],
           [1.6254082 ]])



## Model ignoring climate change


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

Create a vector to plug in for time, t


```python
y_vec = np.zeros((100,1))
```


```python
y_vec[0] = 1
```

plug vector y_vec in the equation for behavior of capital stock per unit of effective labour.


```python
for i in range(1, 100):
    y_vec[i] = (s*y_vec[i-1]**alpha+(1-delta)*y_vec[i-1])/((1+g)*(1+n))
```

### Plot both functions in the same graph


```python
fid, ax = plt.subplots()
ax.plot(x_vec, label="model including climate change")
ax.plot(y_vec, label="model ignoring climate change")
ax.set(xlabel = 'time, t', ylabel ='k(t)')
plt.title('capital stock per unit of effective labour over time')
plt.legend()
plt.show()
```

Both equations converge to a steady state over time. The model including climate change has a lower steady state because it implements climate change as a external cost.

# Steady state calculations 

We now want to find an **analytic an analytical** expression for the steady state, that is: where $ k_{t+1} = k_{t} =k_{ss}$


```python

k = sm.symbols('k')
alpha = sm.symbols('alpha')
delta = sm.symbols('delta')
s = sm.symbols('s')
A = sm.symbols('A')
D = sm.symbols('D')
n = sm.symbols('n')
g = sm.symbols('g')
```


```python
sm.init_printing(use_unicode=True)

f = (1-D)*((k)**alpha)
ss = sm.Eq(k,(s*f+(1-delta)*k)/((1+n)*(1+g)))
kss = sm.solve(ss,k)[0]
kss

```

In the standard Solow model, **without the damage-function**, the steady state capital is given by: 


```python
f1= (k)**alpha
ss1 = sm.Eq(k,(s*f1+(1-delta)*k)/((1+n)*(1+g)))
kss1 = sm.solve(ss1,k)[0]
kss1
```

From these equations it becomes clear that, with D≤1, the steady state capital growth is reduced by including the damage-function, compared to the the model without

We can also find the **specific numerical value** for steady state capital. Usning the same parameter values as before, we have 

$$s=0.2 $$
$$ g=0.02 $$ 
$$ n=0.01 $$
$$ delta=0.1$$
$$ alpha=1/3 $$
$$ D=0,1 $$


```python
ss_func = sm.lambdify((s,g,n,delta,alpha,D),kss)
```

Steady state capital, including the damage function, is then given by:


```python
ss_func(0.2,0.02,0.01,0.1,1/3, 0.1)
```

By comparison, excluding the damage function, the steady state is given by:


```python
ss_func(0.2,0.02,0.01,0.1,1/3, 0)
```

Notice how these values corresponds to steady state levels the figure. 

It becomes clear that the steady state capital excluding the damage function (representing climate change) is higher than the steady state capital in case of climate damage. Climate change is a negative externality for the society. For example, increased temperature causes sea level to rise and drought periods. This again, can increase human conflicts. Today's generation benefits from the use of fossil fuels, which releases greenhouse gases into the atmosphere and in turn causes climate change. But the benefit from using fossil fuels, does not internalize these external costs of climate change. 

# Impact of damage on steady state capital 

In order to visualize the impact of damage due to climate change (ceteris paribus), we plot the steady state as a dependent variable of D.

First, create a list that contains all values we want D to take


```python
Ds = [0, 0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]
```

Define the function that returns the steady state only dependent on D.



```python
def calc_ss(D):
    return ss_func(0.2,0.02,0.01,0.1,1/3, D)

```

Create a while loop so that a new list (nDs)is generated that contains all steady states for the respective D.



```python
index = 0
while index <= len(Ds)-1:
    Ds[index] = calc_ss(Ds[index])
    index = index +1

print(Ds)
```

    [1.903831539231319, 1.6255198652134988, 1.3622709566173723, 1.1150018222907796, 0.8848209414512797, 0.6731060958136442, 0.4816355156188146, 0.31283144391803264, 0.17028386957717148, 0.06020443945235181]



```python
Damage = [0,0.1, 0.2, 0.3, 0.4, 0.5, 0.6, 0.7, 0.8, 0.9]
```

## Plot: 


```python
fid, ax = plt.subplots()
ax.plot(Damage, Ds)
ax.set(xlabel = 'Damage', ylabel ='Steady state')
plt.show()
```

From the plot we can see that as damage increases, the steady state capital decreases. This could be due to higher expenses in order to compensate for i.e. agricultural damage.


```python

```


```python

```

d4e2c821fcb8279c6e902b26709f36d825a50324
ipynb
Jupyter Notebook
modelproject/modelproject.ipynb
NumEconCopenhagen/projects-2019-tba
eea677b973b0205f293272027623ca3c13a3c23e
MIT
modelproject/modelproject.ipynb
NumEconCopenhagen/projects-2019-tba
eea677b973b0205f293272027623ca3c13a3c23e
MIT
2019-04-08T17:01:11.000Z
2019-05-14T18:47:37.000Z
modelproject/modelproject.ipynb
NumEconCopenhagen/projects-2019-tba
eea677b973b0205f293272027623ca3c13a3c23e
MIT
2019-03-22T14:44:02.000Z
2019-03-22T14:44:26.000Z
Qwen/Qwen-72B
1. YES
2. YES