# Overlapping Generation Model with Pay-As-You-Go

Imports and set magics:


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

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

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


## Model description

We consider a production economy populated by agents that live for two periods and as such overlapping generations. Time, $t$, is discrete and infinite. $L_t$ individuals are born i period t and the popultation grows at a constant rate $n$ such that:

$$L_t=L_{t-1}+n)$$

Agents derive utility from consumption while alive:

 $$U_t=u(c_{1t})+\frac{1}{1+\rho}u(c_{2t})$$

where
1. $c_{1t}$ is consumption when young at time $t$
* $c_{2t}$ is consumption when old at time $t$+1
* $\rho$ is the dicount rate and
* $u(x)=log(x)$, i.e. the special case of logagrithmic utility is examined

Agents work when young supplying 1 unit of labour inelastically at wage rate $w_t$. They split the labour income between consumption, $c_{1t}$, and savings, $s_t$. When old, they retire and solely consume their gross savings. We let $r_{t+1}$ denote the interest rate between $t$ and $t+1$. Thus, agents face the following budget constrints in period 1 and 2 respectively:

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

By substituting out $s_t$ the lifetime budget constraint emerges:

$$c_{1t}+\frac{c_{2t+1}}{1+r_{t+1}}=w_t$$

We assume the production function is Cobb-Douglas with standard properties so:

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

## Analytical solution

Using Sympy the analytical solution can be found using a similar approach as presented in the KU course Macroeconomics III.
First we consider the households' problem.

### Households' problem

The households' problem is solved by firstly setting up the lagrangian, finding the first order conditions with respect to consumptions from which the Euler equation can be derived.


```python
#Define variabels
c1 = sm.symbols('c_1t')
c2 = sm.symbols('c_2t+1')
st = sm.symbols('s_t')
rho = sm.symbols('rho')
wt = sm.symbols('w_t')
rt1 = sm.symbols('r_t+1')
la = sm.symbols('lambda')
```


```python
#Set up utility, budget constraints and lagrangian
uc1 = sm.ln(c1)
uc2 = sm.ln(c2)
U = uc1 + 1/(1+rho) * uc2
bc1 = wt-st #budget constraint for period 1
bc2 = (1+rt1)*st ##budget constraint for period 1
bc = c1 + c2 / (1+rt1) - wt #lifetime budget constraint
L = U + la * (-1)*bc

#Finding FOCs
foc1 = sm.Eq(0,sm.diff(L,c1))
foc2 = sm.Eq(0,sm.diff(L,c2))

#Solve for the lagrangian multiplyer
la1 = sm.solve(foc1,la)[0]
la2 = sm.solve(foc2,la)[0]
                
#Euler equation
euler = sm.solve(sm.Eq(la1,la2),c1)[0] #Equate lambdas and solve
eulereq = sm.Eq(euler, c1)
display(eulereq)
```


$\displaystyle \frac{c_{2t+1} \left(\rho + 1\right)}{r_{t+1} + 1} = c_{1t}$


Thus the Euler equation has been found. This equation decribes how marginal utilities of consumption between youth and old age are optimally related.
To characterise optimal savings, we combine the Euler equation and agents' budget constraints for individual periods.


```python
#substitute for savings
euler_sav = eulereq.subs(c1,bc1)
euler_sav1 = euler_sav.subs(c2,bc2)

#Solve for optimal savings

saving = sm.solve(euler_sav1,st)[0] #Solve
opt_saving = sm.Eq(saving, st) #equate with s_t
display(opt_saving)
```


$\displaystyle \frac{w_{t}}{\rho + 2} = s_{t}$


Thus optimal savings have been found. For the log utility case substitution and income effects of the interest rate cancel out. A constant fraction of labour income is saved. 

### Firms' Problem

The firms' problem for the Cobb-Douglas case is solved using first order conditions. The firms maximisation problem is
$$ max      K_{t}^{\alpha}L_{t}^{1-\alpha}-r_tK_t-w_tL_t$$


```python
#Define variables
alpha = sm.symbols('alpha')
K = sm.symbols('K_t')
k = sm.symbols('k_t')
k1 = sm.symbols('k_t+1')
L = sm.symbols('L_t')
n = sm.symbols('n')
rt = sm.symbols('r_t')
wt = sm.symbols('w_t')
Y = sm.symbols('Y')

#Set up production function
Y = K**alpha * L**(1-alpha)

#Find partial derivatives
dyk = sm.diff(Y,K)
dyl = sm.diff(Y,L)

#Solve for interest rate and wage
rate = sm.Eq(rt,dyk) 
wage = sm.Eq(wt,dyl)

#Find per worker expressions
rate_pw = sm.Eq(rt,(alpha*k**(alpha-1)))
wage_pw = sm.Eq(wt, ((1-alpha)*k**(alpha)))
display(rate_pw,wage_pw)
```


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



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


