#!/usr/bin/env python3
"""Reproduce the exploratory Taipei study from three public, frozen CSVs.

python run_econometrics_v9.py --input-dir . --output-dir reproduced
Dependencies: numpy, pandas, scipy, statsmodels, openpyxl; matplotlib for figures.
No network, no missing-value imputation, no data-driven choice of main model.
The historical econometrics.py and its green-tea data are NOT used.
"""
from __future__ import annotations
import argparse, hashlib, itertools, json, math, platform, sys, warnings
from datetime import datetime, timezone
from pathlib import Path
import numpy as np
import pandas as pd
import scipy
from scipy.spatial.distance import cdist
import statsmodels
import statsmodels.api as sm
from statsmodels.stats.multitest import multipletests
from statsmodels.stats.outliers_influence import variance_inflation_factor

SEED = 20260913
INPUTS = ('complete_cases.csv', 'competition_neighbors.csv', 'hotel_neighbors.csv')
COMPETITION = {
    'c_tri_3': ('Triangle · 3 min', 'Σ max(0, 1−t/3)'),
    'c_tri_5': ('Triangle · 5 min', 'Σ max(0, 1−t/5)'),
    'c_tri_10': ('Triangle · 10 min', 'Σ max(0, 1−t/10)'),
    'c_count_3': ('Nombre · 3 min', 'Σ 1(t<3)'),
    'c_count_5': ('Nombre · 5 min', 'Σ 1(t<5)'),
    'c_count_10': ('Nombre · 10 min', 'Σ 1(t<10)'),
    'c_exp_10': ('Exponentiel · 10 min', 'Σ exp(−t/5) × 1(t<10)'),
    'c_gravity_10': ('Gravitaire · 10 min', 'Σ [1/(1+t)] × 1(t<10)'),
    'c_nearest_10': ('Plus proche · 10 min', '10−min(10, min_j t_ij) ; 0 si aucun voisin <10'),
}
HOTELS = {
    'none': ('Sans hôtel', 'Aucun contrôle hôtelier'),
    'h_tri_5': ('Triangle · 5 min', 'ln(1 + Σ max(0, 1−t/5))'),
    'h_tri_10': ('Triangle · 10 min', 'ln(1 + Σ max(0, 1−t/10))'),
    'h_count_5': ('Nombre · 5 min', 'ln(1 + Σ 1(t<5))'),
    'h_count_10': ('Nombre · 10 min', 'ln(1 + Σ 1(t<10))'),
    'h_exp_10': ('Exponentiel · 10 min', 'ln(1 + Σ exp(−t/5) × 1(t<10))'),
    'h_gravity_10': ('Gravitaire · 10 min', 'ln(1 + Σ [1/(1+t)] × 1(t<10))'),
    'h_nearest_10': ('Plus proche · 10 min', '10−min(10, min_h t_ih) ; 0 si aucun hôtel <10'),
}
NONSPATIAL = {
    'rating': ('Note moyenne', 'rating', 'Note /5 de la plateforme documentée.'),
    'log_reviews': ('Nombre d’avis (log)', 'review_count', 'ln(1+nombre exact d’avis) ; réputation et demande réalisées, potentiellement endogènes.'),
    'hours_per_week': ('Heures par semaine', 'hours_per_week', 'Amplitude hebdomadaire des horaires documentés ; pas une mesure de qualité de service.'),
    'log_menu': ('Largeur du menu (log)', 'menu_n_products', 'ln(nombre de produits nommés) ; variantes de taille parfois comptées séparément selon la plateforme.'),
    'menu_has_bubble': ('Perles au menu', 'menu_has_bubble', 'Mention observée de perles/garnitures de ce type.'),
    'menu_has_milk_tea': ('Thé au lait au menu', 'menu_has_milk_tea', 'Mention observée de thé au lait.'),
    'menu_has_fresh_milk': ('Lait frais au menu', 'menu_has_fresh_milk', 'Mention observée de lait frais.'),
    'menu_has_fruit_tea': ('Thé fruité au menu', 'menu_has_fruit_tea', 'Mention observée de boissons au fruit.'),
    'menu_has_plant_milk': ('Lait végétal au menu', 'menu_has_plant_milk', 'Mention observée d’option au lait végétal.'),
    'menu_claims_organic': ('Mention bio', 'menu_claims_organic', 'Allégation dans le menu ; aucune certification indépendante déduite.'),
    'menu_claims_high_mountain': ('Mention haute montagne', 'menu_claims_high_mountain', 'Allégation présente dans le menu ; ne qualifie pas nécessairement le thé de référence.'),
    'menu_brewing_claim': ('Mention de préparation', 'brewing_claims', 'Au moins un terme de préparation déjà relevé dans le menu (par exemple extraction ou infusion à froid).'),
    'menu_origin_claim': ('Mention d’origine', 'tea_origin_claims', 'Au moins une origine géographique déjà relevée dans le menu ; pas nécessairement celle du thé de référence.'),
    'earl_grey': ('Variante Earl Grey', 'reference_tea_family', 'Produit de référence classé Earl Grey, deux cas ; dimension de recette et non preuve de qualité supérieure.'),
}
CONTEXT = {
    'youth_share': ('Part des 15–34 ans', 'village_youth_share_15_34', 'Part, de 0 à 1, du village administratif ; juin 2026.'),
    'log_mrt_distance': ('Distance au métro (log)', 'distance_nearest_mrt_m', 'ln(1+distance géodésique en mètres à l’entrée la plus proche), pas un temps de marche.'),
    'log_mrt_footfall': ('Flux du métro (log)', 'nearest_mrt_daily_footfall', 'ln(1+entrées+sorties quotidiennes moyennes de la station dont le centroïde des entrées est le plus proche), juin 2026 ; pas des visiteurs uniques.'),
    'log_income': ('Revenu du district 2023 (log)', 'district_disposable_income_twd', 'ln(revenu disponible du ménage dans le district), observations de 2023 (année ROC 112) ; proxy agrégé et antérieur à 2026.'),
    'land_index': ('Indice foncier 2012–2013 · historique', 'district_commercial_land_index', 'Indice d’évolution du 01/10/2012 au 31/03/2013 (ROC 1011001~1020331). Très ancien : diagnostic archivé, hors modèle retenu pour 2026 ; ni niveau de loyer ni coût actuel.'),
    'log_osm_food': ('Alimentation OSM (log)', 'osm_food_service_400m', 'ln(1+POI alimentaires OSM à 400 m) ; couverture cartographique inégale, pas une fréquentation observée.'),
    'log_osm_retail': ('Commerces OSM (log)', 'osm_retail_400m', 'ln(1+commerces OSM à 400 m) ; concentration commerciale pouvant être déterminée conjointement.'),
    'log_osm_education': ('Éducation OSM (log)', 'osm_education_400m', 'ln(1+lieux d’éducation OSM à 400 m) ; ni inscriptions ni fréquentation.'),
    'log_osm_tourism': ('Tourisme OSM (log)', 'osm_tourism_400m', 'ln(1+POI touristiques OSM à 400 m), pouvant inclure des hôtels ; proxy de couverture.'),
    'log_osm_markets': ('Marchés OSM (log)', 'osm_markets_400m', 'ln(1+marchés OSM à 400 m).'),
}
MAIN_CONTROLS = ['log_density', 'h_tri_10', 'rating', 'log_menu']
DIRECT_CHANNELS = ('direct_order', 'merchant_menu', 'in_store')


