# Model Project - OLG


```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

```

# Model description



In this project we will first consider a an OLG model with a pay as you go system which we will later compare to an OLG model with a fully funded system.

Individuals live for two periods and we assume constant population:

$L_t=L_{t+1}$

The lifetime utility of a household is defined as:

$U_t=ln(c_{1t})+\frac{1}{(1+\rho)}ln(c_{2t+1}), \ \ \rho>-1$

This is subject to the following constraints:

$c_{1t}=w_t(1-\tau)-s_t$

$c_{2t+1}=s_t(1+r_{t+1})+\tau w_{t+1}$

Thus the model has the following parameters:

$L_t$: The population at time $t$

$U_t$: Lifetime household utility

$\rho$: Consumption discount rate

$c_{1t}$: Consumption of the young

$c_{2t+1}$: Consumption of the old

$w_t$: Wage

$r_t$: Return on capital

$s_t$: Savings


```python
# Define symbols for sympy
K_t = sm.symbols('K_t')
L_t = sm.symbols('L_t')
A = sm.symbols('A')
alpha = sm.symbols('alpha')
delta = sm.symbols('delta')
rho = sm.symbols('rho')
w_t = sm.symbols('w_t')
r_t = sm.symbols('r_t')
k_t = sm.symbols('k_t')
k_plus = sm.symbols('k_t+1')
Pi = sm.symbols('Pi_t')

s_t = sm.symbols('s_t')
U_t = sm.symbols('U_t')
U = sm.symbols('U_t')
r_plus = sm.symbols('r_t+1')
k = sm.symbols('k^*')

d_t = sm.symbols('d_t')
d_plus = sm.symbols('d_t+1')
tau = sm.symbols('tau')
s_payg = sm.symbols('s_payg')
w_plus = sm.symbols('w_t+1')
k_payg = sm.symbols('k_payg')
U_payg = sm.symbols('U_paygt')
kpaygss = sm.symbols('k_payg^*')
```


```python
#In order to derive the wage (w) and return on capital (r) we need to solve the firms profit maximization problem.

#Define the profit function
profit = A*(K_t**alpha)*(L_t**(1-alpha))-w_t*L_t-(r_t+delta)*K_t

#Setup the FOC's
foc1=sm.diff(profit,K_t)
foc2=sm.diff(profit,L_t)

#Solve FOC's
solve_foc1=sm.solve(foc1,r_t)[0]
solve_foc2=sm.solve(foc2,w_t)[0]

#Print eqations
print('The profit of a given firm is defined as:')
display(sm.Eq(Pi,profit))
print('By taking the first order conditions to this, we get the expressions for w an r:')
display(sm.Eq(r_t,solve_foc1))
display(sm.Eq(w_t,solve_foc2))
```

    The profit of a given firm is defined as:



$\displaystyle \Pi_{t} = A K_{t}^{\alpha} L_{t}^{1 - \alpha} - K_{t} \left(\delta + r_{t}\right) - L_{t} w_{t}$


    By taking the first order conditions to this, we get the expressions for w an r:



$\displaystyle r_{t} = \frac{A K_{t}^{\alpha} L_{t} L_{t}^{- \alpha} \alpha}{K_{t}} - \delta$



$\displaystyle w_{t} = A K_{t}^{\alpha} L_{t}^{- \alpha} \left(1 - \alpha\right)$


We now assume full capital depreciation $\delta=1$


```python
#Setup equations
sample_r=alpha*A*k_t**(alpha-1)-1
sample_w=(1-alpha)*A*k_t**alpha

#Display them
print('By using the marginal products of capital and labour we can simplify these expressions to:')
display(sm.Eq(r_t,sample_r))
display(sm.Eq(w_t,sample_w))
```

    By using the marginal products of capital and labour we can simplify these expressions to:



$\displaystyle r_{t} = A \alpha k_{t}^{\alpha - 1} - 1$



$\displaystyle w_{t} = A k_{t}^{\alpha} \left(1 - \alpha\right)$


We now move unto the household optimization problem and determine the savings rate.


