# IPCH alternatif pour la France, janvier 1996 – T2 2026 (script du 19 octobre 2026).
#
# Principe : on reconstruit l'IPCH français à partir de ses postes à trois chiffres (ECOICOP v2, Eurostat
# prc_hicp_minr et prc_hicp_iw), par la même formule que l'IPCH (Laspeyres chaîné annuellement sur décembre),
# puis on modifie trois blocs :
#   1. loyers effectifs (CP041) : dérive ajoutée au secteur libre (scénarios bas, central, haut) ;
#   2. effets qualité : sur les postes à fort ajustement qualité (S1, S2 de effet-qualite/simulation),
#      l'évolution française est remplacée par l'évolution allemande relative au panier témoin ;
#   3. logement des propriétaires occupants, absent de l'IPCH, ajouté selon trois approches :
#      A acquisitions nettes (indice OOH d'Eurostat), B paiements (B1 intérêts, B2 intérêts et capital),
#      C loyers imputés.
# Les écarts avec l'IPCH publié ne viennent ainsi que des choix de méthode.
#
# Lancer : DATA_DIR=~/Dropbox/website/data python3 ipch_alternatif.py
# Sorties : resultats/*.csv
import os, sys, re, pathlib, numpy as np, pandas as pd

ICI = pathlib.Path(__file__).resolve().parent
DATA = pathlib.Path(os.environ.get('DATA_DIR', ICI / '../../website/data')).expanduser()
os.environ.setdefault('EUROSTAT_DIR', str(DATA / 'eurostat'))
sys.path.insert(0, str(ICI / '../effet-qualite/simulation'))
from commun import charger, S1, S2, S3  # noqa: E402

LOC = ICI / 'data'
RES = ICI / 'resultats'; RES.mkdir(exist_ok=True)
DERNIER = '2026-06'                     # fin T2 2026
AMORT = 0.055                           # amortissement annuel du capital, en part de l'encours (5 à 6,5 %)

MOIS = [f'{y}-{m:02d}' for y in range(1996, 2027) for m in range(1, 13) if f'{y}-{m:02d}' <= DERNIER]
def mois(y): return [t for t in MOIS if t.startswith(str(y))]
def base_mois(y): return '1996-01' if y == 1996 else f'{y-1}-12'
def q_de_m(t): return f'{t[:4]}Q{(int(t[5:]) - 1) // 3 + 1}'
def trim_vers_mois(s):
    """Série trimestrielle 'AAAAQn' -> mensuelle (valeur du trimestre pour ses trois mois)."""
    return pd.Series({t: s.get(q_de_m(t), np.nan) for t in MOIS})

# ---------------------------------------------------------------- 1. IPCH par postes
M, W, lab = charger(pays=['FR', 'DE'])
POSTES = [c for c in W.columns if re.fullmatch(r'CP\d{3}', c) and W.loc['FR', c].max() > 0]

def relatifs(g, y, cols):
    return M.loc[g].reindex(mois(y))[cols] / M.loc[(g, base_mois(y)), cols]

# ---------------------------------------------------------------- 2. données logement
ins = lambda f: pd.read_parquet(DATA / 'insee' / f)

# 2a. part du secteur libre dans l'indice des loyers (régression des variations trimestrielles)
L = ins('INDICES_LOYERS.parquet')
L = L[(L.NATURE == 'INDICE') & (L.REF_AREA == 'FM') & (L.CORRECTION == 'BRUT')]
L = L.pivot_table(index='TIME_PERIOD', columns='SECTEURS_LOYERS', values='OBS_VALUE').sort_index()
d = np.log(L).diff().dropna()
y_, x_ = d.TOUS_SECTEURS - d.SECTEUR_SOCIAL, d.SECTEUR_LIBRE - d.SECTEUR_SOCIAL
S_PRIV = float((x_ * y_).sum() / (x_ ** 2).sum())

SCEN_LOYERS = {  # biais sur le secteur libre, en point par an
    'bas':     lambda y: 0.4,
    'central': lambda y: 0.4 if y <= 2019 else 1.0,
    'haut':    lambda y: 1.0 if y <= 2019 else 2.0,
}