def clean(x):
    if isinstance(x, dict): return {str(k): clean(v) for k, v in x.items()}
    if isinstance(x, (list, tuple, np.ndarray)): return [clean(v) for v in x]
    if isinstance(x, np.generic): return clean(x.item())
    if isinstance(x, float) and not math.isfinite(x): return None
    return x


def save_json(path, obj):
    path.write_text(json.dumps(clean(obj), ensure_ascii=False, separators=(',', ':'), allow_nan=False)+'\n')


def aggregate_times(values, prefix):
    t = np.asarray(values, dtype=float) / 60
    if not np.isfinite(t).all(): raise ValueError('A routed contributor is missing; no zero imputation allowed.')
    out = {}
    for limit in ((3, 5, 10) if prefix == 'c' else (5, 10)):
        out[f'{prefix}_tri_{limit}'] = float(np.maximum(0, 1-t/limit).sum())
        out[f'{prefix}_count_{limit}'] = int((t < limit).sum())
    inside = t[t < 10]
    out[f'{prefix}_exp_10'] = float(np.exp(-inside/5).sum())
    out[f'{prefix}_gravity_10'] = float((1/(1+inside)).sum())
    out[f'{prefix}_nearest_10'] = float(10-inside.min()) if inside.size else 0.
    return out


def prepare(input_dir):
    raw = pd.read_csv(input_dir/INPUTS[0])
    assert len(raw) > 0 and raw.shop_id.nunique() == len(raw)
    assert raw.sample_eligible.all() and raw.within_taipei.all()
    assert set(raw.reference_channel) <= set(DIRECT_CHANNELS) | {'delivery'}
    assert (raw.reference_price_twd > 0).all()
    df = raw.copy()
    dictionary = []
    for file, prefix in ((INPUTS[1], 'c'), (INPUTS[2], 'h')):
        neighbors = pd.read_csv(input_dir/file)
        other = 'competitor_id' if prefix == 'c' else 'hotel_id'
        assert not neighbors.duplicated(['shop_id', other]).any()
        if prefix == 'c': assert not (neighbors.shop_id == neighbors.competitor_id).any()
        for ix, row in raw.iterrows():
            ns = neighbors[neighbors.shop_id == row.shop_id]
            values = aggregate_times(ns.walk_s.to_numpy(), prefix)
            for key, value in values.items():
                df.loc[ix, key+'_raw' if prefix == 'h' and 'nearest' not in key else key] = value
            for t in ((3, 5, 10) if prefix == 'c' else (5, 10)):
                expected = row[f'ci_{t}min' if prefix == 'c' else f'hotel_access_{t}min']
                assert abs(values[f'{prefix}_tri_{t}']-expected) < .00011, (row.shop_id, t)
                assert values[f'{prefix}_count_{t}'] == row[f'count_{t}min' if prefix == 'c' else f'hotel_count_{t}min']
    for key in HOTELS:
        if key not in ('none', 'h_nearest_10'): df[key] = np.log1p(df[key+'_raw'])
    df['log_price'] = np.log(df.reference_price_twd)
    df['delivery_price'] = (df.reference_channel == 'delivery').astype(float)
    df['rating_nidin'] = (df.rating_platform == 'nidin').astype(float)
    df['log_density'] = np.log(df.village_population_density_km2)
    df['log_reviews'] = np.log1p(df.review_count)
    df['log_menu'] = np.log(df.menu_n_products)
    df['menu_brewing_claim'] = df.brewing_claims.map(lambda v: float(bool(json.loads(v))))
    df['menu_origin_claim'] = df.tea_origin_claims.map(lambda v: float(bool(json.loads(v))))
    df['earl_grey'] = (df.reference_tea_family == 'earl_grey').astype(float)
    df['size_L'] = (df.reference_size_label == 'L').astype(float)
    df['h_capacity_10'] = np.log1p(df.hotel_capacity_access_10min)
    df['h_rooms_10'] = np.log1p(df.hotel_rooms_access_10min)
    df['log_ci_5'] = np.log1p(df.c_tri_5)
    df['ci_squared'] = df.c_tri_5**2
    df['ci_delivery'] = df.c_tri_5 * df.delivery_price
    for key, (_, source, _) in CONTEXT.items():
        df[key] = np.log1p(df[source]) if key.startswith('log_') else df[source]
    # Positive income is in ln, not ln(1+x), as specified in the plan.
    df['log_income'] = np.log(df.district_disposable_income_twd)
    for key in NONSPATIAL: df[key] = df[key].astype(float)
    def describe(key, label, category, source, definition, role='candidate'):
        s = df[key]
        is_binary = set(s.dropna().unique()).issubset({0, 1, False, True})
        dictionary.append(dict(variable=key, label=label, category=category, source_fields=source,
            definition=definition, n_observed=int(s.notna().sum()), n_unique=int(s.nunique()),
            positive_cases=int(s.sum()) if is_binary else None,
            rare_binary=bool(is_binary and s.nunique() == 2 and min(s.sum(), len(s)-s.sum()) < 5),
            role=role, missing_policy='No imputation; subset explicitly named for extensions'))
    describe('log_price', 'Prix du gobelet (log)', 'outcome', 'reference_price_twd', 'ln(prix régulier affiché en TWD) ; formats non harmonisés.', 'outcome')
    for key, (label, formula) in COMPETITION.items(): describe(key, label, 'competition', INPUTS[1]+': walk_s/60', formula+' ; tous voisins recensés, chaînes incluses.')
    for key, (label, formula) in HOTELS.items():
        if key != 'none': describe(key, label, 'hotels', INPUTS[2]+': walk_s/60', formula+' ; hébergements autorisés, pas des touristes observés.')
    for key, (label, source, definition) in NONSPATIAL.items(): describe(key, label, 'non_spatial', source, definition)
    describe('menu_claims_single_origin', 'Mention single origin', 'non_spatial_constant', 'menu_claims_single_origin', 'Zéro dans le corpus selon le codage lexical existant ; aucune variance exploitable.', 'not_estimated_constant')
    describe('delivery_price', 'Canal de prix livraison', 'measurement', 'reference_channel', 'Prix Foodpanda de livraison=1, prix publié par le marchand=0 ; ne décrit pas l’absence de livraison de la boutique.', 'secondary_measurement_sensitivity_only')
    describe('rating_nidin', 'Note issue de Nidin', 'measurement', 'rating_platform', 'Source Nidin ou Google documentée ; contrôle de source, pas une caractéristique de service.', 'sensitivity')
    describe('reference_size_label', 'Format commercial', 'measurement', 'reference_size_label', 'M, L, 700 ml, 12 oz ou standard ; standard signifie format non publié.', 'sensitivity_M_L_only')
    describe('log_density', 'Densité résidentielle (log)', 'context', 'village_population_density_km2', 'ln(population par km² du village administratif), juin 2026.', 'mandatory_grid')
    for key, (label, source, definition) in CONTEXT.items(): describe(key, label, 'context', source, definition)
    for key, source, label in [('h_capacity_10', 'hotel_capacity_access_10min', 'Capacité hôtelière pondérée'), ('h_rooms_10', 'hotel_rooms_access_10min', 'Chambres hôtelières pondérées')]:
        describe(key, label, 'hotel_extension', source, 'ln(1+somme capacité ou chambres × poids triangulaire à 10 min) ; inconnu si un hôtel contributeur manque.', 'matched_subset')
    assert len(NONSPATIAL) == 14
    # Optional context extensions can be unavailable at a newly collected origin.
    # fit_model reports those entire specifications as missing_values_no_imputation.
    assert all(df[k].notna().all() for k in [*COMPETITION, *NONSPATIAL, 'log_density', *list(HOTELS)[1:]])
    return df, pd.DataFrame(dictionary)