```python
#We solve the household optimization problem and determine savings (s)

#Substitute constraints into utility function
U_t = sm.log(w_t*(1-tau)-s_t)+(1/(1+rho))*sm.log(s_t*(1+r_plus)+tau*w_plus)

#Take the FOC in regards to savings
FOC_U = sm.diff(U_t,s_t)

#Isolate savings
isolate_s = sm.solve(FOC_U,s_t)[0]

#Print results
print('Substitute the constraints into the utility function')
display(sm.Eq(U,U_t))
print('Take the FOC in regards to savings:')
display(sm.Eq(0,FOC_U))
print('Isolate savings:')
display(sm.Eq(s_t,isolate_s))
print('we can simplify this to:')
display(sm.Eq(s_t,-tau*w_plus*(rho+1)/((r_plus+1)*(rho+2))+w_t*(1-tau)/(rho+2)))

```

    Substitute the constraints into the utility function



$\displaystyle U_{t} = \log{\left(- s_{t} + w_{t} \left(1 - \tau\right) \right)} + \frac{\log{\left(s_{t} \left(r_{t+1} + 1\right) + \tau w_{t+1} \right)}}{\rho + 1}$


    Take the FOC in regards to savings:



$\displaystyle 0 = \frac{r_{t+1} + 1}{\left(\rho + 1\right) \left(s_{t} \left(r_{t+1} + 1\right) + \tau w_{t+1}\right)} - \frac{1}{- s_{t} + w_{t} \left(1 - \tau\right)}$


    Isolate savings:



$\displaystyle s_{t} = \frac{- r_{t+1} \tau w_{t} + r_{t+1} w_{t} - \rho \tau w_{t+1} - \tau w_{t} - \tau w_{t+1} + w_{t}}{r_{t+1} \rho + 2 r_{t+1} + \rho + 2}$


    we can simplify this to:



$\displaystyle s_{t} = - \frac{\tau w_{t+1} \left(\rho + 1\right)}{\left(r_{t+1} + 1\right) \left(\rho + 2\right)} + \frac{w_{t} \left(1 - \tau\right)}{\rho + 2}$


# Steady state

We now wish to determine the steady state level of capital $k^*$


```python
# We define capital accumulation using the previuos expressions and solve for the steady state level of capital.

#capital accumulation 
print('capital accumulation is given by')
display(sm.Eq(k_plus,-tau*w_plus*(rho+1)/((r_plus+1)*(rho+2))+w_t*(1-tau)/(rho+2)))
print('We insert the expressions for w and r and get')
display(sm.Eq(k_plus,(A*k_t**(alpha)*(1-alpha)*(1-tau))/(rho+2)-(k_plus**(1-alpha)*tau*w_plus*(rho+1))/(A*alpha*(rho+2))))
```

    capital accumulation is given by



$\displaystyle k_{t+1} = - \frac{\tau w_{t+1} \left(\rho + 1\right)}{\left(r_{t+1} + 1\right) \left(\rho + 2\right)} + \frac{w_{t} \left(1 - \tau\right)}{\rho + 2}$


    We insert the expressions for w and r and get



$\displaystyle k_{t+1} = \frac{A k_{t}^{\alpha} \left(1 - \alpha\right) \left(1 - \tau\right)}{\rho + 2} - \frac{k_{t+1}^{1 - \alpha} \tau w_{t+1} \left(\rho + 1\right)}{A \alpha \left(\rho + 2\right)}$



```python
capital_accumulation = 1/(1+(1+rho)/(2+rho)*(1-alpha)/alpha*tau)*(((1-alpha)*(1-tau))/(2+rho)*A*k_t**alpha)
capital_accumulation
```




$\displaystyle \frac{A k_{t}^{\alpha} \left(1 - \alpha\right) \left(1 - \tau\right)}{\left(1 + \frac{\tau \left(1 - \alpha\right) \left(\rho + 1\right)}{\alpha \left(\rho + 2\right)}\right) \left(\rho + 2\right)}$



In the steady state it is true that 

$k_{t+1}=k_t=k^*$

We can therefore rewrite the above to:


```python
#Define the steady state capital
ssk=sm.Eq(k,(1/(1+(1+rho)/(2+rho)*(1-alpha)/(alpha)*tau)*((1-alpha)*(1-tau))/(2+rho)*A)**(1/(1-alpha)))
display(ssk)

```


