import json
import math
import random
import unittest
from time import time

import pandas as pd
from leximpact_socio_fisca_simu_etat.config import Configuration

# from sklearn.metrics import r2_score
from tqdm.notebook import tqdm

from leximpact_prepare_data.calib_and_copules import (
    copules_to_df,
    get_copules_revkire,
    pandas_to_vaex,
)
from leximpact_prepare_data.monte_carlo import *
from leximpact_prepare_data.toolbase import *

config = Configuration(project_folder="leximpact-prepare-data")
tc = unittest.TestCase()
import numpy as np
import seaborn as sns
sns.set(rc={"figure.figsize": (20, 8)})

Import de la base de données ERFS-FPR

On charge un échantillon de population préparé par l'INSEE.

Pour cela on utilise le fichier qui sort de l'étape précédente.

filepath = config.get("DATA_OUT") + "03_erfs_rfr_ind" + config.get("YEAR_ERFS") + ".h5"
# print(filepath)
erfs_03 = pd.read_hdf(filepath)

Passage des individus aux foyers fiscaux

A noter : le wprm est un facteur indiquant le "poids" mathématique du foyer par rapport à la population française. Il est le même pour tous les individus d'un foyer, il ne faut donc pas le sommer quand on regroupe les individus par foyer, mais le garder tel quel.

sample_pop_ff = individus_to_foyers_fiscaux(erfs_03)

Création de copule des salaire de l'ERFS

Nous allons effectuer sur l'ERFS ce que nous faison normalement sur POTE.

sample_pop_ff.columns
vaex_df = pandas_to_vaex(sample_pop_ff)
# Le module leximpact_prepare_data.calib_and_copules s'attend à une colonne revkire et non rfr
_ = vaex_df.rename("rfr", "revkire")
%%time
# Temps de traitement : 5 secondes.
vaex_df = vaex_df.sort("revkire", ascending=True)
copules_salaire = get_copules_revkire(vaex_df, 50, "salaire_de_base", 50)
 

Utilisation de Monte-Carlo

%%time
# Temps de traitement : 12s
df, errors = apply_Monte_Carlo(
    sample_pop_ff,
    copules_salaire,
    bucket_level_col_name="rfr",
    out_col="monte_carlo_salaire_de_base",
    nb_tirage=20,
    seed=None,
    debug=False,
)

Analyse du résultat

print(
    f"Somme salaire_de_base : {df.salaire_de_base.sum():,}, somme monte_carlo_salaire_de_base : {df.monte_carlo_salaire_de_base.sum():,}. Ecart {abs((df.monte_carlo_salaire_de_base.sum()-df.salaire_de_base.sum())*100)/df.salaire_de_base.sum():.2f}%"
)
print(
    f"Moyenne salaire_de_base : {df.salaire_de_base.mean():,.2f}, moyenne monte_carlo_salaire_de_base : {df.monte_carlo_salaire_de_base.mean():,.2f}. Ecart {abs((df.monte_carlo_salaire_de_base.mean()-df.salaire_de_base.mean())*100)/df.salaire_de_base.mean():.2f}%"
)
print(
    f"Ecart-type salaire_de_base : {df.salaire_de_base.std():,.2f}, ecart-type monte_carlo_salaire_de_base : {df.monte_carlo_salaire_de_base.std():,.2f}. Ecart {abs((df.monte_carlo_salaire_de_base.std()-df.salaire_de_base.std())*100)/df.salaire_de_base.std():.2f}%"
)

Le résultat de la méthode de Monte-Carlo est très proche de la réalité, avec une erreur inférieure à 1% avec nb_tirage=20.

dfsub = df.query(
    "rfr < 250_000 and salaire_de_base < 250_000 and monte_carlo_salaire_de_base < 250_000"
)
ax = sns.scatterplot(
    data=dfsub,
    x="rfr",
    y="monte_carlo_salaire_de_base",
    alpha=0.2,
    label="monte_carlo_salaire_de_base",
)
ax = sns.scatterplot(
    data=dfsub, x="rfr", y="salaire_de_base", alpha=0.06, label="salaire_de_base", ax=ax
)
_ = ax.set_title(
    f"Evolution des salaires < 250 000 € en fonction du RFR sur l'ERFS-FPR {config.get('YEAR_ERFS')}",
    fontsize=18,
)