def effect(beta, amount):
    return float(100*np.expm1(beta*amount))


def fit_model(df, spatial, controls, sds, model_id, family, label='', outcome='log_price', extra_spatial=()):
    spatial_terms = [spatial, *extra_spatial]
    features = spatial_terms + list(controls)
    row = dict(model_id=model_id, family=family, label=label, spatial=spatial, hotel=next((v for v in controls if v in HOTELS or v in ['h_capacity_10','h_rooms_10']), 'none'),
        controls='|'.join(controls), spatial_terms='|'.join(spatial_terms), outcome=outcome, n=len(df),
        parameters=len(features)+1, status='not_estimable', reason='')
    if len(set(features)) != len(features): raise ValueError('Duplicated feature')
    X = sm.add_constant(df[features].astype(float), has_constant='add')
    y = df[outcome].astype(float)
    if not np.isfinite(X.to_numpy()).all() or not np.isfinite(y).all():
        row['reason'] = 'missing_values_no_imputation'; return row, [], None
    if len(df)-X.shape[1] < 6:
        row['reason'] = 'fewer_than_six_residual_df'; return row, [], None
    if np.linalg.matrix_rank(X) != X.shape[1]:
        row['reason'] = 'rank_deficient'; return row, [], None
    ols = sm.OLS(y, X).fit()
    h = ols.get_influence().hat_matrix_diag
    if np.max(h) >= 1-1e-9:
        row['reason'] = 'unit_leverage_HC3_undefined'; return row, [], None
    robust = ols.get_robustcov_results(cov_type='HC3', use_t=True)
    if not np.isfinite(robust.cov_params()).all():
        row['reason'] = 'nonfinite_HC3'; return row, [], None
    baseX = X.drop(columns=spatial_terms)
    base = sm.OLS(y, baseX).fit()
    bh = base.get_influence().hat_matrix_diag
    coef = []
    for term, beta, se, p, interval in zip(X.columns, robust.params, robust.bse, robust.pvalues, robust.conf_int()):
        coef.append(dict(model_id=model_id, family=family, term=term, beta=float(beta), se_hc3=float(se), p_hc3=float(p), ci_low=float(interval[0]), ci_high=float(interval[1])))
    c = coef[1]; sd = float(sds[spatial])
    nonlinear = bool(extra_spatial)
    row.update(status='estimated', df_resid=int(ols.df_resid), beta=c['beta'], se_hc3=c['se_hc3'], p_hc3=c['p_hc3'], ci_low=c['ci_low'], ci_high=c['ci_high'],
        spatial_sd_main=sd, effect_sd_pct=effect(c['beta'], sd) if outcome=='log_price' and not nonlinear else None,
        effect_sd_low_pct=effect(c['ci_low'], sd) if outcome=='log_price' and not nonlinear else None,
        effect_sd_high_pct=effect(c['ci_high'], sd) if outcome=='log_price' and not nonlinear else None,
        effect_5_units_pct=effect(c['beta'], 5) if outcome=='log_price' and not nonlinear else None,
        effect_sd_twd=c['beta']*sd if outcome!='log_price' else None,
        r2=float(ols.rsquared), r2_adjusted=float(ols.rsquared_adj), r2_without_ci=float(base.rsquared),
        delta_r2=float(ols.rsquared-base.rsquared), partial_r2=float((base.ssr-ols.ssr)/base.ssr),
        loo_rmse_log=float(np.sqrt(np.mean((ols.resid/(1-h))**2))) if outcome=='log_price' else None,
        loo_rmse_without_ci_log=float(np.sqrt(np.mean((base.resid/(1-bh))**2))) if outcome=='log_price' else None,
        max_leverage=float(h.max()), max_cooks_distance=float(ols.get_influence().cooks_distance[0].max()),
        n_non_spatial=sum(v in NONSPATIAL for v in controls), rare_controls='|'.join(v for v in controls if v in NONSPATIAL and df[v].nunique()==2 and min(df[v].sum(), len(df)-df[v].sum())<5),
        condition_number=float(ols.condition_number))
    if nonlinear:
        R = np.zeros((len(spatial_terms), X.shape[1]))
        for i, term in enumerate(spatial_terms): R[i, list(X.columns).index(term)] = 1
        row['p_joint_spatial_hc3'] = float(robust.f_test(R).pvalue)
    return row, coef, ols