$\displaystyle k^{*} = \left(\frac{A \left(1 - \alpha\right) \left(1 - \tau\right)}{\left(1 + \frac{\tau \left(1 - \alpha\right) \left(\rho + 1\right)}{\alpha \left(\rho + 2\right)}\right) \left(\rho + 2\right)}\right)^{\frac{1}{1 - \alpha}}$


## Analytical solution


```python
#Create a function to determine capital level in steady state given a set of parameters
ssk_solve = sm.solve(ssk,k)
ss_func = sm.lambdify((alpha,tau,rho,A),ssk_solve)

# Assign arbitrary values to the parameters
value_alpha = 0.20
value_tau = 0.25
value_rho = 0.1
value_A = 2

#Run the function with the values
res=ss_func(value_alpha,value_tau,value_rho,value_A)

#Display results
print('Thus we find that the ss capital level, given the parameters, is: '+str(k)+' = '+str(res[0]))
```

    Thus we find that the ss capital level, given the parameters, is: k^* = 0.29345335876374135


We can also check this by doing it manually.


```python
#Define function
def ss_alt(alpha, tau, rho, A):
    #Basically the equation for the steady state capital
    return((1/(1+(1+rho)/(2+rho)*(1-alpha)/(alpha)*tau)*((1-alpha)*(1-tau))/(2+rho)*A)**(1/(1-alpha)))
#Print results with given parameters
print(f'The steady state capital in the pay as you go model is: '+str(ss_alt(value_alpha,value_tau,value_rho,value_A)))

```

    The steady state capital in the pay as you go model is: 0.29345335876374135


# Numerical solution


```python
#Define a lambda function for the capital accumulation function and optimize in regards to k
obj_kpayg = lambda kpayg: kpayg - (1/(1+(1+value_rho)/(2+value_rho)*(1-value_alpha)/value_alpha*value_tau)*(((1-value_alpha)*(1-value_tau))/(2+value_rho)*value_A*kpayg**value_alpha))

r_payg = optimize.root_scalar(obj_kpayg, bracket=[0.1, 100], method='brentq')

#Print results
print(r_payg)
print(f'k^*={r_payg.root:.4f}')
```

          converged: True
               flag: 'converged'
     function_calls: 8
         iterations: 7
               root: 0.2934533587636173
    k^*=0.2935


# OLG Fully Funded system

Next an OLG model with a fully funded system is solved.



# Model description

The lifetime utility of a household is the same as in the pay as you go model:

$U_t=ln(c_{1t})+\frac{1}{(1+\rho)}ln(c_{2t+1}), \ \ \rho>-1$

The constraint now looks as following, the government takes part of their wage and invest it. In period 2 they receive this amount plus the interest earned, plus any additional savings they have made:

$c_{1t}=w_t(1-\tau)-s_t$

$c_{2t+1}=(s_t+w_t\tau)(1+r_{t+1})$

There is still full depreciation and a constant population. The firms maximisation problem does not change. 

The approach for solving this model is the same as for the pay as you go model.


```python
#We solve the household optimization problem and determine savings (s)

#Substitute constraints into utility function
U_t2 = sm.log(w_t*(1-tau)-s_t)+(1/(1+rho))*sm.log((s_t+w_t*tau)*(1+r_plus))

#Take the FOC in regards to savings
FOC_U2 = sm.diff(U_t2,s_t)

#Isolate savings
isolate_s2 = sm.solve(FOC_U2,s_t)[0]

#Print results
print('Substitute the constraints into the utility function')
display(sm.Eq(U,U_t2))
print('Take the FOC in regards to savings:')
display(sm.Eq(0,FOC_U2))
print('Isolate savings:')
display(sm.Eq(s_t,isolate_s2))

#Reducing
print('Reduce this to:')
s2 = w_t*(1/(2+rho)-tau)
s2 = sm.Eq(s_t,s2)
display(s2)

print('Substitute in wage:')
display(s2.subs(w_t, sample_w))

#Substitute in k_t
cap=s2.subs(w_t, sample_w)

#Display the cap accum
print('From this we can get the capital accumulation expression')
display(cap.subs(s_t, k_plus))
```

    Substitute the constraints into the utility function