Thus the firms' problem has been solved for per worker optimal values.

### Capital Accumulation per Worker

Capital per worker in period $t+1$ is equal to the savings in the previous period:
$$k_{t+1}(1+n)=s_t$$
Optimal savings are substituted in and the equation is solved for the per worker capital in period $t+1$.


```python
#Substitute for optimal savings
savings = opt_saving.subs(wt, ((1-alpha)*k**(alpha)))

#Use definition of aggregate capital
cap = k1 * (1+n)
agg_cap = savings.subs(st,cap)

#Solve for k_t+1 to obtain law of-motion
lom = sm.solve(agg_cap,k1)[0]
lom = sm.Eq(k1,lom)
display(lom)
```


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


Thus the law of motion for capital accumulation per worker has been found.The function is concave, so there will be a unique steady state (SS) and capital stock will converge to it
To find the steady state, we use that in SS $k_t=k_{t+1}=k^*$:


```python
#Define SS capital variable 
kstar = sm.symbols('k^*')

#Use rule for SS
lom1 = lom.subs(k1, kstar)
lom2 = lom1.subs(k, kstar)

#Solve for k*
lom3 = sm.solve(lom2, kstar)[0]
lom4 = sm.Eq(kstar,lom3)
display(lom4)
```


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


Thus the algebraic SS solution has been found.
For later use, we turn the solution into a python function


```python
ss_func = sm.lambdify((n, rho, alpha),lom3)
```

## Numerical solution

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

$$0 = \frac{(k_t^*)^\alpha (1-\alpha)}{n\rho+2n+\rho+2}-k_{t+1}^*$$

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


```python
#Define parameter values
alpha = 1/3
n = 0.01
rho = 0.1

#Define objective function i.e. the low of motion as seen above
f = lambda k:k**alpha #input for the objective function
obj_kss = lambda kss: kss - (f(kss)*(1-alpha))/(n*rho+2*n+rho+2) #the objective functiin is the law of motion

#Optimize using the bisect method
result = optimize.root_scalar(obj_kss,bracket=[0.1,100],method='bisect')

print(f'The analytical SS solution is: {ss_func(n, rho, alpha):.3f}')
print('The numerical SS soltion is:',result.root)
```

    The analytical SS solution is: 0.176
    The numerical SS soltion is: 0.1762187469065779


The analytical and numarical results are virtually the same.