# 2b. comptes nationaux : dépense de consommation hors SIFIM et hors loyers imputés (champ proche de l'IPCH)
C = ins('T_CONSO_EFF_FONCTION.parquet')
C = C[C.sheet == 'MEURcour'].pivot_table(index='year', columns='code', values='value')
DEP_IPCH = C.CPDEPHSI - C.CP042
IMPUTES = C.CP042

# 2c. acquisitions nettes : FBCF en construction des ménages hors EI x taux de propriétaires / part des
#     postes « neufs » dans l'indice OOH (méthode de la BCE, Statistics Paper n° 47, éq. 6-8)
F = ins('CNA-2020-FBCF-SI.parquet')
F = F[(F['SECT-INST'] == 'S14B') & (F.UNIT_MEASURE == 'EUROS_COURANTS')]
F = F.pivot_table(index='TIME_PERIOD', columns='CNA_PRODUIT', values='OBS_VALUE'); F.index = F.index.astype(int)
FZ = F['A17-FZ'].copy()
for a in range(int(FZ.last_valid_index()) + 1, 2026):               # prolongé par la FBCF totale des ménages
    FZ[a] = FZ[a - 1] * F.NNTOTAL[a] / F.NNTOTAL[a - 1]
OR = pd.read_parquet(DATA / 'eurostat/ilc_lvho02.parquet')
OR = OR[(OR.geo == 'FR') & (OR.tenure == 'OWN') & (OR.hhcomp == 'TOTAL') & (OR.rskpovth == 'TOTAL')]
OR = OR.assign(time=OR.time.astype(int)).set_index('time')['values'] / 100
OOW = pd.read_parquet(DATA / 'eurostat/prc_hpi_ooinw.parquet'); OOW = OOW[OOW.geo == 'FR']
OOW = OOW.assign(time=OOW.time.astype(int)).pivot_table(index='time', columns='expend', values='values')
def val(s, a):  # valeur de l'année a, sinon la plus proche
    s = s.dropna(); return s.get(a, s.iloc[0] if a < s.index.min() else s.iloc[-1])
E_OOH = pd.Series({a: FZ[a] * val(OR, a) / (val(OOW.DW_ACQ_NEW, a) / 1000) for a in range(1995, 2026)})

# 2d. prix OOH : Eurostat (2009 T1 – 2026 T2), rétropolé par un indice composite avant 2009
O = pd.read_parquet(DATA / 'eurostat/prc_hpi_ooq.parquet')
O = O[(O.geo == 'FR') & (O.unit == 'I15_Q')].pivot_table(index='time', columns='expend', values='values')
P = ins('IPLA-IPLNA-2015.parquet'); P = P[(P.CORRECTION == 'BRUT') & (P.REF_AREA == 'FM')]
P = P.pivot_table(index='TIME_PERIOD', columns='INDICATEUR', values='OBS_VALUE'); P.index = P.index.str.replace('-', '')
icc = pd.read_csv(LOC / 'icc.csv', index_col=0).icc; icc.index = icc.index.str.replace('-', '')
hicp_q = lambda p: M.loc['FR', p].groupby(lambda t: q_de_m(t)).mean()
neufs = P.IPLN.copy()                                                 # prix des logements neufs, 2000 T1+
neufs = neufs.combine_first(P.IPLA_E * neufs['2000Q1'] / P.IPLA_E['2000Q1'])  # logements anciens avant 2000
w09 = OOW.loc[2009]
comp = pd.DataFrame({'NEWP': neufs, 'NEWSB': icc, 'OTH': hicp_q('TOTAL'), 'OWN': hicp_q('CP043')}).loc['1995Q1':'2009Q1']
comp = comp / comp.loc['2009Q1']
proxy = (comp.NEWP * w09.DW_ACQ_NEWP + comp.NEWSB * w09.DW_ACQ_NEWSB + comp.OTH * w09.DW_ACQ_OTH + comp.OWN * w09.DW_OWN) / 1000
P_OOH = O.TOTAL.combine_first(proxy * O.TOTAL['2009Q1'])