$\displaystyle U_{t} = \log{\left(- s_{t} + w_{t} \left(1 - \tau\right) \right)} + \frac{\log{\left(\left(r_{t+1} + 1\right) \left(s_{t} + \tau w_{t}\right) \right)}}{\rho + 1}$


    Take the FOC in regards to savings:



$\displaystyle 0 = - \frac{1}{- s_{t} + w_{t} \left(1 - \tau\right)} + \frac{1}{\left(\rho + 1\right) \left(s_{t} + \tau w_{t}\right)}$


    Isolate savings:



$\displaystyle s_{t} = \frac{w_{t} \left(- \rho \tau - 2 \tau + 1\right)}{\rho + 2}$


    Reduce this to:



$\displaystyle s_{t} = w_{t} \left(- \tau + \frac{1}{\rho + 2}\right)$


    Substitute in wage:



$\displaystyle s_{t} = A k_{t}^{\alpha} \left(1 - \alpha\right) \left(- \tau + \frac{1}{\rho + 2}\right)$


    From this we can get the capital accumulation expression



$\displaystyle k_{t+1} = A k_{t}^{\alpha} \left(1 - \alpha\right) \left(- \tau + \frac{1}{\rho + 2}\right)$




# Steady State

In the steady state it is true that 

$k_{t+1}=k_t=k^*$

We can therefore rewrite the above to:


```python
#Define the steady state capital expression
ssk2=sm.Eq(k,((A*(alpha*rho*tau+2*alpha*tau-alpha-rho*tau-2*tau+1))/(rho+2))**(-1/(alpha-1)))
#display it
display(ssk2)

```


$\displaystyle k^{*} = \left(\frac{A \left(\alpha \rho \tau + 2 \alpha \tau - \alpha - \rho \tau - 2 \tau + 1\right)}{\rho + 2}\right)^{- \frac{1}{\alpha - 1}}$


# Analytical Solution


```python
#Create a function to determine capital level in steady state given a set of parameters
ssk_solve2 = sm.solve(ssk2,k)
ss_func2 = sm.lambdify((alpha,tau,rho,A),ssk_solve2)

#Run the function with the values
res2=ss_func2(value_alpha,value_tau,value_rho,value_A)

#Display results
print('Thus we find that the ss capital level, given the parameters, is: '+str(k)+' = '+str(res2[0]))

```

    Thus we find that the ss capital level, given the parameters, is: k^* = 0.28070029686156567


# Nummerical Solution

We now run the nummerical analysis


```python
display(cap.subs(s_t, k_plus))
display((A*k_t**alpha*(1-alpha))*(-tau+1/(rho+2)))
```


$\displaystyle k_{t+1} = A k_{t}^{\alpha} \left(1 - \alpha\right) \left(- \tau + \frac{1}{\rho + 2}\right)$



$\displaystyle A k_{t}^{\alpha} \left(1 - \alpha\right) \left(- \tau + \frac{1}{\rho + 2}\right)$



```python
#Define a lambda function for the capital accumulation function and optimize in regards to k
obj_kff = lambda kff: kff - ((value_A*kff**value_alpha*(1-value_alpha))*(-value_tau+1/(value_rho+2)))

r_ff = optimize.root_scalar(obj_kff, bracket=[0.1, 100], method='brentq')

#Print results
print(r_ff)
print(f'k^*={r_ff.root:.4f}')
```

          converged: True
               flag: 'converged'
     function_calls: 8
         iterations: 7
               root: 0.2807002968610902
    k^*=0.2807


# Comparision of the fully funded system and the payg system

From this we see that in the payg system there was a capital level of 0.293 and in the fully funded it was 0.281. It would therefore seem that with the given parameters that the capital level is higher in the payg-model.

One could furthermore try to change the parameters, like the tax level to observe the changes which we will do now.