```python
def transition(k,alpha,rho,n):
    w = k**alpha * (1-alpha)
    s = w / (rho+2)
    k1 = s / (n+1)
    return k1

def transition_curve(alpha,rho,n):
    
    #grids
    k_1 = np.linspace(1e-20,6,1000)
    k_2 = np.empty(1000)
    
    #construct array of next periods capita
    for i,k in enumerate(k_1):
        k_plus = transition(k,alpha,rho,n)
        k_2[i] = k_plus
        
    return k_1, k_2

def fig(rho):
    k_1, k_2 = transition_curve(alpha,rho,n)
    fig = plt.figure()
    ax = fig.add_subplot(1,1,1)
    ax.plot(k_1,k_2)
    ax.plot(k_1,k_1, '--', color='grey')
    plt.text(0.35,0.26,'$k_t+1$')
    plt.text(0.35,0.32,'$k_t=k_t+1$')
    plt.text(0.2,0.22,'$k^*$')
    ax.set_xlabel('$k_t$')
    ax.set_ylabel('$k_t+1$')
    ax.set_title('Transition curve')
    ax.set_xlim([0,0.4])
    ax.set_ylim([0,0.4]);
    return

import ipywidgets as widgets
widgets.interact(fig, 
    rho = widgets.FloatSlider(description='rho', min=0, max=16, step=0.01, value=0.1),
);
```


    interactive(children=(FloatSlider(value=0.1, description='rho', max=16.0, step=0.01), Output()), _dom_classes=…


## Pay-As-You-Go Extenstion

In a pay-as-you-go (PAYG) system the governemnt raises contributions $d_t$ from current young and pays them out as benefit to current old:
$$b_t=(1+n)d_t$$
Thus the budget contraints change from the basic example. When young in period $t$ th ebudget constraint is
$$c_{1r}+s_t+d_t=w_t$$
while when old in period $t+1$ the budget constraint is
$$c_{2t+1}=(1+r_{t+1})s_t+(1+n)d_{t+1}$$
Thus the lifetime budget constraint is 
$$c_{2t+1}=(1+r_{t+1})s_t+dw_t$$

We start by setting up a modified lagrangian, solving first order conditions and solving for the new Euler equation.


```python
#Define new variables
d = sm.symbols('d_t')

#New life time budget constraint
new_bc = c1+c2/(1+rt1)-d/(1+rt1)*wt-(1-d)*wt

#New lagrangian
L1 = U + la * (-1)*new_bc

#First order conditions
foc3 = sm.Eq(0,sm.diff(L1,c1))
foc4 = sm.Eq(0,sm.diff(L1,c2))

#Solve first order conditions
lam1 = sm.solve(foc3,la)[0]
lam2 = sm.solve(foc4,la)[0]

#Solve for Euler equation
euler1 = sm.solve(sm.Eq(lam1,lam2),c1)[0] 
euler2 = sm.Eq(euler1, c1)
display(euler2)
```


$\displaystyle \frac{c_{2t+1} \left(\rho + 1\right)}{r_{t+1} + 1} = c_{1t}$


It is immidiately noticable that the Euler equation has not changed. This is due to contributions are taken as a given.
Again, we substitute out consumption with budget constraints and solve for savings.


```python
# substituting
new_euler3 = euler2.subs(c2,((1+rt1)*st+d*wt))
new_euler4 = new_euler3.subs(c1,(wt-st-d*wt))

# solving for savings
new_saving = sm.solve(new_euler4,st)[0]
opt_saving = sm.Eq(new_saving,st)
opt_saving
```




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



A new optimal individual saving has been found.
While the firms' problem remains unchanged we immidiately proceed to finding the new capital accumulation.
We proceed by creating and plotting the new capital accumulation curve.

It should be noticed that this does not exhibit the same standard properties at would be expected from analysing this scenario, which basis for believing that this might be incorrect. Specifically, individual savings are not positively correlated with the the future benefits.


```python
sav1 = opt_saving.subs(rt1,(alpha*k**(alpha-1)))
sav2 = sav1.subs(wt,(1-alpha)*k**(alpha))
sav_func = sm.lambdify((k,alpha,d,rho),new_saving)

def transition_payg(k,alpha,rho,d,n):
    w = k**alpha * (1-alpha)
    s = sav_func(k,alpha,d,rho)
    k1 = s / (n+1)
    return k1

#Something goes wrong here. This will be fixed as quickly as possible.
```


```python

#initializing parameters again
d = 0.3
alpha = 0.3
n = 0.01
rho = 0.1

def fig_payg(d):
    k_1_payg, k_2_payg = transition_curve_payg(alpha,rho,n,d,rt)
    k_1, k_2 = transition_curve(alpha,rho,n)
    fig = plt.figure()
    ax = fig.add_subplot(1,1,1)
    ax.plot(k_1_payg,k_2_payg,label="PAYG")
    ax.plot(k_1,k_2,label="Baseline")
    ax.plot(k_1_payg,k_1_payg, '--', color='grey')
    ax.set_xlabel('$k_t$')
    ax.set_ylabel('$k_t+1$')
    ax.set_title('Transition curve')
    ax.set_xlim([0,0.2])
    ax.set_ylim([0,0.2]);
    ax.legend()
    return

import ipywidgets as widgets
widgets.interact(fig_payg, 
    d = widgets.FloatSlider(description='d', min=0, max=1, step=0.01, value=0.1),
);
```


    interactive(children=(FloatSlider(value=0.1, description='d', max=1.0, step=0.01), Output()), _dom_classes=('w…


Obviously this plot is not working, due to the above root cause. The plot when functioning would show that for logarithmic utility functions a PAYG system will decrease capital accumulation.

# Conclusion

It can be seen that solving the PAYG OLG model both analytically and numerically provides the same useful insights into the funding of a social security system where the young pays taxes which are provided to the old in terms of benefits. Capital accumulation is smaller in an economy with PAYG.

ce188ce8643cc3999c4dc17b29751b78f878947c
ipynb
Jupyter Notebook
modelproject/OLG model with PAYG.ipynb
NumEconCopenhagen/projects-2020-jrf330
f15c4e0ff4686d44c9a2b1fd293e4f669ef25bef
MIT
modelproject/OLG model with PAYG.ipynb
NumEconCopenhagen/projects-2020-jrf330
f15c4e0ff4686d44c9a2b1fd293e4f669ef25bef
MIT
2020-04-17T13:19:18.000Z
2020-05-11T12:23:58.000Z
modelproject/OLG model with PAYG.ipynb
NumEconCopenhagen/projects-2020-jrf330
f15c4e0ff4686d44c9a2b1fd293e4f669ef25bef
MIT
Qwen/Qwen-72B
1. YES
2. YES