# 2e. paiements : encours et taux des crédits à l'habitat (Banque de France), prix des logements
B = pd.read_parquet(DATA / 'bdf/immobilier.parquet'); B['t'] = pd.to_datetime(B.date).dt.strftime('%Y-%m')
B = B.pivot_table(index='t', columns='variable', values='value')
ENCOURS = B['BSI1.M.FR.N.R.A220Z.A.1.U6.2254FR.Z01.E']                # millions d'euros
MIR = pd.read_parquet(DATA / 'bdf/MIR.parquet'); MIR = MIR[MIR.variable == 'MIR.M.FR.B.A22.J.R.A.2250.EUR.O']
TAUX = MIR.assign(t=pd.to_datetime(MIR.date).dt.strftime('%Y-%m')).set_index('t').value.sort_index()
# taux moyen sur l'encours des crédits à l'habitat de plus de cinq ans, depuis janv. 2003
oat = pd.read_parquet(DATA / 'eurostat/irt_lt_mcby_m.parquet'); oat = oat[oat.geo == 'FR']
oat = oat.assign(t=oat.time.str.replace('M', '-')).set_index('t')['values'].sort_index()
lisse = oat.ewm(halflife=60).mean()                                   # l'encours se renouvelle lentement
avant = lisse[lisse.index < TAUX.index[0]]
TAUX = pd.concat([TAUX.iloc[0] + avant - lisse[TAUX.index[0]], TAUX])  # rétropolation 1996-2002 par l'OAT lissée
HPI = P.IPLNA.combine_first(P.IPLA_E * P.IPLNA['2000Q1'] / P.IPLA_E['2000Q1'])
HPI_M = trim_vers_mois(HPI)
TAUX_M = TAUX.reindex(MOIS).ffill()
P_INT = TAUX_M * HPI_M                                                # convention de l'indice RPI britannique
moy = lambda s, a: s[[t for t in s.index if t.startswith(str(a))]].mean()
INTERETS = pd.Series({a: moy(ENCOURS, a) * moy(TAUX, a) / 100 for a in range(1995, 2026)})
CAPITAL = pd.Series({a: AMORT * moy(ENCOURS, a) for a in range(1995, 2026)})

# poids (pour mille de l'IPCH), année y : dépenses de y-1 rapportées à la dépense de y-1
POIDS = pd.DataFrame({
    'A_acquisitions_nettes': 1000 * E_OOH / DEP_IPCH,
    'B1_interets': 1000 * INTERETS / DEP_IPCH,
    'B2_capital': 1000 * CAPITAL / DEP_IPCH,
    'C_loyers_imputes': 1000 * IMPUTES / DEP_IPCH,
}).loc[1995:2025]
POIDS.index = POIDS.index + 1
POIDS['IPCH_loyers_effectifs'] = W.loc['FR', 'CP041'].reindex(POIDS.index)
POIDS.round(1).to_csv(RES / 'poids_logement.csv')

P_OOH_M = trim_vers_mois(P_OOH)
def rel_serie(s, y):
    return s.reindex(mois(y)) / s[base_mois(y)]