Q : Pourquoi la taille des bucket (zones qui forment un rectangle bleu) ne semble pas constante mais être de plus en plus grande ?

R : Parceque ils sont en nombre de personnes identiques et pas en écart de montant identique : Plus on monte dans les revenus plus il faut des fourchettes larges pour avoir le même nombre de personnes.

import matplotlib.ticker as ticker

dfsub = df.query("salaire_de_base > 0 and monte_carlo_salaire_de_base > 0")
ax = sns.histplot(
    data=dfsub[["salaire_de_base", "monte_carlo_salaire_de_base"]],
    kde=True,
    log_scale=True,
    bins=50,
)
xticks = np.logspace(0, 6.45, num=31)
xlabels = [f"{j:,.0f} €" for j in xticks]
_ = ax.set_xticks(xticks)
_ = ax.set_xticklabels(xlabels, rotation=75)
_ = ax.set_xlabel("Salaire", fontsize=16)
_ = ax.set_ylabel("Nombre de foyers", fontsize=16)
_ = ax.set_title(
    f"Histogramme des salaires sur l'ERFS-FPR {config.get('YEAR_ERFS')}\n(échelle logarithmique)",
    fontsize=18,
)

L'histogramme montre un bon respect de la distribution.

On a maintenant plus confiance dans notre méthode, voyons ce que cela donne avec les données de POTE sur les salaires.

Test avec les salaires de POTE

Import des copules de POTE

# file = "CalibPote-2019-rev_salaire.json"
file = (
    config.get("COPULES")
    + "/"
    + "ExportCopule-"
    + config.get("YEAR_ERFS")
    + "-rev_salaire.json"
)
# file ='CalibPote-1000-2019-revkire.json'
with open(file) as f:
    copule_var = json.load(f)
df_copules = copules_to_df(copule_var)
if config.get("YEAR_POTE") == "2018":
    tc.assertEqual(df_copules.sum_tranche_var.sum(), 649_078_215_971)
elif config.get("YEAR_POTE") == "2019":
    tc.assertEqual(df_copules.sum_tranche_var.sum(), 649_078_215_971)
else:
    raise Exception("Ajouter contrôles !")
pop_total = 0
for bucket in copule_var:
    pop_total += bucket["nb_foyer"]["zero"]
    pop_total += bucket["nb_foyer"]["nonzero"]
if config.get("YEAR_POTE") == "2008":
    assert pop_total == 38_487_937
sample_pop_ff["salaire_wprm"] = sample_pop_ff.wprm * sample_pop_ff.salaire_de_base
print(
    f'Somme des salaires des foyers {sample_pop_ff["salaire_wprm"].sum():,}€ Somme des wprm : {sample_pop_ff.wprm.sum():,}'
)
print("Somme des salaires dans POTE 2018 : 649,078,215,971€")

Ajout des salaires

data_to_process = [
    {
        "file": "ExportCopule-2018-rev_salaire.json",
        "column_name": "salaire_de_base_monte_carlo",
        "aggregat": 649078215971,
    }
]


def integration_data_ff(data_to_process, nb_tirage):
    for data in data_to_process:
        variable_name = data["column_name"]
        # On charge les copules
        with open(config.get("COPULES") + data["file"]) as fichier:
            dict_copules_RFR_Data = json.load(fichier)
        print(
            f"Nombre de bucket : {len(dict_copules_RFR_Data)} Nombre de sous-bucket : {len(dict_copules_RFR_Data[1]['buckets'])}"
        )
        # Chercher la distribution optimale
        # On ajoute une distribution de data dans notre échantillon de population
        # Pour chaque ligne de l'ERFS-FPR on appel calcul_distrib_from_copules en lui donnant le RFR de la ligne en cours.
        df, errors = apply_Monte_Carlo(
            df=sample_pop_ff,
            copules=dict_copules_RFR_Data,
            bucket_level_col_name="rfr",
            out_col=data["column_name"],
            nb_tirage=nb_tirage,
        )
        # Calcul de l'erreur
        value = (sample_pop_ff[variable_name] * sample_pop_ff["wprm"]).sum()
        expected = data["aggregat"]
        # print(f"On trouve {value:,} € de {variable_name} alors qu'en 2019 c'était {expected:,} €, soit {((expected-value)/expected)*100}%")
        data["error"] = ((expected - value) / expected) * 100

    return sample_pop_ff, data_to_process
