# Rasmus 
YOUR PROJECT TITLE

> **Note the following:** 
> 1. This is *not* meant to be an example of an actual **model analysis project**, just an example of how to structure such a project.
> 1. Remember the general advice on structuring and commenting your code from [lecture 5](https://numeconcopenhagen.netlify.com/lectures/Workflow_and_debugging).
> 1. Remember this [guide](https://www.markdownguide.org/basic-syntax/) on markdown and (a bit of) latex.
> 1. Turn on automatic numbering by clicking on the small icon on top of the table of contents in the left sidebar.
> 1. The `modelproject.py` file includes a function which could be used multiple times in this notebook.

Imports and set magics:


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

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

# local modules
import modelproject
```

    The autoreload extension is already loaded. To reload it, use:
      %reload_ext autoreload


The general Solow model consist of the followning equations:

\begin{align}
(1)&\:\:\:Y_t = K_t^\alpha (A_t L_t)^{1-\alpha}\\
(2)&\:\:\:K_{t+1} = s Y_t +(1+\delta)K_t\\
(3)&\:\:\:L_{t+1} = (1+n)L_t\\
(4)&'\:\:\:A_{t+1} = (1+g)A_t$
\end{align}


$Y_t$, $K_t$, $A_t$, and $L_t$ are the production, capital, technology level, and population in year $t$ respectively. $s$ are the saving rate in physical and human capital respectively. $\alpha$ contribution variables of production and $n$ and $g$ are growth in tecnology and population respectively. Finally we have $\delta$ as a paramater for the depreciation.

Since we have that capital is both accumulating and depreciating, there exist a stable level of capital, known as the steady state. But because the tecknology level and population is constantly growing, the steady state exist in a tecnology adjusted capital per capita. This translate into a steady state in technology adjusted BNP per capita. We define these terms as

\begin{align}
\tilde{y}\equiv \frac{Y_t}{A_t L_t}\\            
\tilde{k}\equiv \frac{K_t}{A_t L_t}
\end{align}

By dividing both sides of $Y_t = K_t^\alpha (A_t L_t)^{1-\alpha}$ with $A_t L_t$ we get that $\tilde{y}_t = \tilde{k}^\alpha_t$. If we insert equation (1),(3), and (4) into (2), we get a transistion equation for $\tilde{k}$:

\begin{equation}
    \tilde{k}_{t+1} = \frac{s \tilde{k}_t^\alpha + (1-\delta) \tilde{k}_t}{(1+n)(1+g)}
\end{equation}

I start out by solving the steady states for $\tilde{y}$ and $\tilde{k}$. First symbolically


```python
ktilde = sm.symbols('k_tilde')
alpha = sm.symbols('alpha')
delta = sm.symbols('delta')
s = sm.symbols('s')
g = sm.symbols('g')
n = sm.symbols('n')
```


```python
ktilde1 = (s*ktilde**(alpha) + (1-delta)*ktilde)/( (1+n)*(1+g) )
kss = sm.solve(sm.Eq(ktilde,( s*ktilde**(alpha) + (1-delta)*ktilde)/( (1+n)*(1+g) )),ktilde)[0]
kss
```




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




```python
ytilde = ktilde**alpha
yss =ytilde.subs(ktilde,kss)
yss
```




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




```python

```


```python
ss_ktilde_func = sm.lambdify((s,g,n,delta,alpha),kss)
ss_ytilde_func = sm.lambdify((s,g,n,delta,alpha),yss)
ktilde1_func = sm.lambdify((ktilde,s,g,n,delta,alpha),ktilde1)
```

We could also have solved steady states it numerically using a rootfinder


```python
alpha =0.3
delta =0.05
s = 0.2
n = 0.01
g = 0.02
x0 = 1

def ssk_num_func(s,g,n,delta,alpha):
    def obj(ktilde):
        return ktilde-(s*ktilde**(alpha) + (1-delta)*ktilde)/((1+n)*(1+g))
    result = optimize.root(obj, x0)
    return result.x[0]

def ssy_num_func(s,g,n,delta,alpha):
    return ssk_num_func(s,g,n,delta,alpha)**alpha
print(ssk_num_func(s,g,n,delta,alpha))
print(ss_ktilde_func(s,g,n,delta,alpha))
print(ssy_num_func(s,g,n,delta,alpha))
print(ss_ytilde_func(s,g,n,delta,alpha))

```

    3.689237466703929
    3.689237466703934
    1.4793842241482769
    1.4793842241482773


# transition


```python
ktilde_vec = np.linspace(0,5,100)

trans = ktilde1_func(ktilde_vec,s,g,n,delta,alpha)
plt.plot(ktilde_vec, trans)
plt.plot(ktilde_vec,ktilde_vec)
```


```python

```


```python
ktilde_vec = np.linspace(0,7,10)
trans = ktilde1fun(ktilde_vec,0.2,0.3,0.05,0.01,0.02)
print(trans)
plt.plot(ktilde_vec,trans)
plt.plot(ktilde_vec,ktilde_vec)
plt.axvline(x=ss_ktilde_func(0.2,0.02,0.01,0.05,0.3))
```


```python
def tjeckconc(alpha):
    def plot_ss_func(n,g,delta,s):
        ktilde_vec = np.linspace(0,ss_ktilde_func(s,g,n,delta,alpha),100)
        def depriciation(ktilde,n,g,delta):
            return (n+g+delta+n*g)*ktilde
        def growthk(ktilde,s,alpha):
            return s*ktilde**alpha
        s = 0.2
        depri = depriciation(ktilde_vec,n,g,delta)
        grow = growthk(ktilde_vec,s,alpha)
        plt.plot(ktilde_vec,depri)
        plt.plot(ktilde_vec, grow)
    widgets.interact(plot_ss_func,
    n=widgets.FloatSlider(description="n", min=0.005, max=0.6, step=0.005, value=0.01),
    g=widgets.FloatSlider(description="g", min=0.005, max=0.6, step=0.005, value=0.02),
    delta=widgets.FloatSlider(description="$\delta$", min=0.01, max=0.6, step=0.01, value=0.05),
    s=widgets.FloatSlider(description="$s$", min=0.05, max=0.6, step=0.01, value=0.2),
    );

tjeckconc(alpha)
```


    interactive(children=(FloatSlider(value=0.01, description='n', max=0.6, min=0.005, step=0.005), FloatSlider(va…



```python
def plot_ss(g,n,delta,alpha):
    def plot_ss_func(s,new_s):
        max_x = max(ss_ktilde_func(s,g,n,delta,alpha),ss_ktilde_func(new_s,g,n,delta,alpha))
        ktilde_vec = np.linspace(0,max_x*1.2,100)
        def depriciation(ktilde,n,g,delta):
            return (n+g+delta+n*g)*ktilde
        def growthk(ktilde,s,alpha):
            return s*ktilde**alpha
        depri = depriciation(ktilde_vec,n,g,delta)
        grow = growthk(ktilde_vec,s,alpha)
        plt.plot(ktilde_vec,depri)
        plt.plot(ktilde_vec, grow)
        plt.axvline(x=ss_ktilde_func(s,g,n,delta,alpha), linestyle='--', color = 'grey')

        grow = growthk(ktilde_vec,new_s,alpha)
        plt.plot(ktilde_vec, grow)
        plt.axvline(x=ss_ktilde_func(new_s,g,n,delta,alpha), linestyle='--', color = 'black')

    widgets.interact(plot_ss_func,
    s=widgets.FloatSlider(description="s", min=0.05, max=0.6, step=0.02, value=0.15),
    new_s=widgets.FloatSlider(description="New s", min=0.05, max=0.6, step=0.02, value=0.25),
    );

plot_ss(0.02,n,delta,alpha)
```


    interactive(children=(FloatSlider(value=0.15, description='s', max=0.6, min=0.05, step=0.02), FloatSlider(valu…



```python

```


```python

```

# Increase in saving rate

As $\tilde{y}^*$ is constant, we have that $y^*$ is growing at at contanst rate along with $A_t$, $g$, $\tilde{y}^* = A_t y^* = A_0(1+g)^t y^*$ insteady state. If some of the parameter chance, the growth in $y$ will change until the economy again has found its steady state. If 


```python

```


```python

```


```python
#initial parameters
alpha =0.3
delta =0.05
s = 0.2
n = 0.01
g = 0.02
A0 = 10   
def plotconvergens(new_s,new_s_period=5,T=50):
    
    y_vec = []
    ktilde0 = ss_ktilde_func(s,g,n,delta,alpha)
    y0 = ss_ytilde_func(s,g,n,delta,alpha)*A0

    for i in range(0,new_s_period):
        y_vec.append(y0*(1+g)**i)

    ktilde = ktilde0
    for i in range(new_s_period,T):
        ktilde = ktilde1_func(ktilde,new_s,g,n,delta,alpha)
        y = ktilde**alpha*A0*(1+g)**i
        y_vec.append(y)

    ln_y = np.log(y_vec)
    g_y_vec = np.diff(ln_y)

    fig = plt.figure(dpi = 130) 
    ax = fig.add_subplot(2,1,1)
    ax.plot(ln_y)

    ax = fig.add_subplot(2,1,2)
    plt.plot(g_y_vec)
    plt.show()                
```


```python
widgets.interact(plotconvergens,
    #alphas=widgets.fixed(new_s_period=5,T=100),
    new_s=widgets.FloatSlider(description="New s", min=0.01, max=0.99, step=0.05, value=0.25),
    new_s_period=widgets.fixed(5),
    T=widgets.fixed(50)
);
```


    interactive(children=(FloatSlider(value=0.25, description='New s', max=0.99, min=0.01, step=0.05), Output()), …


# Model description

The Solow model with human capital consist of the followning equations:

\begin{align}
Y_t &= K_t^\alpha  H_t^\varphi (A_t L_t)^{1-\alpha-\varphi}\\
K_{t+1} &= s_K Y_t +(1+\delta)K_t\\
H_{t+1} &= s_H Y_t + (1+\delta)H_t\\
L_{t+1} &= (1+n)L_t\\
A_{t+1} &= (1+g)A_t
\end{align}


$Y_t$, $K_t$, $H_t$, $A_t$, and $L_t$ are the production, physical capital, human capital, technology level, and population in year $t$ respectively. $s_K$ and $s_H$ are the saving rate in physical and human capital respectively. $\alpha$ and $\varphi$ are contribution variables of production and $n$ and $g$ are growth in tecnology and population respectively. Finally we have $\delta$ as a paramater for the depreciation in both physical and human capital.

$\tilde{y}\equiv \frac{Y_t}{A_t L_t}$              $\tilde{k}\equiv \frac{K_t}{A_t L_t}$              $\tilde{h}\equiv \frac{H_t}{A_t L_t}$

## 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. 




```python
# Defining symbols for sympy 

# period t
K = sm.symbols("K_t")
H = sm.symbols("H_t")
A = sm.symbols("A_t")
L = sm.symbols("L_t")

phi = sm.symbols("varphi")
sK = sm.symbols("s_K")
sH = sm.symbols("s_H")
alpha = sm.symbols('alpha')
delta = sm.symbols('delta')
varphi = sm.symbols('varphi')
g = sm.symbols('g')
n = sm.symbols('n')

kstar = sm.symbols("ktilde^*")
hstar = sm.symbols("htilde^*")
ystar = sm.symbols("ytilde^*")

ktilde = sm.symbols("ktilde_t")
htilde = sm.symbols("htilde_t")
ytilde = sm.symbols("ytilde_t")
```


```python
Y = K**alpha*H**phi*(A*L)**(1-alpha-phi)

# Physical capital
Kt1 = sK*Y +(1-delta)*K

# Human capital
Ht1 = sH*Y +(1-delta)*H

# Technology level and labor force
At1 = (1+g)*A
Lt1 = (1+n)*L
Y
```




$\displaystyle H_{t}^{\varphi} K_{t}^{\alpha} \left(A_{t} L_{t}\right)^{- \alpha - \varphi + 1}$



transition equations


```python

ktilde1 = ( sK*ktilde**alpha*htilde**varphi + (1-delta)*ktilde) /((1+n)*(1+g))
htilde1 = ( sH*ktilde**alpha*htilde**varphi + (1-delta)*htilde) /((1+n)*(1+g))
```


```python

```


```python
redundant = ((g+1)*(n+1))
ksolow = sm.simplify(ktilde1-ktilde)
ksolowss = ksolow*redundant
ssk = sm.simplify(sm.Eq(ksolowss,0))
ssk
```




$\displaystyle - \tilde{h}_{t}^{\varphi} \tilde{k}_{t}^{\alpha} s_{K} + \tilde{k}_{t} \left(\delta - 1\right) + \tilde{k}_{t} \left(g + 1\right) \left(n + 1\right) = 0$




```python
hsolow = sm.simplify(htilde1-htilde)
hsolowss = hsolow*redundant
ssh = sm.Eq(hsolowss,0)
ssh
```




$\displaystyle - \tilde{h}_{t} \left(\delta - 1\right) - \tilde{h}_{t} \left(g + 1\right) \left(n + 1\right) + \tilde{h}_{t}^{\varphi} \tilde{k}_{t}^{\alpha} s_{H} = 0$




```python

# Finding steady state for physical capital
ktildesub = sm.solve(ssh, ktilde)[0]
ktildesub
#kSS = sm.Eq(kstar, sm.solve(ssk, ktilde)[0])
#kSS
```




$\displaystyle \left(\frac{\tilde{h}_{t}^{1 - \varphi} \left(\delta + g n + g + n\right)}{s_{H}}\right)^{\frac{1}{\alpha}}$




```python
hss1 = ssh.subs(ktilde, ktildesub)
sm.simplify(hss1)
n
```




$\displaystyle \tilde{h}_{t} \left(\delta - 1\right) + \tilde{h}_{t} \left(g + 1\right) \left(n + 1\right) - \tilde{h}_{t}^{\varphi} s_{H} \left(\left(\frac{\tilde{h}_{t}^{1 - \varphi} \left(\delta + g n + g + n\right)}{s_{H}}\right)^{\frac{1}{\alpha}}\right)^{\alpha} = 0$




```python
sej = sm.solve(hss1)
sej
scipy.optimize.root
```

## Numerical solution

You can always solve a model numerically. 

Define first the set of parameters you need. 

Then choose one of the optimization algorithms that we have gone through in the lectures based on what you think is most fitting for your model.

Are there any problems with convergence? Does the model converge for all starting values? Make a lot of testing to figure these things out. 


```python
alpha =0.3
varphi = 0.3
delta =0.05
sH = 0.2
sK = 0.3
n = 0.01
g = 0.02

htilde = 13.3

def obj(ktilde, alpha, sK, sH, varphi, delta, n, g):
    return ktilde-(sK*ktilde**alpha*htilde**varphi + (1-delta)*ktilde) /((1+n)*(1+g))


```

# Further analysis


```python

x_vec = np.linspace(0.0000000001,20,100)
obj(x_vec, alpha, sK, sH, varphi, delta, n, g)
res = optimize.root(obj,20,args=(alpha, sK, sH, varphi, delta, n, g),method='hybr')
res.x
```




    array([19.95934363])




```python

```

Make detailed vizualizations of how your model changes with parameter values. 

Try to make an extension of the model. 

# Conclusion

Add concise conclusion. 


```python
# The complete model

def funYt(Kt,Ht,At,Lt, alpha, phi):
    Yt = Kt**alpha * Ht**phi * (At*Lt)**(1-alpha-phi)
    return Yt

def funKt1(sK, Yt, delta, Kt):
    Kt1 = sK*Yt - delta*Kt + Kt
    return Kt1

def funHt1(sH, Yt, delta, Ht):
    Ht1 = sH*Yt - delta*Ht + Ht
    return Ht1

def funLt1(Lt,n):
    Lt1 = (1+n)*Lt
    return Lt1

def funAt1(At, g):
    At1 = (1+g)*At
    return At1
```


```python
def funkttilde(Kt,At,Lt):
    kttilde = Kt/(At*Lt)
    return kttilde

def funhttilde(Ht,At,Lt):
    httilde = Ht/(At*Lt)
    return httilde

def funyttilde(kttilde, httilde, alpha, phi):
    yttilde = kttilde**alpha * httilde**phi
    return yttilde

def funkttilde1(yttilde, sK, delta, ktilde, n, g):
    kttilde1 = (sK*yttilde + (1-delta)*kttilde) / ((1+n)+(1+g))
    return kttilde1

def funhttilde1(yttilde, sH, delta, htilde, n, g):
    httilde1 = (sH*yttilde + (1-delta)*httilde) / ((1+n)+(1+g))
    return httilde1

def funobj_sskttilde(kttilde, kttilde1):
    res = kttilde-kttilde1
    return res
```


```python
Yt = 10
Kt = 10
Ht = 20
At = 8
Lt = 100
alpha =0.3
phi = 0.3
delta =0.05
sH = 0.2
sK = 0.3
n = 0.01
g = 0.02
```


```python
httilde = funkttilde(Kt,At,Lt)  
kttilde = funhttilde(Ht,At,Lt)
yttilde = funyttilde(kttilde, httilde, alpha, phi)
```


```python
kttilde1 = funkttilde1(yttilde, sK, delta, kttilde, n, g)
funobj_sskttilde(kttilde, kttilde1)
```




    0.013744841152236894




```python
def deltak(sK, alpha, phi, n, g, delta, httilde):
    return ( (n+g+delta+n*g)/(sK*httilde**phi) )**(1/(alpha-1))

def deltah(sH, alpha, phi, n, g, delta, kttilde):
    return ( (n+g+delta+n*g)/(sH*kttilde**alpha) )**(1/(phi-1))
```


```python
httilde = np.linspace(0.01,30,1000)
kttilde = np.linspace(0.01,30,1000)
```


```python
h_vec = deltak(sK, alpha, phi, n, g, delta, httilde)
k_vec = deltah(sH, alpha, phi, n, g, delta, kttilde)
plt.plot(httilde, h_vec)
plt.plot(k_vec,httilde)
```


```python
kstar = ( (sK**(1-phi)*sH**(phi))/(n+g+delta+n*g))**(1/(1-alpha-phi))
hstar = ( (sK**(alpha)*sH**(1-alpha))/(n+g+delta+n*g))**(1/(1-alpha-phi))
print(hstar, kstar)
```

    13.310902819474649 19.966354229211966


Microtest


```python
from types import SimpleNamespace
par = SimpleNamespace()
par.alpha = 0.5
par.w = 10
par.L = 10
par.p1 = 5
par.beta = 0.8


def util(x1):
    x2 = budget(par)-par.p1*x1
    return -x1**par.alpha*x2**(1-par.alpha)

def budget(par):
    return par.w*par.L

def budgetcon(x1):
    x2 = budget(par)-par.p1*x1
    return budget(par)-par.p1*x1+-x2

```


```python
def demand(par):
    result = optimize.minimize(util,6,
                               method='Nelder-Mead',
                                )
    return result.x
demand(par)
```




    array([9.99997559])




```python
x1_vec = np.linspace(0,budget(par)/par.p1,1000)
plt.plot(x1_vec,-util(x1_vec))
```


```python
x2 = budget(par)-par.p1*x1
print(budget(par))
print(budgetcon(x1))
print(-util(x1))
```

    100
    0
    19.364916731037088



```python
def production(L):
    return par.L**par.beta*20

def profit(par):
    Pi = par.p1*demand(par)-par.L*par.w
    return Pi
```


```python
print(production(L))
profit(par)
```

    126.19146889603867





    array([-50.00012207])




```python

```

25b314029e120c6efeaef15b632e41f190527ca6
ipynb
Jupyter Notebook
modelproject/modelproject.ipynb
NumEconCopenhagen/projects-2022-ludvigsen
b19d3081d988d41e0acf8be0903e28043a7762bc
MIT
modelproject/modelproject.ipynb
NumEconCopenhagen/projects-2022-ludvigsen
b19d3081d988d41e0acf8be0903e28043a7762bc
MIT
2022-03-29T12:16:49.000Z
2022-03-31T11:53:27.000Z
modelproject/modelproject.ipynb
NumEconCopenhagen/projects-2022-ludvigsen
b19d3081d988d41e0acf8be0903e28043a7762bc
MIT
Qwen/Qwen-72B
1. YES
2. YES