# ---------------------------------------------------------------- 3. chaînage
def chaine(loyers=None, qualite=None, ooh=None, forfait=0.0):
    """loyers : scénario de SCEN_LOYERS ; qualite : liste de postes remplacés par l'Allemagne (relatif) ;
    ooh : None, 'A', 'A2', 'B1', 'B2' ou 'C' ; forfait : dérive qualité sur l'indice d'ensemble, en point par an
    (borne de effet-qualite.qmd : 0,3 point pour les remplacements en cours d'année, environ 0,6 avec le renouvellement)."""
    niveau, out = 100.0, {}
    for y in range(1996, int(DERNIER[:4]) + 1):
        mm = mois(y); R = relatifs('FR', y, POSTES).copy(); w = W.loc[('FR', y), POSTES].fillna(0)
        k = pd.Series([int(t[5:]) - (1 if y == 1996 else 0) for t in mm], index=mm)
        r041 = R['CP041'].copy()
        derive = np.exp(SCEN_LOYERS[loyers](y) / 100 * k / 12) if loyers else pd.Series(1.0, index=mm)
        if qualite:
            RD = relatifs('DE', y, POSTES); wD = W.loc[('DE', y), POSTES].fillna(0)
            tem = [p for p in POSTES if p not in S3 and R[p].notna().all() and RD[p].notna().all()]
            adj = ((R[tem] * w[tem]).sum(axis=1) / w[tem].sum()) / ((RD[tem] * wD[tem]).sum(axis=1) / wD[tem].sum())
            for p in qualite:
                R[p] = (RD[p] * adj).where(RD[p].notna(), R[p])
        R['CP041'] = r041 * (1 + S_PRIV * (derive - 1))
        ok = R.notna().all()
        num = (R.loc[:, ok] * w[ok]).sum(axis=1); den = w[ok].sum()
        if forfait: num = num * np.exp(forfait / 100 * k / 12)   # sur le panier de l'IPCH seulement
        extra = []
        if ooh == 'A': extra = [(POIDS.A_acquisitions_nettes[y], rel_serie(P_OOH_M, y))]
        if ooh == 'A2': extra = [(POIDS.A_acquisitions_nettes[y], rel_serie(HPI_M, y))]  # prix des logements, terrain compris
        if ooh in ('B1', 'B2'): extra = [(POIDS.B1_interets[y], rel_serie(P_INT, y))]
        if ooh == 'B2': extra += [(POIDS.B2_capital[y], rel_serie(HPI_M, y))]
        if ooh == 'C': extra = [(POIDS.C_loyers_imputes[y], r041 * derive)]  # biais entier : loyers de marché
        for wy, Ry in extra:
            num = num + wy * Ry; den += wy
        agg = num / den
        for t, v in agg.items(): out[t] = niveau * v
        niveau = out[mm[-1]]
    return pd.Series(out)

PUBLIE = M.loc['FR', 'TOTAL'].reindex(MOIS); PUBLIE = PUBLIE / PUBLIE['1996-01'] * 100

VARIANTES = {
    'IPCH reconstruit': {},
    'Loyers, scénario bas': dict(loyers='bas'),
    'Loyers, scénario central': dict(loyers='central'),
    'Loyers, scénario haut': dict(loyers='haut'),
    'Qualité, S1 (Allemagne relatif)': dict(qualite=S1),
    'Qualité, S2 avec véhicules': dict(qualite=S2),
    'Logement A, acquisitions nettes': dict(ooh='A'),
    'Logement A2, acquisitions au prix des logements': dict(ooh='A2'),
    'Logement B1, paiements (intérêts)': dict(ooh='B1'),
    'Logement B2, paiements (intérêts et capital)': dict(ooh='B2'),
    'Logement C, loyers imputés': dict(ooh='C'),
}
for x, nom in [(None, 'sans logement des propriétaires'), ('A', 'A acquisitions nettes'),
               ('A2', 'A2 acquisitions au prix des logements'), ('B1', 'B1 intérêts'),
               ('B2', 'B2 intérêts et capital'), ('C', 'C loyers imputés')]:
    VARIANTES[f'Ensemble bas, {nom}'] = dict(loyers='bas', qualite=S1, ooh=x)
    VARIANTES[f'Ensemble central, {nom}'] = dict(loyers='central', qualite=S2, ooh=x)
    VARIANTES[f'Ensemble haut, {nom}'] = dict(loyers='haut', qualite=S2, ooh=x)

for q in (0.4, 0.6):
    VARIANTES[f'Qualité {q:.1f} pt seule'.replace('.', ',')] = dict(forfait=q)

# Grille de l'IPCH alternatif : loyers du scénario central ; trois hypothèses de qualité ; six traitements
# du logement des propriétaires. Codes des colonnes du fichier ipch-alternatif.csv entre parenthèses.
QUALITES = {'Allemagne': ('DE', dict(qualite=S2)), '0,4': ('04', dict(forfait=0.4)), '0,6': ('06', dict(forfait=0.6))}
LOGEMENTS = {'aucun': None, 'A': 'A', 'A2': 'A2', 'B1': 'B1', 'B2': 'B2', 'C': 'C'}
GRILLE = {}
for qn, (qc, qd) in QUALITES.items():
    for ln, lx in LOGEMENTS.items():
        nom = f'Alternatif, qualité {qn}, logement {ln}'
        VARIANTES[nom] = dict(loyers='central', ooh=lx, **qd)
        GRILLE[nom] = f'alt_q{qc}_{ln}'