sample_pop_ff, data_to_process = integration_data_ff(data_to_process, nb_tirage=10)
# sample_pop_ff
# data_to_process

Analyse des résultats

 
 
sample_pop_ff[["salaire_de_base", "salaire_de_base_monte_carlo"]].describe()
nbzero_avant = sample_pop_ff[sample_pop_ff["salaire_de_base"] == 0][
    "salaire_de_base"
].count()
nbzero_monte_carlo = sample_pop_ff[sample_pop_ff["salaire_de_base_monte_carlo"] == 0][
    "salaire_de_base_monte_carlo"
].count()
print(
    f"Nombre d'individus sans salaires avant {nbzero_avant:,}. Nombre avec Monte-Carlo {nbzero_monte_carlo:,}"
)
df_copules[["seuil_var_supp", "sum_tranche_var", "mean_tranche_var"]].describe()
def var_max(df_erfs, df_calib_pote, var_name, var_erfs, var_erfs_apres, var_pote):
    df_tmp = df_erfs.groupby(["idfoy"]).sum()
    erfs_max = df_tmp[var_erfs].max()
    erfs_max_apres = df_tmp[var_erfs_apres].max()
    print(
        f"{var_name} max dans ERFS {config.get('YEAR_ERFS')} avant traitement : {erfs_max:,.0f}€, après Monte-Carlo : {erfs_max_apres:,.0f}€,"
        + f" {var_name} max dans POTE {config.get('YEAR_ERFS')} : {df_calib_pote[var_pote].max():,.0f}€, écart {(df_calib_pote[var_pote].max()-erfs_max_apres)/df_calib_pote[var_pote].max()*100:.2f}%"
    )


def var_mean(df_erfs, df_calib_pote, var_name, var_erfs, var_erfs_apres, var_pote):
    mean_pote = df_calib_pote[var_pote].mean()
    df_erfs["tmp_var_pop_avant"] = df_erfs[var_erfs] * df_erfs["wprm"]
    erfs_mean_avant = df_erfs["tmp_var_pop_avant"].sum() / df_erfs["wprm"].sum()
    df_erfs["tmp_var_pop_apres"] = df_erfs[var_erfs_apres] * df_erfs["wprm"]
    erfs_mean_apres = df_erfs["tmp_var_pop_apres"].sum() / df_erfs["wprm"].sum()
    print(
        f"{var_name} moyen dans ERFS {config.get('YEAR_ERFS')} avant traitement : {erfs_mean_avant:,.0f}€, après Monte-Carlo : {erfs_mean_apres:,.0f}€,"
        + f" {var_name} moyen dans POTE {config.get('YEAR_ERFS')} : {mean_pote:,.0f}€, écart {(mean_pote-erfs_mean_apres)/mean_pote*100:.2f}%"
    )


def var_sum(df_erfs, df_calib_pote, var_name, var_erfs, var_erfs_apres, var_pote):
    df_erfs["tmp_var_pop_avant"] = df_erfs[var_erfs] * df_erfs["wprm"]
    sum_erfs = df_erfs["tmp_var_pop_avant"].sum()
    df_erfs["tmp_var_pop_apres"] = df_erfs[var_erfs_apres] * df_erfs["wprm"]
    erfs_sum_apres = df_erfs["tmp_var_pop_apres"].sum()
    print(
        f"{var_name} total dans ERFS {config.get('YEAR_ERFS')} avant traitement : {sum_erfs:,.0f}€, après Monte-Carlo : {erfs_sum_apres:,.0f}€,"
        + f" {var_name} total dans POTE {config.get('YEAR_ERFS')} : {df_calib_pote[var_pote].sum():,.0f}€, écart {(df_calib_pote[var_pote].sum()-erfs_sum_apres)/df_calib_pote[var_pote].sum()*100:.2f}%"
    )