def holm(rows, family_size):
    p = [r.get('p_hc3', 1.) if r['status']=='estimated' else 1. for r in rows]
    assert len(p) == family_size
    adjusted = multipletests(p, alpha=.05, method='holm')[1]
    for r, value in zip(rows, adjusted): r['p_holm_family'] = float(value)


def summarize_grid(rows):
    valid = [r for r in rows if r['status']=='estimated']
    a = pd.DataFrame(valid)
    def one(g):
        return dict(n=len(g), negative=int((g.beta<0).sum()), negative_share=float((g.beta<0).mean()),
            p_below_05=int((g.p_hc3<.05).sum()), holm_below_05=int((g.p_holm_family<.05).sum()),
            effect_sd_min=float(g.effect_sd_pct.min()), effect_sd_median=float(g.effect_sd_pct.median()), effect_sd_max=float(g.effect_sd_pct.max()),
            partial_r2_min=float(g.partial_r2.min()), partial_r2_median=float(g.partial_r2.median()), partial_r2_max=float(g.partial_r2.max()),
            loo_improved=int((g.loo_rmse_log < g.loo_rmse_without_ci_log).sum()))
    return dict(candidate_models=len(rows), estimated=len(valid), skipped=len(rows)-len(valid), total=one(a),
        by_competition=[dict(spatial=key, label=COMPETITION[key][0], **one(g)) for key,g in a.groupby('spatial',sort=False)],
        by_hotel=[dict(hotel=key, label=HOTELS[key][0], **one(g)) for key,g in a.groupby('hotel',sort=False)])