SERIES = pd.DataFrame({'IPCH publié': PUBLIE, **{k: chaine(**v) for k, v in VARIANTES.items()}})
SERIES.round(3).to_csv(RES / 'series_mensuelles.csv', index_label='mois')

# ---------------------------------------------------------------- 3 bis. fichier de l'IPCH alternatif
# Indices mensuels, base 100 en moyenne 2025 (comme l'IPCH d'Eurostat), et glissements annuels.
B25 = SERIES[['IPCH publié', *GRILLE]].copy()
B25 = B25 / B25.loc[[t for t in B25.index if t.startswith('2025')]].mean() * 100
B25 = B25.rename(columns={'IPCH publié': 'ipch_publie', **GRILLE})
B25.round(3).to_csv(ICI / 'ipch-alternatif.csv', index_label='mois')
GA = (B25 / B25.shift(12) - 1) * 100
GA.dropna().round(2).to_csv(ICI / 'ipch-alternatif-glissement-annuel.csv', index_label='mois')
MA = B25.groupby(lambda t: t[:4]).mean(); MA = MA[MA.index <= '2025']
pd.concat({'indice': MA.round(3), 'inflation_pct': (MA.pct_change() * 100).round(2)}, axis=1) \
  .to_csv(ICI / 'ipch-alternatif-annuel.csv', index_label='annee')
pd.DataFrame({'code': list(GRILLE.values()), 'variante': list(GRILLE)}).to_csv(ICI / 'ipch-alternatif-codes.csv', index=False)

# ---------------------------------------------------------------- 4. tableaux
FIN = ('T2 2026', ['2026-04', '2026-05', '2026-06'], 2026.375)
DEBUTS = [('1996', mois(1996), 1996.5), ('2000', mois(2000), 2000.5), ('2008', mois(2008), 2008.5),
          ('T2 2017', ['2017-04', '2017-05', '2017-06'], 2017.375),
          ('T4 2019', ['2019-10', '2019-11', '2019-12'], 2019.875),
          ('T4 2021', ['2021-10', '2021-11', '2021-12'], 2021.875)]
lignes = []
for nom, mm, t0 in DEBUTS:
    n = FIN[2] - t0
    for v in SERIES.columns:
        c = SERIES[v][FIN[1]].mean() / SERIES[v][mm].mean()
        c0 = SERIES['IPCH publié'][FIN[1]].mean() / SERIES['IPCH publié'][mm].mean()
        lignes.append(dict(debut=nom, variante=v, annees=round(n, 2), cumul_pct=round((c - 1) * 100, 1),
                           par_an_pct=round((c ** (1 / n) - 1) * 100, 2),
                           ecart_pt_par_an=round((c ** (1 / n) - c0 ** (1 / n)) * 100, 2)))
T = pd.DataFrame(lignes); T.to_csv(RES / 'tableau_periodes.csv', index=False)

ann = SERIES.groupby(lambda t: t[:4]).mean()
ann = ann[ann.index <= '2025']
(ann.pct_change() * 100).round(2).dropna().to_csv(RES / 'inflation_annuelle.csv', index_label='annee')

repro = ((SERIES['IPCH reconstruit'] / SERIES['IPCH publié'] - 1) * 100)
pd.Series({'part_secteur_libre': round(S_PRIV, 3), 'ecart_max_reproduction_pct': round(repro.abs().max(), 3),
           'ecart_fin_reproduction_pct': round(repro.iloc[-1], 3)}).to_csv(RES / 'parametres.csv', header=['valeur'])

if __name__ == '__main__':
    pd.set_option('display.width', 250); pd.set_option('display.max_rows', 500)
    print('Part du secteur libre :', round(S_PRIV, 3), '| reproduction, écart max', round(repro.abs().max(), 3), '%')
    print(POIDS.round(1).iloc[::5].to_string())
    print(T.pivot_table(index='variante', columns='debut', values='par_an_pct', sort=False)
          [[d[0] for d in DEBUTS]].to_string())