Fertility¶
AIM: match the population data with the fertility rate, to be able to extrapolate for next years
The fertility rate ($fr$) is the average number of children per women, and has to be above 2 to keep the population stable.
ISTAT computes the fertility rate in DCIS_FECONDITA1, but I have to do some data reconcilation to make sure that I have a function I can use to compute newborns from the population data.
In [1]:
import numpy as np
import pandas as pd
import json
import matplotlib
import matplotlib.pyplot as plt
import plotly.express as px
import plotly.graph_objects as go
import plotly.io as pio
from plotly.subplots import make_subplots
import wbgapi as wb # Italian data are probably copied from ISTAT anyway
from copy import deepcopy
import sdmx
import warnings
pio.renderers.default = 'vscode+notebook'
pd.options.plotting.backend = "plotly"
warnings.filterwarnings("ignore", category=FutureWarning)
def get_colors(n, cmap_name="rainbow"):
"""Get colors for px colors_discrete argument, given the number of colors needed, n."""
cmap = matplotlib.colormaps[cmap_name]
colors = [cmap(i) for i in np.linspace(0, 1, n)] # Generate colors
colors_str = [f"rgba({int(color[0]*250)}, {int(color[1]*250)}, {int(color[2]*250)}, 1.0)" for color in colors]
return colors_str
In [2]:
dfp=pd.read_csv("../data/pop_by_age_year.csv", index_col=0).rename(columns=int) # From Notebook#14
In [ ]:
df6 = sdmx.to_pandas(sdmx.Client("ISTAT").data(
resource_id="25_326_DF_DCIS_FECONDITA1_5",
key={
"REF_AREA": "IT" # entire Italy
}
)).reset_index().astype({"TIME_PERIOD": int})
df6
Out[ ]:
In [4]:
px.line(
df6,
x="TIME_PERIOD",
y="value",
color="CITIZENSHIP",
).update_layout(
title="Italian fertility rate by citizenship",
width=1000,
).show()
In [ ]:
df7 = sdmx.to_pandas(sdmx.Client("ISTAT").data(
resource_id="25_326_DF_DCIS_FECONDITA1_7",
key={
"REF_AREA": "IT" # entire Italy
}
)).reset_index()
df7
Out[ ]:
In [6]:
# clean the dataset
fert_by_age = (
df7
.query("CITIZENSHIP == 'TOTAL'")
.assign(year = pd.to_datetime(df7["TIME_PERIOD"]).dt.year)
.assign(age = df7["MOTHER_AGE"].str.extract(r"(\d+)").astype(int))
.pivot(index="age", columns="year", values="value")
)
fert_by_age
Out[6]:
In [7]:
fig = px.line(
data_frame=fert_by_age,
x=fert_by_age.index,
y=fert_by_age.columns,
color_discrete_sequence=get_colors(n=fert_by_age.columns.size),
labels={"x": "Year", "y": "Fertility rate"},
).update_layout(
xaxis_title="Age",
yaxis_title="Rel. fertility by age",
yaxis_range=[0, 100],
legend_title="Scenarios",
margin=dict(l=10, r=10, t=20, b=10),
width=780,
height=320,
)
fig.write_html("../images_output/fertility_by_age.html", auto_play=False)
fig.show()
In [8]:
# Get a cumulative fertility rate by age
fert_by_age_cum = fert_by_age.div(fert_by_age.sum(axis=0), axis=1) # first normalize the fertility rate by age
fert_by_age_cum = fert_by_age_cum.cumsum() # then get the cumulative fertility rate by age
fert_by_age_cum.loc[17] = 0 # enforcing 0 at age 17 to set constraints on sigmoid
fert_by_age_cum.sort_index(inplace=True)
px.line(
data_frame=fert_by_age_cum,
x=fert_by_age_cum.index,
y=fert_by_age_cum.columns,
color_discrete_sequence=get_colors(n=fert_by_age_cum.columns.size),
labels={"x": "Year", "y": "Fertility rate"},
).update_layout(
width=1000
).show()
In [9]:
# fit the sigmoids and extrapolate for next years
from scipy.optimize import curve_fit
def constrained_sigmoid(x: float, x0: float, k: float) -> float:
"""
Compute a constrained sigmoid function that is 0 at the start and 1 at the end.
Parameters:
x (float): The input value.
x0 (float): The x-value of the sigmoid's midpoint.
k (float): The steepness of the curve.
Returns:
float: The output of the constrained sigmoid function.
"""
y = 1 / (1 + np.exp(-k * (x - x0)))
y = y - y[0]
y = y / y[-1]
return y
def fit_sigmoid(x, y):
p0 = [np.mean(x), 1]
popt, pcov = curve_fit(constrained_sigmoid, x, y, p0, method='dogbox')
return popt
fig, axs = plt.subplots(5,6, figsize=[15, 10])
coeff_by_year = {} # {year: [x0, k]}
for iplot, year in enumerate(fert_by_age_cum.columns):
x = fert_by_age_cum.index
y = fert_by_age_cum.loc[:, year]
coeff_by_year[year] = fit_sigmoid(x, y)
y_fit = constrained_sigmoid(x, *coeff_by_year[year])
ax = axs.flatten()[iplot]
ax.scatter(x, y, color="black", marker="o", s=2)
ax.set_title(f"Year: {year}")
ax.plot(x,y_fit, color="red", lw=0.5)
ax.grid()
plt.tight_layout()
plt.show()
In [10]:
df_coeff = pd.DataFrame(coeff_by_year).T
df_coeff.columns = ["x0", "k"]
fig = make_subplots(rows=1, cols=2).update_layout(
width=1000,
title="Sigmoid parameters over time: x0 (center) an k (steepness)"
)
fig.add_trace(go.Scatter(x=df_coeff.index, y=df_coeff["x0"], mode="markers", name="x0"),row=1, col=1)
fig.add_trace(go.Scatter(x=df_coeff.index, y=df_coeff["k"], mode="markers", name="k"), row=1, col=2)
# I will keep x0 linear and k constant to have just one parameter to control (x0, which is average age of fertility)
x0_coeff = np.polyfit(df_coeff.index, df_coeff["x0"], 1)
k = df_coeff["k"].mean()
fig.add_trace(go.Scatter(x=df_coeff.index, y=np.polyval(x0_coeff, df_coeff.index), mode="lines", name="x0 linear fit"), row=1, col=1)
fig.add_trace(go.Scatter(x=df_coeff.index, y=[k]*df_coeff.index.size, mode="lines", name="k mean"), row=1, col=2)
fig.show()
print("x0_coeff", x0_coeff)
print("k", k)
Visually, it looks ok to model the cumulative sum of fertility over mothers' age like this: a sigmoid, function of x0 and k, where:
x0is a linear function of the year of observation (i.e., motherhood age is shifing ahead with years)kis the steepness of the sigmoid, it is going up and down, but there is not a significant error if I keep it constant (otherwise I should guess whether to fit it with a parabola or a sinusoid!)
Now that I know the disribution of the fertility rate, I can distribute the number of new borns accordingly.
In [11]:
total_newborns_by_year = dfp.iloc[0,:].shift(-1).dropna().to_dict() # newborns @ year of observation (shifting by -1 year beacause if there is a <1 year old on January 1st, he was born in the previous year)
fertiliy_by_year = {}
for year, tot_nb in total_newborns_by_year.items():
women = dfp.loc[17:50, year]/2 # assuming half of the population are women, which is reasonable globally but may not be the case for each age group 18-50
x0 = np.polyval([ 0.09458, -159.9], year)
k = -0.296
sigmoid = constrained_sigmoid(women.index, x0, k)
distrib = np.diff(sigmoid) # turning the sigmoid to a normaized distribution from 18 to 50
mothers = sum(distrib * women.values[1:])
fertiliy_by_year[year] = tot_nb / mothers
# REMEMBER: the fertility_by_year I computed here will be very similar to the official one by ISTAT,
# but slightly different because there may be diffences in the calculations (e.g., exact age of the mothers, infancy mortality, etc.).
# However, it is consistent to generate the correct number of newborns.
Store and recover data¶
In [12]:
# Store the coefficients
fertility_json = {
"last_observed_year": max(fertiliy_by_year.keys()),
"last_observed_value": fertiliy_by_year[max(fertiliy_by_year.keys())],
"sigmoid_x0": list(x0_coeff),
"sigmoid_k": k,
"assumed_women_ratio": 0.5
}
with open("../data/fertility_fitted.json", "w") as f:
json.dump(fertility_json, f, indent=4)
In [13]:
# make a function that recovers the data and gives the nexborns from the population
fdata = json.load(open("../data/fertility_fitted.json"))
print("Fertility coefficients:")
display(fdata)
def get_newborns_next_year(dfp, year, fertility_rate, fdata):
"""Fom population and fertility rate for a given year,
get the number of newborns for the NEXT year (fist January of year+1).
"""
def constrained_sigmoid(x0, k):
"""Get sigmoid function, which is a cumsum of a distribution."""
x = np.arange(17, 50+1)
y = 1 / (1 + np.exp(-k * (x - x0)))
y = y - y[0]
y = y / y[-1]
return y
pop = dfp[year]
mothers_potential = pop.loc[18:50].to_numpy() * fdata["assumed_women_ratio"]
x0 = np.polyval(fdata["sigmoid_x0"], year)
distrib = np.diff(constrained_sigmoid(x0, fdata["sigmoid_k"]))
newborns = sum(mothers_potential * distrib) * fertility_rate
return int(newborns)
get_newborns_next_year(dfp=dfp, year=2013, fertility_rate=1, fdata=fdata)
Out[13]:
Make plots compared to official data¶
In [14]:
# Get also the values from the World Bank that date back to 1960 (but should be the same as ISTAT)
ss_wb_fr = wb.data.DataFrame("SP.DYN.TFRT.IN", "ITA").T.dropna()["ITA"]
ss_wb_fr.index = [ int(year[-4:]) for year in ss_wb_fr.index ]
ss_wb_fr
Out[14]:
In [15]:
# compare all the fertility rates
fig = go.Figure().add_trace(go.Scatter(
x=df6.query("CITIZENSHIP == 'TOTAL'").TIME_PERIOD,
y=df6.query("CITIZENSHIP == 'TOTAL'").value,
mode="lines",
name="From ISTAT 'TOTAL'",
)).add_trace(go.Scatter(
x=fert_by_age.columns.astype(int),
y=fert_by_age.sum(axis=0)/1000,
mode="lines",
name="From ISTAT age-by-age summed",
)).add_trace(go.Scatter(
x=ss_wb_fr.index,
y=ss_wb_fr.values,
mode="lines",
name="From World Bank",
)).add_trace(go.Scatter(
x=list(fertiliy_by_year.keys()),
y=list(fertiliy_by_year.values()),
mode="lines",
name="From calculations",
)).update_layout(
width=1000,
title="Comparison of Fertility Rates from different sources and calculations",
xaxis_title="Year",
yaxis_title="Fertility Rate"
)
fig.show()
Note that there is still some discrepancy between the fertility rate I compute and the ISTAT one, that may be due to:
- approximating exactly 50% of the population as women
- approximation in the sigmoid model shape with no mothers below 18 and above 50 years old
- the indirect way I'm computing the
mothersfrom a distribution fitted on fertility rate
In [16]:
# https://esploradati.istat.it/databrowser/#/it/dw/categories/IT1,POP,1.0/POP_DEMOPROJ/DCIS_PREVDEM1/IT1,165_889_DF_DCIS_PREVDEM1_3,1.0
df3 = sdmx.to_pandas(sdmx.Client("ISTAT").data(
resource_id="165_889_DF_DCIS_PREVDEM1_3",
key={
"REF_AREA": "IT", # entire Italy
"DATA_TYPE": "TFR" # Total Fertility Rate
}
)).reset_index().astype({"TIME_PERIOD": int})
df3
Out[16]:
In [17]:
# pivot data to have year, and OBS_VALUE of INT_PREV="PROJMED" PROJLOW80 and PROJUPP80
df3p = df3.pivot(index="TIME_PERIOD", columns="FORECAST_INTERVAL", values="value")
df3p.head()
Out[17]:
In [18]:
# plot the data with INT_PREV="PROJMED" as line and INT_PREV="PROJLOW80" and "PROJUP80" as area of the confidence interval
go.Figure(
).add_trace(go.Scatter(
x=ss_wb_fr.index,
y=ss_wb_fr.values,
mode="lines",
name="From World Bank",
)
).add_trace(go.Scatter(
x=df3p.index,
y=df3p["PROJMED"],
mode="lines",
name="Projection",
line=dict(color="black")
)
).add_trace(go.Scatter(
x=df3p.index,
y=df3p["PROJLOW80"],
mode="lines",
name="Projection 80% confidence (lower bound)",
fill=None,
line=dict(color="gray"), # Change color and make it more transparent
)
).add_trace(go.Scatter(
x=df3p.index,
y=df3p["PROJUPP80"],
mode="lines",
name="Projection 80% confidence (upper bound)",
fill="tonexty",
line=dict(color="gray"), # Change color and make it more transparent
)
).update_layout(
width=1000,
title="Projection of Fertility Rate from ISTAT",
xaxis_title="Year",
yaxis_title="Fertility Rate",
).show()
Assuming 3 scenarios for future fertility rate¶
In [19]:
last_year = max(fertiliy_by_year.keys())
end_year = 2100
fr_projections = {}
fr_projections['fertility: ↓ 0 in 2100'] = pd.Series(np.linspace(fertiliy_by_year[last_year], 0, end_year-last_year+1), index=range(last_year, end_year+1))
fr_projections['fertility: eq. to last'] = pd.Series([fertiliy_by_year[last_year]] * (end_year-last_year+1), index=range(last_year, end_year+1))
fr_projections['fertility: ↑ 2 in 2100'] = pd.Series(np.linspace(fertiliy_by_year[last_year], 2, end_year-last_year+1), index=range(last_year, end_year+1))
fig1 = deepcopy(fig)
for key, fr in fr_projections.items():
fig1.add_trace(
go.Scatter(
x=fr.index, # Need to shift by one because they refer to "end of y" while I refer to "start of y+1" for the same quantity
y=fr,
name=f"SCENARIO - {key}",
line=dict(dash="dot")
)
)
# add projection from ISTAT
fig1.add_trace(go.Scatter(
x=df3p.index,
y=df3p["PROJMED"],
mode="lines", name="ISTAT Projection ±80 percentile", line=dict(color="black")))
fig1.add_trace(go.Scatter(
x=df3p.index,
y=df3p["PROJLOW80"],
mode="lines", showlegend=False, fill=None,line=dict(color="gray")))
fig1.add_trace(go.Scatter(
x=df3p.index,
y=df3p["PROJUPP80"],
mode="lines", showlegend=False, fill="tonexty",line=dict(color="gray")))
fig1.update_layout(
title=None,
xaxis_title="Year",
legend_title="Data Source",
margin=dict(l=0, r=0, t=20, b=0),
width=780,
height=280,
)
fig1.write_html("../images_output/fertility_scenarios.html")
last_istat_value = df6.query("CITIZENSHIP == 'TOTAL'").iloc[-1]['value']
print(f"Fertilyty from last year ({last_year}) - my estimate: {fertiliy_by_year[last_year]:.2f}, ISTAT: {last_istat_value:.2f}")
print("Fertility rate future projection: 3 scenarios")
fig1.show()
I can not compute the future newborn because I need to know the exact population in the 18-50 range to infer the number of mothers.
Conclusions¶
- the fertility rate is decreasing
- the age of motherhood is shifting ahead and it is approximable to a normal distribution over the 18-50 years range
- considering past data and future scenario the discrepancy I get between my estimate of fertility rate and the ISTAT one is negligible: much of the gap I had in provious calculation was due to neglecting the distribution of mothers' age
- I can not undertand why ISTAT projection are applying some sort of "mean reversion" instead of "trend following": is that an opinionated choice?
Follow-up¶
- Use the sigmoid model of fertility to estimate the number of new borns in the future from the number of 18-50y women
- Investigate the motivations behind ISTAT projection, and how they propagate to "too optimistic" model of the future population that make the INPS projection looking better!