def diagnostics(df, sds, main_fit):
    rows, coef = [], []
    for kind, groups in [('shop', [(sid, [sid]) for sid in df.shop_id]), ('district', [(d, df.loc[df.district==d,'shop_id'].tolist()) for d in df.district.unique()])]:
        for value, ids in groups:
            sub = df[~df.shop_id.isin(ids)]
            r, c, _ = fit_model(sub, 'c_tri_5', MAIN_CONTROLS, sds, f'omit_{kind}_{value}', 'omission')
            r.update(omission_type=kind, omitted=value, omitted_n=len(ids))
            rows.append(r); coef.extend(c)
    infl = main_fit.get_influence()
    influence = df[['shop_id','name','district','reference_price_twd','c_tri_5']].copy()
    influence['residual_log'] = main_fit.resid
    influence['fitted_log'] = main_fit.fittedvalues
    influence['leverage'] = infl.hat_matrix_diag
    influence['cooks_distance'] = infl.cooks_distance[0]
    influence['studentized_residual'] = infl.resid_studentized_external
    loo = pd.DataFrame([r for r in rows if r['omission_type']=='shop']).set_index('omitted')
    influence['beta_ci_without_shop'] = influence.shop_id.map(loo.beta)
    influence['effect_sd_pct_without_shop'] = influence.shop_id.map(loo.effect_sd_pct)
    # Leave one entire district out: evaluate both models on exactly the same held-out shops.
    cv = []
    for district in df.district.unique():
        test = df[df.district==district]; train = df[df.district!=district]
        for ci in (False, True):
            features = (['c_tri_5'] if ci else []) + MAIN_CONTROLS
            X = sm.add_constant(train[features], has_constant='add')
            assert np.linalg.matrix_rank(X)==X.shape[1]
            fit = sm.OLS(train.log_price, X).fit()
            pred = fit.predict(sm.add_constant(test[features], has_constant='add'))
            for ix, value in pred.items():
                cv.append(dict(shop_id=df.loc[ix,'shop_id'], district=district, model='with_ci' if ci else 'without_ci', observed_log=float(df.loc[ix,'log_price']), predicted_log=float(value), squared_error=float((df.loc[ix,'log_price']-value)**2)))
    cvdf = pd.DataFrame(cv)
    cvsummary = {key:dict(rmse_log=float(np.sqrt(g.squared_error.mean())), mae_log=float(np.abs(g.observed_log-g.predicted_log).mean()), n=len(g)) for key,g in cvdf.groupby('model')}
    # Symmetric binary k-nearest-neighbor matrix; no row standardization.
    xy = df[['lat','lng']].to_numpy(dtype=float)
    radians = np.radians(xy)
    lat, lng = radians[:,0], radians[:,1]
    a = np.sin((lat[:,None]-lat[None,:])/2)**2 + np.cos(lat[:,None])*np.cos(lat[None,:])*np.sin((lng[:,None]-lng[None,:])/2)**2
    distances = 2*6371000*np.arcsin(np.sqrt(np.clip(a,0,1)))
    np.fill_diagonal(distances, np.inf)
    W = np.zeros((len(df),len(df)))
    for i in range(len(df)): W[i, np.argsort(distances[i],kind='stable')[:4]] = 1
    W = np.maximum(W,W.T)
    z = np.asarray(main_fit.resid); z = z-z.mean()
    observed = float(len(z)/W.sum()*(z@W@z)/(z@z))
    rng = np.random.default_rng(SEED)
    permuted = np.asarray([rng.permutation(z) for _ in range(9999)])
    Is = len(z)/W.sum()*np.einsum('bi,ij,bj->b',permuted,W,permuted)/(z@z)
    center = -1/(len(z)-1)
    p = float((1+np.sum(np.abs(Is-center)>=abs(observed-center)))/(len(Is)+1))
    moran = dict(I=observed, null_expectation=center, p_two_sided_permutation=p, permutations=9999, seed=SEED, weights='binary symmetric union of directed 4-nearest shops by geodesic distance; diagonal zero', caveat='Exploratory residual diagnostic. Does not adjust HC3 intervals or prove spatial independence.')
    X = main_fit.model.exog
    vif = [dict(term=name, vif=float(variance_inflation_factor(X,i))) for i,name in enumerate(main_fit.model.exog_names) if name!='const']
    return rows, coef, influence, cvdf, dict(spatial_cv=cvsummary, moran=moran, vif=vif)


def compare_blocks(df, main_fit):
    """Post-result, explicitly disclosed explanatory comparison requested by the user."""
    blocks={'competition':('Concurrence',['c_tri_5']),
        'demand':('Demande locale',['log_density','h_tri_10']),
        'positioning':('Note et largeur du menu',['rating','log_menu'])}
    rows=[];full=main_fit;robust=full.get_robustcov_results(cov_type='HC3',use_t=True)
    full_rmse=float(np.sqrt(np.mean((full.resid/(1-full.get_influence().hat_matrix_diag))**2)))
    for key,(label,removed) in blocks.items():
        columns=[c for c in full.model.exog_names if c not in removed and c!='const']
        reduced=sm.OLS(df.log_price,sm.add_constant(df[columns],has_constant='add')).fit()
        R=np.zeros((len(removed),len(full.model.exog_names)))
        for i,term in enumerate(removed):R[i,full.model.exog_names.index(term)]=1
        rows.append(dict(block=key,label=label,removed_variables='|'.join(removed),n=len(df),full_r2=float(full.rsquared),
            reduced_r2=float(reduced.rsquared),delta_r2=float(full.rsquared-reduced.rsquared),
            partial_r2=float((reduced.ssr-full.ssr)/reduced.ssr),p_joint_hc3=float(robust.f_test(R).pvalue),
            loo_rmse_full_log=full_rmse,loo_rmse_without_block_log=float(np.sqrt(np.mean((reduced.resid/(1-reduced.get_influence().hat_matrix_diag))**2)))))
    for row,p in zip(rows,multipletests([r['p_joint_hc3'] for r in rows],method='holm')[1]):row['p_holm_three_blocks']=float(p)
    return rows