```python
# Assign a new tax value (tau)
value_tau = 0.35

#Numerical for payg with new tax level
obj_kpayg = lambda kpayg: kpayg - (1/(1+(1+value_rho)/(2+value_rho)*(1-value_alpha)/value_alpha*value_tau)*(((1-value_alpha)*(1-value_tau))/(2+value_rho)*value_A*kpayg**value_alpha))

r_payg = optimize.root_scalar(obj_kpayg, bracket=[0.1, 100], method='brentq')

#Print results
print('For a tax value of 35% we observe in the payg-system: ')
print(f'k^*={r_payg.root:.4f}')

#numerical for ff with new tax level
obj_kff = lambda kff: kff - ((value_A*kff**value_alpha*(1-value_alpha))*(-value_tau+1/(value_rho+2)))

r_ff = optimize.root_scalar(obj_kff, bracket=[0.1, 100], method='brentq')

#Print results
print('And in the fully funded system we observe: ')
print(f'k^*={r_ff.root:.4f}')
```

    For a tax value of 35% we observe in the payg-system: 
    k^*=0.2089
    And in the fully funded system we observe: 
    k^*=0.1353



```python
# Assign a new tax value (tau)
value_tau = 0.05

#Numerical for payg with new tax level
obj_kpayg = lambda kpayg: kpayg - (1/(1+(1+value_rho)/(2+value_rho)*(1-value_alpha)/value_alpha*value_tau)*(((1-value_alpha)*(1-value_tau))/(2+value_rho)*value_A*kpayg**value_alpha))

r_payg = optimize.root_scalar(obj_kpayg, bracket=[0.1, 100], method='brentq')

#Print results
print('For a tax value of 5% we observe in the payg-system: ')
print(f'k^*={r_payg.root:.4f}')

#numerical for ff with new tax level
obj_kff = lambda kff: kff - ((value_A*kff**value_alpha*(1-value_alpha))*(-value_tau+1/(value_rho+2)))

r_ff = optimize.root_scalar(obj_kff, bracket=[0.1, 100], method='brentq')

#Print results
print('And in the fully funded system we observe: ')
print(f'k^*={r_ff.root:.4f}')
```

    For a tax value of 5% we observe in the payg-system: 
    k^*=0.5894
    And in the fully funded system we observe: 
    k^*=0.6197



```python
tau_vec = np.linspace(0,0.4,100)
obj_kpayg = lambda kpayg, value_tau: kpayg - (1/(1+(1+value_rho)/(2+value_rho)*(1-value_alpha)/value_alpha*value_tau)*(((1-value_alpha)*(1-value_tau))/(2+value_rho)*value_A*kpayg**value_alpha))

r_payg_vec = [optimize.root_scalar(obj_kpayg, bracket=[0.0001, 100], method='brentq',args=(tau_val)).root for tau_val in tau_vec ]



obj_kff = lambda kff, value_tau: kff - ((value_A*kff**value_alpha*(1-value_alpha))*(-value_tau+1/(value_rho+2)))

r_ff_vec = [optimize.root_scalar(obj_kff, bracket=[0.0001, 100], method='brentq',args=(tau_val)).root for tau_val in tau_vec]
```


```python
plt.plot(tau_vec,r_payg_vec,color='red')
plt.plot(tau_vec,r_ff_vec);
```

# Conclusion

From these results it would seem that capital levels in the fully funded system are actually higher given that the tax is low.
It would therefore seem that by increasing taxes the capital levels in both systems fall, but that the fully funded system falls at a greater rate than the payg-system.

In conclusion taxes seems to affect the fully funded system more severely than the payg-system.


```python

```

57817743a0c57bdd768772bee8dd9fafece666be
ipynb
Jupyter Notebook
modelproject/modelproject.ipynb
AskerNC/projects-2021-stressfaktoren
c2767869d98c8bfbe6d71ad5dd950e94f0878f6c
MIT
modelproject/modelproject.ipynb
AskerNC/projects-2021-stressfaktoren
c2767869d98c8bfbe6d71ad5dd950e94f0878f6c
MIT
modelproject/modelproject.ipynb
AskerNC/projects-2021-stressfaktoren
c2767869d98c8bfbe6d71ad5dd950e94f0878f6c
MIT
Qwen/Qwen-72B
1. YES
2. YES