# Model Project of the Optimal Number of Children for Household#

**Importing libraries**

Paste pip install git+https://github.com/elben10/pydst in the Terminal before you run the code, in order to have pydst Package.


```python
import numpy as np
import scipy as sp
from scipy import linalg
from scipy import optimize
from scipy.optimize import minimize_scalar
from scipy import interpolate
import sympy as sm
import math as ma

%matplotlib inline
import matplotlib.pyplot as plt
from matplotlib import cm
from mpl_toolkits.mplot3d import Axes3D

from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
import os
import pydst
import pandas as pd
import numpy as np
dst = pydst.Dst(lang='da')
```

## 1. Model description ##

In this project, we try to find the optimal number of children for household by combining utility function and budget constraint. In order to do so, we firstly study and visualize the utility function. The second, we find the general expression of the optimal number of children. Then we set up parameters to get the nummeral result of b and compare it with the result obtain from the general expression. Finally, we extend the model with the cost of having children.

**The model is:**

<div style="font-size: 16px">

$ max\,U(c, n)= {\gamma}\,\log\,n + (1-{\gamma})\,\log\,c $

$\text{s.t.}\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,\,c\,=\,(1-{\tau}{\phi}b)\,y\,-\,{\rho}b$

**where**

- c is the household consumption level
- n is the number of surviving children of the household
- ${\gamma}$ is the preference parameter between consumption and children. Thus, $0 < {\gamma} < 1$.
- y is the total income of the household
- ${\tau}$ is the share of income dedicated to grow children (per surviving child). Thus, $0 < {\tau} < 1$.
- ${\phi}$ is ratio of survival children to the total fertility. Thereby, $n\,=\,{\phi}\,b$. Thus, $0 < {\phi} < 1$.
- ${\rho}$ is cost of having one child born regardless of life expectancy.

Galor assumes ${\rho} = 0$. This leads to the conclusion that optimal number of the surviving children of the household is not affected by mortality. Following explanations and codes will prove the conclusion of the model by applying the algorithms from this course.

## 2. Analyze the utility function ##

###  a. Visualize the utility function in a 3D form### 

We will start with the visulaization of the utility function in the form of 3D. The exogenous variables of the funtion is ${n}$ and ${c}$. Each of these variables will be the x-axis and y-axis for the 3D graph, and z-axis will show the amount of utility. In order to show the effect which ${\gamma}$ has as an endogenous variable of the function, we will use Slider widget. By moving the dot on the scroll bar located upon the 3D graph, you can see how the utility function is changed by ${\gamma}$.