def figures(out, df, grid, named, info, influence):
    import matplotlib
    matplotlib.use('Agg')
    import matplotlib.pyplot as plt
    from matplotlib.ticker import PercentFormatter
    plt.rcParams.update({'font.family':'DejaVu Sans','font.size':10,'axes.spines.top':False,'axes.spines.right':False,'svg.fonttype':'none','savefig.bbox':'tight'})
    def save(fig, name):
        for ext in ['svg','png','pdf']: fig.savefig(out/(name+'.'+ext),dpi=180)
        plt.close(fig)
    colors={'direct_order':'#285b4b','merchant_menu':'#a17a3f'}
    fig, ax = plt.subplots(figsize=(8.8,4.6))
    for ch,label in [('direct_order','Commande directe (20)'),('merchant_menu','Menu propre du marchand (1)')]:
        a=df[df.reference_channel==ch];ax.scatter(a.c_tri_5,a.reference_price_twd,label=label,s=42,c=colors[ch],edgecolor='white',linewidth=.7)
    ax.set(xlabel='Concurrence triangulaire à 5 min à pied',ylabel='Prix affiché du gobelet (TWD)',title=f'{len(df)} prix directs · boutiques présumées indépendantes · septembre 2026')
    ax.legend(frameon=False); ax.grid(axis='y',alpha=.15);save(fig,'econometric_prices')
    selected=[r for r in named if r['model_id'] in ['M0','M1','M2','M3','M4','M5','M6']]
    fig,ax=plt.subplots(figsize=(9,4.5));ys=np.arange(len(selected))[::-1]
    for y,r in zip(ys,selected):
        ax.errorbar(r['effect_sd_pct'],y,xerr=[[r['effect_sd_pct']-r['effect_sd_low_pct']],[r['effect_sd_high_pct']-r['effect_sd_pct']]],fmt='o',color='#285b4b',capsize=3)
    ax.set_yticks(ys,[r['label'] for r in selected]);ax.axvline(0,color='#999',lw=1)
    ax.set(xlabel='Écart de prix associé à +1 écart-type de concurrence (%)',title='Séquence de modèles · intervalles HC3 à 95 %, ponctuels')
    ax.grid(axis='x',alpha=.15);save(fig,'econometric_models')
    g=pd.DataFrame([r for r in grid if r['status']=='estimated']).sort_values(['effect_sd_pct','model_id'])
    fig,axs=plt.subplots(2,1,figsize=(9,6.7),gridspec_kw={'height_ratios':[2.3,1]})
    x=np.arange(1,len(g)+1);axs[0].fill_between(x,g.effect_sd_low_pct,g.effect_sd_high_pct,color='#e2e9e3',label='Intervalles HC3 à 95 %, non simultanés')
    axs[0].plot(x,g.effect_sd_pct,color='#285b4b',lw=1.5,label='Estimation par spécification')
    axs[0].axhline(0,color='#888',lw=.8);axs[0].set(xlabel='Spécifications classées par estimation, sans sélection',ylabel='Écart associé (%)',title=f'{len(g):,} modèles estimés · même échantillon de {len(df)} prix directs'.replace(',',' '));axs[0].legend(frameon=False,fontsize=8)
    data=[g.loc[g.spatial==k,'effect_sd_pct'] for k in COMPETITION]
    axs[1].boxplot(data,orientation='horizontal',showfliers=False,tick_labels=[v[0] for v in COMPETITION.values()],patch_artist=True,boxprops={'facecolor':'#e2e9e3'},medianprops={'color':'#285b4b'})
    axs[1].set(xlabel='Distribution des estimations par indice (%)')
    axs[1].axvline(0,color='#888',lw=.8);axs[1].tick_params(axis='y',labelsize=8)
    fig.tight_layout();save(fig,'econometric_specifications')
    fig, ax = plt.subplots(figsize=(8.5,4.3))
    ordered=influence.sort_values('effect_sd_pct_without_shop')
    ax.scatter(np.arange(len(ordered)),ordered.effect_sd_pct_without_shop,color='#285b4b',s=23)
    main=next(r for r in named if r['model_id']=='M3');ax.axhline(main['effect_sd_pct'],color='#a17a3f',label=f'Modèle avec les {len(df)} boutiques')
    ax.axhline(0,color='#888',lw=.8);ax.set(xlabel='Omission d’une boutique à la fois, classée par résultat',ylabel='Écart associé à +1 écart-type (%)',title='Influence des observations sur le modèle de référence')
    ax.legend(frameon=False,fontsize=9);save(fig,'econometric_influence')


def workbook(out, df, variables, grid, context_grid, named, coef, influence, cv, omissions, info, all_df):
    from openpyxl.styles import Font, PatternFill
    metadata = pd.DataFrame([{'item':k,'value':json.dumps(clean(v),ensure_ascii=False) if isinstance(v,(list,dict)) else v} for k,v in info.items() if k not in ['grid_summary','context_summary','models','diagnostics']])
    tables={'Lire avant analyse':metadata,'Prix directs principaux':df,'Dossiers complets tous prix':all_df,'Dictionnaire':variables,'Modeles presentes':pd.DataFrame(named),'Apport des groupes':pd.DataFrame(info['block_comparison']),'Grille 7632':pd.DataFrame(grid),'Contexte 720':pd.DataFrame(context_grid),'Tous les coefficients':pd.DataFrame(coef),'Influence boutiques':influence,'Omissions':pd.DataFrame(omissions),'Validation districts':cv,'Resume indices':pd.DataFrame(info['grid_summary']['by_competition']),'Resume hotels':pd.DataFrame(info['grid_summary']['by_hotel'])}
    with pd.ExcelWriter(out/'taipei_econometric_study.xlsx',engine='openpyxl') as writer:
        for sheet,table in tables.items():
            # Excel treats leading '=' as formulas; extracted names remain text.
            safe=table.map(lambda x: "'"+x if isinstance(x,str) and x.startswith(('=','+','-','@')) else x)
            safe.to_excel(writer,sheet_name=sheet,index=False)
            ws=writer.sheets[sheet];ws.freeze_panes='A2';ws.auto_filter.ref=ws.dimensions
            for c in ws[1]:c.font=Font(bold=True,color='FFFFFF');c.fill=PatternFill('solid',fgColor='285B4B')
            for column in ws.iter_cols(min_row=1,max_row=1):
                col=column[0];ws.column_dimensions[col.column_letter].width=min(42,max(14,len(str(col.value))+2))