var_max(
    sample_pop_ff,
    df_copules,
    "Salaire",
    "salaire_de_base",
    "salaire_de_base_monte_carlo",
    "seuil_var_supp",
)
var_mean(
    sample_pop_ff,
    df_copules,
    "Salaire",
    "salaire_de_base",
    "salaire_de_base_monte_carlo",
    "mean_tranche_var",
)
var_sum(
    sample_pop_ff,
    df_copules,
    "Salaire",
    "salaire_de_base",
    "salaire_de_base_monte_carlo",
    "sum_tranche_var",
)

Analyse du résultat

Analyse des erreurs

%%time
# Estimation de la variabilité des résultats due au random state
error_table = []
for i in range(3):
    sample_pop_ff, data_to_process = integration_data_ff(data_to_process, nb_tirage=5)
    errors = {}
    for data in data_to_process:
        errors[data["column_name"]] = data["error"]
    error_table.append(errors)

error_table = pd.DataFrame(error_table)
error_table[data["column_name"]].tolist()
error_table.salaire_de_base_monte_carlo.std()
# 0.40 @ 10 essais avec nb_tirage=10
# 0.40 @ 10 essais avec nb_tirage=20
# 0.37 @ 10 essais avec nb_tirage=120

Analyse graphique

sample_pop_ff
df = sample_pop_ff
salaire_max = 250_000
dfsub = df.query(
    "rfr < @salaire_max and salaire_de_base < @salaire_max and monte_carlo_salaire_de_base < @salaire_max and salaire_de_base_monte_carlo < @salaire_max"
)
ax = sns.scatterplot(
    data=dfsub,
    x="rfr",
    y="monte_carlo_salaire_de_base",
    alpha=0.2,
    label="monte_carlo_salaire_de_base (ERFS)",
)
ax = sns.scatterplot(
    data=dfsub,
    x="rfr",
    y="salaire_de_base_monte_carlo",
    alpha=0.2,
    label="salaire_de_base_monte_carlo (POTE)",
)
ax = sns.scatterplot(
    data=dfsub, x="rfr", y="salaire_de_base", alpha=0.06, label="salaire_de_base", ax=ax
)
_ = ax.set_title(
    f"Evolution des salaires < 250 000 € en fonction du RFR sur l'ERFS-FPR {config.get('YEAR_ERFS')}",
    fontsize=18,
)

Q : Pourquoi la taille des buckets venant de POTE (orange) ne resemble pas à celle de ceux venant de l'ERFS (zones qui forment un rectangle bleu) ?

R :

import matplotlib.ticker as ticker

dfsub = df.query(
    "salaire_de_base > 0 and monte_carlo_salaire_de_base > 0 and salaire_de_base_monte_carlo > 0"
)
ax = sns.histplot(
    data=dfsub[
        [
            "salaire_de_base",
            "monte_carlo_salaire_de_base",
            "salaire_de_base_monte_carlo",
        ]
    ],
    kde=True,
    log_scale=True,
    bins=50,
)
xticks = np.logspace(0, 6.45, num=31)
xlabels = [f"{j:,.0f} €" for j in xticks]
_ = ax.set_xticks(xticks)
_ = ax.set_xticklabels(xlabels, rotation=75)
_ = ax.set_xlabel("Salaire", fontsize=16)
_ = ax.set_ylabel("Nombre de foyers", fontsize=16)
_ = ax.set_title(
    f"Histogramme des salaires sur l'ERFS-FPR {config.get('YEAR_ERFS')}\n(échelle logarithmique)",
    fontsize=18,
)
sample_pop_ff_sorted = sample_pop_ff.sort_values(by=["rfr"]).reset_index()
# sns.scatterplot(data=sample_pop_ff_sorted, x=sample_pop_ff_sorted.index, y="salaire_de_base", alpha=0.9)
sns.scatterplot(
    data=sample_pop_ff_sorted[["salaire_de_base", "salaire_de_base_monte_carlo"]],
    alpha=0.2,
)
# sns.scatterplot(data=sample_pop_ff, x=sample_pop_ff.index, y="rfr")
bins = 100
sample_pop_ff_sorted_filtered = sample_pop_ff_sorted.query(
    "10_000 < salaire_de_base < 100_000 and 10_000 <  salaire_de_base_monte_carlo < 100_000"
)
_ = sample_pop_ff_sorted_filtered.salaire_de_base.hist(bins=bins)
_ = sample_pop_ff_sorted_filtered.salaire_de_base_monte_carlo.hist(bins=bins, alpha=0.8)