```python
# Defining all symbols:
c = sm.symbols('c')
n = sm.symbols('n')
y = sm.symbols('y')
b = sm.symbols('b')
gamma = sm.symbols('gamma')
tau = sm.symbols ('tau')
phi = sm.symbols ('phi')

# Setting up the basic for to make the slider of gamma:
def g(gamma):
    
# Defining function:
    u = gamma*sm.ln(n)+(1-gamma)*sm.ln(c)
    _u = sm.lambdify((n,c),u)

# Plotting the graph    
# a. grids
    n_vec = np.linspace(1,10,10)
    c_vec = np.linspace(1,100000,500)
    n_grid,c_grid = np.meshgrid(n_vec,c_vec,indexing='ij')
    f_grid = _u(n_grid,c_grid)

# b. main
    fig = plt.figure()
    ax = fig.add_subplot(1,1,1,projection='3d')
    cs = ax.plot_surface(n_grid,c_grid,f_grid,cmap=cm.jet)

# c. add labels
    ax.set_xlabel('$n$')
    ax.set_ylabel('$c$')
    ax.set_zlabel('$u$')

# d. invert and fix xaxis
    ax.invert_xaxis()
    plt.xlim([0,10])


# e. remove background
    ax.xaxis.pane.fill = False
    ax.yaxis.pane.fill = False
    ax.zaxis.pane.fill = False

# f. add colorbar
    fig.colorbar(cs);
    
# e. title and size
    plt.title("Utility function of the household \n with regard to the consumption \n and the number of surviving children", fontsize=15)
    plt.figure(figsize=(1000,1000))
    
# Finishing the codes for the slider:
slider = widgets.FloatSlider(
    value = 0.5,
    min = 0,
    max = 1,
    step = 0.01,
    description = '${\gamma}$ :',
    disabled = False,
    continuous_update = True,
    orientation = 'horizontal',
    readout = True,
    readout_format = 'd'
)
    

interact(g, gamma= slider);


```


    interactive(children=(FloatSlider(value=0.74, description='${\\gamma}$:', max=1.0, readout_format='d', step=0.…


###  b. Visualize the utility function in a 2D form### 

In order to better show how the utility changes with regard to the exogenous variables ${n}$, ${c}$, we will also plot the function on the 2D graph. Each of the graph lines shown on the graph below has the same utility. Also, by changing the value of gamma with slider, you can see the overall utility function. Please refer to the color bar to check how much utility is allocated to each of the graphs.


```python
# Setting up the basic for to make the slider of gamma:
def g(gamma):
# Defining function:    
    u = gamma*sm.ln(n)+(1-gamma)*sm.ln(c)
    _u = sm.lambdify((n,c),u)
#Plotting the graph
# a. grids
    n_vec = np.linspace(1,10000,10)
    c_vec = np.linspace(1,1000000,500)
    n_grid,c_grid = np.meshgrid(n_vec,c_vec,indexing='ij')
    f_grid = _u(n_grid,c_grid)

# b. Adding info on axis
    fig = plt.figure()
    ax = fig.add_subplot(1,1,1)
    levels = [0,1,2,3,4,5,6,7,8,9,10,11,12]
    
# c. Main
    cs = ax.contour(n_grid,c_grid,f_grid,levels=levels,cmap=cm.jet)
    
# d. add labels
    ax.set_xlabel('$n$')
    ax.set_ylabel('$c$')

# e. add color bar
    clb = fig.colorbar(cs)
    clb.ax.set_title('u')
    
# f. title and size
    plt.title("Utility function of the household \n with regard to the consumption \n and the number of surviving children", fontsize=15)
    plt.figure(figsize=(1000,1000))

# Finishing the codes for the slider:
    slider = widgets.FloatSlider(
    value = 0.5,
    min = 0,
    max = 1,
    step = 0.01,
    description = '${\gamma}$ :',
    disabled = False,
    continuous_update = True,
    orientation = 'horizontal',
    readout = True,
    readout_format = 'd'
)
    

interact(g, gamma= slider);
```


    interactive(children=(FloatSlider(value=0.74, description='${\\gamma}$:', max=1.0, readout_format='d', step=0.…


### c. Numerical Calculations of the utility function with the given ${\gamma}$ values:  ###

In order to show how the changes of the number of children affect the overall utility for a given level of consumption, we print multiple evaluations of the utility function.


```python
def u_func(n,c,gamma = 0.50): #define the utility function with gamma = 0.5 
    return gamma * np.log (n) + gamma * np.log (c)

n_list = [1,2,3,4,5,7,8,9,10] # The number of children ranges from 1 to 10.
c = 10000 # The consumption level is 10000

for n in n_list: # loop through each element in n_list
    u = u_func(n, c, gamma = 0.5)
    print(f'n = {n:.3f}, c = {c:.3f} -> u = {u:.3f}')
```

    n = 1.000, c = 10000.000 -> u = 4.605
    n = 2.000, c = 10000.000 -> u = 4.952
    n = 3.000, c = 10000.000 -> u = 5.154
    n = 4.000, c = 10000.000 -> u = 5.298
    n = 5.000, c = 10000.000 -> u = 5.410
    n = 7.000, c = 10000.000 -> u = 5.578
    n = 8.000, c = 10000.000 -> u = 5.645
    n = 9.000, c = 10000.000 -> u = 5.704
    n = 10.000, c = 10000.000 -> u = 5.756


## 3. Solve the utility function subject to the budget constraint ##

### a. Solve the utility function with the budget constraint to get the optimal value of b ###

In this section we maximize the utility function under the budget constraint. By substituting exogenous variables of the utility function regarding the variable b, we will get the optimal value of b under the maximized utility. This will lead to the constitution of the basis of the conclusion of the theory.


```python
# define all symbols
c = sm.symbols('c')
n = sm.symbols('n')
y = sm.symbols('y')
b = sm.symbols('b')
gamma = sm.symbols('gamma')
tau = sm.symbols ('tau')
phi = sm.symbols ('phi')

# define the utility function
objective = (1-gamma) * sm.ln(c) + gamma * sm.ln(phi*b)

# define the budget constraint
budget_constraint = sm.Eq((c/(1-(tau*phi*b))), y)

# islate c in the budget constraint
c_from_con = sm.solve(budget_constraint, c)
c_from_con

# substitute c as a function of b 
objective_subs = objective.subs(c, c_from_con[0])
objective_subs

# take the derivative of the function from above 
foc = sm.diff(objective_subs, b)

# solve the FOC 
solve_this = sm.Eq(foc,0)

# find the general expression of the optimal value of b
final_b = sm.solve(solve_this,b)
print(objective_subs)
print(final_b)
```

    gamma*log(b*phi) + (-gamma + 1)*log(y*(-b*phi*tau + 1))
    [gamma/(phi*tau)]


### b. Analysis on the conclusion of the simplified model of Galor ### 

As it is shown above, at the maximized utility level, ${b}$ $=$ $\frac {\gamma} {{\phi}{\tau}}$. In other words, ${n} = {\phi}{b} =\frac {\gamma} {\tau}$. So as the model is saying, the optimal number of the surviving children is not affected by the mortality rate. In order to see how exactly the ${n}$ can be changed as the endogenous variable changes, we will take a look at the graph below. For better understanding, the graph is 2D with ${\tau}$ as x-axis and ${n}$ as y-axis. By changing ${\gamma}$ with the slider, you can see how the overall value of ${n}$ can be changed.  


```python
dst = pydst.Dst(lang='da')

# Set up for gamma slider
def u(gamma):
    
#main
    x = np.linspace(0.1,1,100)
    u = gamma/x
    fig = plt.figure()
    plt.plot(x, u)

# Label work        
    plt.ylim([0,5])
    ax = fig.add_subplot(1,1,1)
    ax.set_xlabel('$tau$')
    ax.set_ylabel('$n$')

# title and size
    plt.title("Optimal number of the surviving children \n with regard to gamma and tau ", fontsize=15)
    plt.figure(figsize=(1000,1000))

# Slider code sum up
slider = widgets.FloatSlider(
    value = 0.5,
    min = 0,
    max = 1,
    step = 0.01,
    description = '${\gamma}$:',
    disabled = False,
    continuous_update = True,
    orientation = 'horizontal',
    readout = True,
    readout_format = 'd'
)
    

interact(u, gamma = slider);
```


    interactive(children=(FloatSlider(value=0.5, description='${\\gamma}$:', max=1.0, readout_format='d', step=0.0…


As the ${\gamma}$ increases the overall graph shifts upwards, meaning ${n}$ increases. Also, as ${\tau}$ increases, ${n}$ decreases. The reason is that ${\gamma}$ indicates the preference of having children over consumption, while ${\tau}$ shows the share of the income dedicated to growing up a child.

### c. Find the optimal value of b for given value of parameters using the optimize method directly ###

We give some numerical values to the parameters in order to prove the process done above is numerically correct.


```python
# redefine the utility function without c
utility_func = sm.lambdify((b, c, gamma, phi, tau, y), objective_subs)

# set up numeral value for the parameters
gamma = 0.8
tau = 0.2
phi = 0.8 
y = 10000

#grids
c_vec = np.empty(1000)
b_vec = np.linspace(1, 10, 10)

#solve for each b in grid
for b in enumerate(b_vec):
    obj = lambda b: -utility_func(b, c, gamma, phi, tau, y)
    b0 = 1

# find the optimal value of b
result = optimize.minimize_scalar(obj,b0, method='bounded',bounds=[0,10])
   
result
```




         fun: -2.6292159808041085
     message: 'Solution found.'
        nfev: 12
      status: 0
     success: True
           x: 4.999998342740495




```python
# make use of the general expression b to calculate b subject to the set up of the parameters

# repeat the set-up of the parameters
gamma = 0.8
tau = 0.2
phi = 0.8 

# calculate b using a function of gamma, phi and tau
def b_func(gamma, phi, tau):
    b_optimal = gamma/(phi*tau)
    return b_optimal

b_func(gamma, tau, phi)
```




    4.999999999999999



*The value of the b is the same when it is obtained from the numerical optimization and when it is from the equation on b generated from the utility maximization. This indicates the overall analysis on the model above is less likely to be wrong.*

## 4. Extension of the model ##

For the extension of the model we will consider not only the opportunity cost of having children, but also, the cost of giving a child birth. Galor assumes ${\rho} = 0$. In the extension we relax Galor's assumption as ${\rho} > 0$.


```python
# define all symbols
c = sm.symbols('c')
n = sm.symbols('n')
y = sm.symbols('y')
b = sm.symbols('b')
gamma = sm.symbols('gamma')
tau = sm.symbols('tau')
phi = sm.symbols('phi')
rho = sm.symbols('rho')

# define the utility function
objective = (1-gamma) * sm.ln(c) + gamma * sm.ln(phi*b)

# define the budget constraint
budget_constraint = sm.Eq(((c+(rho*b))/(1-(tau*phi*b))), y)

# islate c in the budget constraint
c_from_con = sm.solve(budget_constraint, c)
c_from_con

# substitute c as a function of b 
objective_subs = objective.subs(c, c_from_con[0])
objective_subs

# take the derivative of the function from above 
foc = sm.diff(objective_subs, b)

# solve the FOC 
solve_this= sm.Eq(foc,0)

# find the general expression of the optimal value of b
final_b= sm.solve(solve_this,b)
print(objective_subs)
print(final_b)
```

    gamma*log(b*phi) + (-gamma + 1)*log(-b*phi*tau*y - b*rho + y)
    [gamma*y/(phi*tau*y + rho)]


### Find the optimal value of b for given value of parameters considering the consumption of children ###


```python
gamma=0.8
tau=0.2
phi=0.8 
rho=400
y=10000

# calculate b using a function of gamma, phi and tau
def b_func(gamma, phi, tau,rho):
    b_optimal = gamma*y/(phi*tau*y+rho)
    return b_optimal

b_func(gamma, tau, phi,rho)
```




    3.9999999999999996



## 5. Conclusions ##

The optimal number of children per household is independent to the household income and mortality index, under the simplified model of Galor. In this case the number surviving children per household is rather affected by the share of income dedicated to growing up a child, and the preference on the children. The extension model contains the cost of giving a child a birth as a positive value, so this leads to income level and mortality index affecting the decision on fertility.

ee321c9c970f53a022c16fa70c69d91ac25034d2
ipynb
Jupyter Notebook
modelproject/model project_final.ipynb
NumEconCopenhagen/projects-2019-leg
afdc76c770f29d1affe819db856d742e1b1ececa
MIT
modelproject/model project_final.ipynb
NumEconCopenhagen/projects-2019-leg
afdc76c770f29d1affe819db856d742e1b1ececa
MIT
2019-04-09T12:00:18.000Z
2019-05-23T12:50:23.000Z
modelproject/model project_final.ipynb
NumEconCopenhagen/projects-2019-leg
afdc76c770f29d1affe819db856d742e1b1ececa
MIT
2019-05-14T13:37:01.000Z
2019-05-14T13:37:01.000Z
Qwen/Qwen-72B
1. YES
2. YES