def run(input_dir, out, make_figures=True, make_workbook=True):
    out.mkdir(parents=True,exist_ok=True)
    all_df, variables=prepare(input_dir)
    df=all_df[all_df.reference_channel.isin(DIRECT_CHANNELS)].copy()
    assert len(df) > 0 and df.shop_id.nunique()==len(df)
    assert not df.delivery_price.any()
    # Availability and rare categories are described on the actual main sample.
    for ix, item in variables.iterrows():
        values=df[item.variable]; binary=set(values.dropna().unique()).issubset({0,1,False,True})
        variables.loc[ix,'n_observed']=int(values.notna().sum())
        variables.loc[ix,'n_unique']=int(values.nunique())
        variables.loc[ix,'positive_cases']=int(values.sum()) if binary else np.nan
        variables.loc[ix,'rare_binary']=bool(binary and values.nunique()==2 and min(values.sum(),len(values)-values.sum())<5)
    sds={k:float(df[k].std(ddof=1)) for k in [*COMPETITION,'log_ci_5']}
    grid,context_grid,named,allcoef=[],[],[],[]
    subsets=[combo for size in range(3) for combo in itertools.combinations(NONSPATIAL,size)]
    print(f'Frozen inputs: {len(all_df)} complete dossiers; {len(df)} direct merchant prices. 9 competition / 7 hotel measures; 14 non-spatial variables.',flush=True)
    for spatial in COMPETITION:
        for hotel in HOTELS:
            for combo in subsets:
                controls=['log_density']+([] if hotel=='none' else [hotel])+list(combo)
                r,c,_=fit_model(df,spatial,controls,sds,f'G{len(grid)+1:04d}','grid')
                grid.append(r);allcoef.extend(c)
    assert len(grid)==7632
    holm(grid,7632)
    print('7,632 specifications processed; computing contextual extensions.',flush=True)
    for spatial in COMPETITION:
        for hotel in HOTELS:
            for context in CONTEXT:
                controls=['log_density']+([] if hotel=='none' else [hotel])+['rating','log_menu',context]
                r,c,_=fit_model(df,spatial,controls,sds,f'Z{len(context_grid)+1:03d}','context_grid')
                context_grid.append(r);allcoef.extend(c)
    assert len(context_grid)==720
    holm(context_grid,720)
    for r,p in zip(grid+context_grid,multipletests([r.get('p_hc3',1.) if r['status']=='estimated' else 1. for r in grid+context_grid],method='holm')[1]):r['p_holm_both_grids']=float(p)
    fits={}
    def named_fit(model_id,label,controls,sub=None,spatial='c_tri_5',outcome='log_price',extra_spatial=()):
        r,c,f=fit_model(df if sub is None else sub,spatial,controls,sds,model_id,'named',label,outcome,extra_spatial)
        named.append(r);allcoef.extend(c);fits[model_id]=f
        return r
    named_fit('M0','Concurrence seule',[])
    named_fit('M1','+ densité résidentielle',['log_density'])
    named_fit('M2','+ exposition hôtelière',['log_density','h_tri_10'])
    named_fit('M3','Référence : + note et menu',MAIN_CONTROLS)
    named_fit('M4','Référence + avis',MAIN_CONTROLS+['log_reviews'])
    named_fit('M5','Référence + horaires',MAIN_CONTROLS+['hours_per_week'])
    named_fit('M6','Référence + avis et horaires',MAIN_CONTROLS+['log_reviews','hours_per_week'])
    named_fit('S_transactional','Commande directe seulement',MAIN_CONTROLS,all_df[all_df.reference_channel=='direct_order'])
    named_fit('S_no_past_network','Sans les trois anciens réseaux non confirmés',MAIN_CONTROLS,df[~df.shop_id.isin(['places_yuanlu','fp_rb6a','fp_f5o9'])])
    named_fit('S_pooled','Tous prix · sans indicateur de canal',MAIN_CONTROLS,all_df)
    named_fit('S_pooled_channel','Tous prix · sensibilité de mesure du canal',MAIN_CONTROLS+['delivery_price'],all_df)
    named_fit('S_delivery_demand','Livraison seule · effectif insuffisant',['log_density','h_tri_10'],all_df[all_df.delivery_price==1])
    named_fit('S_black','Sans Earl Grey',MAIN_CONTROLS,df[df.earl_grey==0])
    named_fit('S_rating_source','Contrôle de source de note',MAIN_CONTROLS+['rating_nidin'])
    named_fit('S_ML','Formats M/L · demande et taille',['log_density','h_tri_10','size_L'],df[df.reference_size_label.isin(['M','L'])])
    for variable,suffix in [('h_capacity_10','capacity'),('h_rooms_10','rooms')]:
        sub=df[df[variable].notna()]
        named_fit('S_'+suffix+'_matched','Référence · cohorte '+suffix,MAIN_CONTROLS,sub)
        named_fit('S_'+suffix,'Hôtels pondérés · '+suffix,[variable if v=='h_tri_10' else v for v in MAIN_CONTROLS],sub)
    named_fit('S_levels','Prix en TWD · référence',MAIN_CONTROLS,outcome='reference_price_twd')
    named_fit('S_logCI','Concurrence ln(1+CI)',MAIN_CONTROLS,spatial='log_ci_5')
    named_fit('S_quadratic','Relation quadratique',MAIN_CONTROLS,extra_spatial=('ci_squared',))
    named_fit('S_interaction','Tous prix · interaction concurrence × canal',MAIN_CONTROLS+['delivery_price'],all_df,extra_spatial=('ci_delivery',))
    holm(named, len(named))  # Separate conservative family for all named sensitivities.
    omissions,omission_coef,influence,cv,diag=diagnostics(df,sds,fits['M3']);allcoef.extend(omission_coef)
    main=next(r for r in named if r['model_id']=='M3')
    main_grid=next(r for r in grid if r['spatial']=='c_tri_5' and r['controls']=='|'.join(MAIN_CONTROLS))
    main['grid_model_id']=main_grid['model_id'];main['p_holm_grid']=main_grid['p_holm_family']
    info=dict(status='exploratory_association',build_id='taipei-direct-prices-20260913-v9',computed_at=datetime.now(timezone.utc).isoformat(),
        research_question='To what extent does spatial competition contribute to retail prices in the bubble tea market in Taipei in 2026?',
        data_wave='2026-09-13',source_build='taipei-direct-prices-20260913-v9',n=len(df),n_complete_dossiers=len(all_df),channel_counts=df.reference_channel.value_counts().to_dict(),all_channel_counts=all_df.reference_channel.value_counts().to_dict(),sample_ids=df.shop_id.tolist(),secondary_sample_ids=all_df.shop_id.tolist(),primary_price_scope='Merchant direct order or merchant posted menu; no delivery-platform price in the main sample; counter parity and equal volume not established',documented_volumes=int(df.reference_volume_ml.notna().sum()),
        input_sha256={name:hashlib.sha256((input_dir/name).read_bytes()).hexdigest() for name in INPUTS},
        versions={'python':platform.python_version(),'numpy':np.__version__,'pandas':pd.__version__,'scipy':scipy.__version__,'statsmodels':statsmodels.__version__},
        inference='OLS; HC3; Student reference with n-k degrees of freedom; pointwise intervals; exploratory, not causal; HC3 does not correct spatial dependence',
        n_imputed_cells=0,main_model_id='M3',primary_competition='c_tri_5',primary_hotel='h_tri_10',spatial_sd_main=sds,
        competition_definitions={k:dict(label=v[0],formula=v[1]) for k,v in COMPETITION.items()},hotel_definitions={k:dict(label=v[0],formula=v[1]) for k,v in HOTELS.items()},
        non_spatial_definitions={k:dict(label=v[0],source=v[1],definition=v[2]) for k,v in NONSPATIAL.items()},context_definitions={k:dict(label=v[0],source=v[1],definition=v[2]) for k,v in CONTEXT.items()},
        grid_summary=summarize_grid(grid),context_summary=summarize_grid(context_grid),models=named,diagnostics=diag,
        omission_summary={kind:dict(n=len(a),effect_sd_min=min(r['effect_sd_pct'] for r in a),effect_sd_max=max(r['effect_sd_pct'] for r in a),positive=sum(r['beta']>0 for r in a),negative=sum(r['beta']<0 for r in a)) for kind in ['shop','district'] if (a:=[r for r in omissions if r['omission_type']==kind and r['status']=='estimated'])})
    selected_cols=['reference_price_twd','c_tri_5','c_count_5','h_tri_10_raw','h_count_10_raw','rating','review_count','hours_per_week','menu_n_products','village_population_density_km2']
    descriptive=df[selected_cols].describe(percentiles=[.25,.5,.75]).T.rename_axis('variable').reset_index()
    by_channel=df.groupby('reference_channel').agg(n=('shop_id','size'),mean_price=('reference_price_twd','mean'),median_price=('reference_price_twd','median'),min_price=('reference_price_twd','min'),max_price=('reference_price_twd','max')).reset_index()
    info['descriptives']=descriptive.to_dict('records');info['descriptives_by_channel']=by_channel.to_dict('records')
    info['block_comparison']=compare_blocks(df,fits['M3'])
    dataframes={'econometric_data.csv':df,'econometric_pooled_data.csv':all_df,'econometric_variables.csv':variables,'econometric_specifications.csv':pd.DataFrame(grid),'econometric_context_specifications.csv':pd.DataFrame(context_grid),'econometric_models.csv':pd.DataFrame(named),'econometric_coefficients.csv':pd.DataFrame(allcoef),'econometric_influence.csv':influence,'econometric_omissions.csv':pd.DataFrame(omissions),'econometric_spatial_cv.csv':cv,'econometric_descriptives.csv':descriptive}
    for name,frame in dataframes.items():frame.to_csv(out/name,index=False,encoding='utf-8-sig')
    pd.DataFrame(info['block_comparison']).to_csv(out/'econometric_blocks.csv',index=False,encoding='utf-8-sig')
    info['variable_inventory']=variables.to_dict('records')
    info['amendment_sha256']=hashlib.sha256((input_dir/'ECONOMETRIC_AMENDMENT.md').read_bytes()).hexdigest()
    info['sample_policy']='All complete presumed-independent dossiers with direct merchant prices; inclusion does not depend on the sign, price level or p-value.'
    save_json(out/'econometric_study.json',info)
    save_json(out/'econometric_input_manifest.json',{'input_sha256':info['input_sha256'],'amendment_sha256':info['amendment_sha256'],'main_sample_ids':info['sample_ids'],'secondary_sample_ids':info['secondary_sample_ids'],'versions':info['versions']})

    # Small records used by the interactive module; full coefficients remain in CSV/XLSX.
    save_json(out/'econometric_explorer.json',[{k:r.get(k) for k in ['model_id','spatial','hotel','controls','status','reason','beta','p_hc3','p_holm_family','effect_sd_pct','effect_sd_low_pct','effect_sd_high_pct','partial_r2','delta_r2','loo_rmse_log','loo_rmse_without_ci_log','rare_controls']} for r in grid])
    if make_figures:figures(out,df,grid,named,info,influence)
    if make_workbook:workbook(out,df,variables,grid,context_grid,named,allcoef,influence,cv,omissions,info,all_df)
    print(json.dumps(clean({'main':main,'grid':info['grid_summary']['total'],'context':info['context_summary']['total'],'diagnostics':diag,'omissions':info['omission_summary']}),ensure_ascii=False,indent=2),flush=True)
    return info


if __name__=='__main__':
    parser=argparse.ArgumentParser(description=__doc__)
    parser.add_argument('--input-dir',type=Path,default=Path(__file__).resolve().parent)
    parser.add_argument('--output-dir',type=Path,default=Path('econometrics_reproduced'))
    parser.add_argument('--no-figures',action='store_true')
    parser.add_argument('--no-workbook',action='store_true')
    args=parser.parse_args()
    run(args.input_dir,args.output_dir,not args.no_figures,not args.no_workbook)