"""
01_wb_data_acquisition.py
World Bank数据获取模块 - 资源重配置效率(RE)指数

获取构建RE指数所需的World Bank Worldwide Governance Indicators (WGI)
和其他相关指标，覆盖G20国家2000-2023年。
"""

import pandas as pd
import numpy as np
import requests
import json
import time
import os
from pathlib import Path

OUTPUT_DIR = Path("/home/user/re_index/raw")
OUTPUT_DIR.mkdir(parents=True, exist_ok=True)

# G20国家 ISO3代码（19个主权国家 + 欧盟成员德法意分别处理；欧盟整体单列）
G20 = {
    "ARG": "Argentina", "AUS": "Australia", "BRA": "Brazil", "CAN": "Canada",
    "CHN": "China", "FRA": "France", "DEU": "Germany", "IND": "India",
    "IDN": "Indonesia", "ITA": "Italy", "JPN": "Japan", "KOR": "Korea, Rep.",
    "MEX": "Mexico", "RUS": "Russia", "SAU": "Saudi Arabia", "ZAF": "South Africa",
    "TUR": "Turkey", "GBR": "United Kingdom", "USA": "United States",
    "EUU": "European Union"
}

# 核心WGI指标 + 其他WB指标
INDICATORS = {
    # WGI - 治理效能 (related to F: bureaucratic friction)
    "GE.EST": "Government Effectiveness",
    "RQ.EST": "Regulatory Quality",
    "CC.EST": "Control of Corruption",
    "RL.EST": "Rule of Law",
    "VA.EST": "Voice and Accountability",
    "PV.EST": "Political Stability",
    # 财政响应 (V: Velocity)
    "GC.XPN.TOTL.GD.ZS": "General government final consumption expenditure (% of GDP)",
    "GC.DOD.TOTL.GD.ZS": "Central government debt (% of GDP)",
    "GC.REV.XGRT.GD.ZS": "Revenue, excluding grants (% of GDP)",
    # 技术转换 (C: Conversion)
    "GB.XPD.RSDV.GD.ZS": "R&D expenditure (% of GDP)",
    "IP.PAT.RESD": "Patent applications, residents",
    "IT.NET.USER.ZS": "Individuals using the Internet (% of population)",
    "TX.VAL.TECH.MF.ZS": "High-technology exports (% of mfg exports)",
    # 经济基础
    "NY.GDP.MKTP.CD": "GDP (current US$)",
    "NY.GDP.PCAP.CD": "GDP per capita (current US$)",
    "NY.GDP.MKTP.KD.ZG": "GDP growth (annual %)",
    "SP.POP.TOTL": "Population, total",
}

def fetch_indicator(code, country_codes, start=2000, end=2023):
    """从World Bank API获取单个指标"""
    countries = ";".join(country_codes)
    url = f"https://api.worldbank.org/v2/country/{countries}/indicator/{code}"
    params = {"format": "json", "date": f"{start}:{end}", "per_page": 20000}
    try:
        r = requests.get(url, params=params, timeout=60)
        if r.status_code != 200:
            print(f"  ✗ {code}: HTTP {r.status_code}")
            return None
        data = r.json()
        if len(data) < 2 or data[1] is None:
            print(f"  ✗ {code}: No data")
            return None
        rows = []
        for obs in data[1]:
            rows.append({
                "iso3": obs["countryiso3code"],
                "country": obs["country"]["value"],
                "year": int(obs["date"]),
                "value": obs["value"],
                "indicator": code
            })
        df = pd.DataFrame(rows)
        return df
    except Exception as e:
        print(f"  ✗ {code}: {e}")
        return None

def main():
    print("=" * 70)
    print("World Bank Data Acquisition for RE Index")
    print("=" * 70)
    print(f"Countries: {len(G20)} (G20)")
    print(f"Indicators: {len(INDICATORS)}")
    print(f"Period: 2000-2023")
    print()
    
    all_data = []
    log = []
    for code, name in INDICATORS.items():
        print(f"Fetching {code} ({name[:50]})...")
        df = fetch_indicator(code, list(G20.keys()))
        if df is not None and len(df) > 0:
            df["indicator_name"] = name
            all_data.append(df)
            n_obs = df["value"].notna().sum()
            n_countries = df.loc[df["value"].notna(), "iso3"].nunique()
            print(f"  ✓ {len(df)} rows, {n_obs} non-null, {n_countries} countries")
            log.append({"code": code, "name": name, "rows": len(df), 
                       "non_null": int(n_obs), "countries": int(n_countries),
                       "status": "OK"})
        else:
            log.append({"code": code, "name": name, "rows": 0,
                       "non_null": 0, "countries": 0, "status": "FAILED"})
        time.sleep(0.2)  # gentle rate-limiting
    
    if all_data:
        full = pd.concat(all_data, ignore_index=True)
        # Long format
        out_long = OUTPUT_DIR / "wb_data_long.csv"
        full.to_csv(out_long, index=False)
        print(f"\n✓ Long format saved: {out_long} ({len(full)} rows)")
        
        # Wide format (pivot)
        wide = full.pivot_table(index=["iso3","country","year"], 
                                columns="indicator", values="value").reset_index()
        out_wide = OUTPUT_DIR / "wb_data_wide.csv"
        wide.to_csv(out_wide, index=False)
        print(f"✓ Wide format saved: {out_wide} ({len(wide)} rows, {len(wide.columns)} cols)")
    
    log_df = pd.DataFrame(log)
    log_df.to_csv(OUTPUT_DIR / "wb_acquisition_log.csv", index=False)
    print(f"\n=== Acquisition Summary ===")
    print(log_df.to_string(index=False))

if __name__ == "__main__":
    main()
"""
01b_wgi_acquisition.py
WGI数据专项获取（使用source=3和正确指标代码）
"""
import pandas as pd
import requests
import time
from pathlib import Path

OUTPUT_DIR = Path("/home/user/re_index/raw")

G20 = ["ARG","AUS","BRA","CAN","CHN","FRA","DEU","IND","IDN","ITA",
       "JPN","KOR","MEX","RUS","SAU","ZAF","TUR","GBR","USA","EUU"]

WGI = {
    "GOV_WGI_GE.EST": "Government Effectiveness",
    "GOV_WGI_RQ.EST": "Regulatory Quality",
    "GOV_WGI_CC.EST": "Control of Corruption",
    "GOV_WGI_RL.EST": "Rule of Law",
    "GOV_WGI_VA.EST": "Voice and Accountability",
    "GOV_WGI_PV.EST": "Political Stability",
}

def fetch_wgi(code, countries, start=2000, end=2023):
    iso = ";".join(countries)
    # WGI is in source=3
    url = f"https://api.worldbank.org/v2/country/{iso}/indicator/{code}"
    params = {"format":"json","date":f"{start}:{end}","per_page":20000,"source":3}
    r = requests.get(url, params=params, timeout=60)
    if r.status_code != 200:
        return None
    data = r.json()
    if len(data) < 2 or data[1] is None:
        return None
    rows = []
    for obs in data[1]:
        rows.append({
            "iso3": obs["countryiso3code"],
            "country": obs["country"]["value"],
            "year": int(obs["date"]),
            "value": obs["value"],
            "indicator": code.replace("GOV_WGI_","")
        })
    return pd.DataFrame(rows)

all_data = []
for code, name in WGI.items():
    print(f"Fetching {code}...")
    df = fetch_wgi(code, G20)
    if df is not None:
        df["indicator_name"] = name
        all_data.append(df)
        print(f"  ✓ {len(df)} rows, non-null={df['value'].notna().sum()}, countries={df.loc[df['value'].notna(),'iso3'].nunique()}")
    time.sleep(0.3)

if all_data:
    full = pd.concat(all_data, ignore_index=True)
    full.to_csv(OUTPUT_DIR / "wgi_data_long.csv", index=False)
    print(f"\n✓ WGI data saved: {len(full)} rows")
    wide = full.pivot_table(index=["iso3","country","year"], columns="indicator", values="value").reset_index()
    wide.to_csv(OUTPUT_DIR / "wgi_data_wide.csv", index=False)
    print(f"✓ WGI wide saved: {wide.shape}")
    # Quick check
    print("\nLatest 2022 GE.EST sample:")
    sample = full[(full.indicator=="GE.EST")&(full.year==2022)&(full.value.notna())][["iso3","value"]].head(20)
    print(sample.to_string(index=False))
"""
02_extract_qog_vars.py
从Quality of Government (QoG) Standard Dataset中提取RE指数相关变量
"""
import pandas as pd
import numpy as np
from pathlib import Path

RAW = Path("/home/user/re_index/raw")
PROC = Path("/home/user/re_index/processed")
PROC.mkdir(exist_ok=True)

# Read header first to scan available variables
print("Loading QoG header to identify variables...")
header = pd.read_csv(RAW/"qog_std_ts.csv", nrows=0)
cols = list(header.columns)
print(f"Total variables: {len(cols)}")

# Key variable prefixes we care about
KEY_PATTERNS = {
    'V-Dem (vdem)':       ['vdem_'],
    'ICRG':               ['icrg_'],
    'Fraser EFW':         ['fi_'],
    'WGI':                ['wbgi_'],
    'Bertelsmann (bti/sgi)': ['bti_','sgi_'],
    'Hanson-Sigman SCI':  ['hsl_'],
    'World Bank':         ['wdi_'],
    'IMF':                ['imf_'],
    'Heritage':           ['hf_'],
    'OECD':               ['oecd_','pmr_'],
    'WIPO/Innovation':    ['gii_','wipo_'],
}

print("\n=== Variable categories ===")
matched = {}
for label, prefixes in KEY_PATTERNS.items():
    found = [c for c in cols if any(c.startswith(p) for p in prefixes)]
    matched[label] = found
    print(f"{label}: {len(found)} variables (sample: {found[:3]})")

# Specifically pick RE-relevant variables
RE_VARS = {
    # === V (Velocity/Reaction): WGI GE, Fiscal indicators ===
    'wbgi_gee':      'WGI Government Effectiveness Estimate',  
    'wbgi_rqe':      'WGI Regulatory Quality Estimate',
    'wbgi_rle':      'WGI Rule of Law Estimate',
    'wbgi_pse':      'WGI Political Stability Estimate',
    'wbgi_vae':      'WGI Voice & Accountability Estimate',
    'wbgi_cce':      'WGI Control of Corruption Estimate',
    # === V: Fiscal responsiveness ===
    'wdi_gge':       'WDI General government final consumption expenditure (% GDP)',
    'wdi_taxrev':    'WDI Tax revenue (% GDP)',
    # === C: Tech conversion ===
    'wdi_internet':  'WDI Individuals using Internet (% population)',
    'wdi_expgnfsh':  'WDI High-tech exports',  
    'wdi_rnd':       'WDI R&D expenditure (% GDP)',
    # === F: Bureaucratic friction ===
    'icrg_qog':      'ICRG Quality of Government index',  
    # === Fraser: Regulation ===
    'fi_regulation': 'Fraser Regulatory Freedom',
    'fi_legalsystem':'Fraser Legal System',
    'fi_index':      'Fraser Economic Freedom overall',
    # === V-Dem: Bureaucratic ===
    'vdem_corr':     'V-Dem political corruption index',
    'vdem_libdem':   'V-Dem liberal democracy index',
    'vdem_egaldem':  'V-Dem egalitarian democracy',
    # === Hanson-Sigman state capacity ===
    'hsl_scfac':     'Hanson-Sigman State Capacity',
    # === Heritage ===
    'hf_score':      'Heritage Economic Freedom Score',
    'hf_busfree':    'Heritage Business Freedom',
}

# Discover actual variable names in QoG
print("\n=== Searching for actual variable names ===")
ID_COLS = ['cname','year','ccodealp','ccode']
selected = list(ID_COLS)
found_map = {}
for target_key, desc in RE_VARS.items():
    # Find direct match or pattern
    matches = [c for c in cols if c == target_key or c.startswith(target_key.replace('_',''))]
    if matches:
        for m in matches[:2]:
            selected.append(m)
            found_map[m] = desc
        print(f"  ✓ {target_key} → {matches[:2]}")
    else:
        # Try fuzzy
        parts = target_key.split('_')
        fuzzy = [c for c in cols if all(p.lower() in c.lower() for p in parts)]
        if fuzzy:
            selected.append(fuzzy[0])
            found_map[fuzzy[0]] = desc
            print(f"  ~ {target_key} → {fuzzy[0]} (fuzzy)")
        else:
            print(f"  ✗ {target_key} not found")

selected = list(dict.fromkeys(selected))  # dedup
print(f"\nSelected columns: {len(selected)}")

# Load only selected columns
print("\nLoading filtered data...")
df = pd.read_csv(RAW/"qog_std_ts.csv", usecols=selected, low_memory=False)
print(f"Loaded: {df.shape}")

# Filter G20 countries
G20_ISO3 = ["ARG","AUS","BRA","CAN","CHN","FRA","DEU","IND","IDN","ITA",
            "JPN","KOR","MEX","RUS","SAU","ZAF","TUR","GBR","USA"]
df_g20 = df[df['ccodealp'].isin(G20_ISO3)].copy()
df_g20 = df_g20[(df_g20.year>=2000)&(df_g20.year<=2023)].sort_values(['ccodealp','year']).reset_index(drop=True)
print(f"G20 filtered: {df_g20.shape}")

# Save
df_g20.to_csv(PROC/"qog_g20_re_vars.csv", index=False)
print(f"✓ Saved: {PROC/'qog_g20_re_vars.csv'}")

# Variable mapping
mapping = pd.DataFrame([(k,v) for k,v in found_map.items()], columns=['variable','description'])
mapping.to_csv(PROC/"qog_variable_mapping.csv", index=False)
print(f"✓ Mapping saved: {PROC/'qog_variable_mapping.csv'}")

# Quick coverage check
print("\n=== Coverage check (% non-null by variable) ===")
data_vars = [c for c in df_g20.columns if c not in ID_COLS]
coverage = df_g20[data_vars].notna().mean().sort_values(ascending=False)
print(coverage.to_string())
"""
03_full_data_extraction.py
完整RE指数变量提取 - 从QoG数据集 + World Bank API补充
"""
import pandas as pd
import numpy as np
import requests
import time
from pathlib import Path

RAW = Path("/home/user/re_index/raw")
PROC = Path("/home/user/re_index/processed")

G20_ISO3 = ["ARG","AUS","BRA","CAN","CHN","FRA","DEU","IND","IDN","ITA",
            "JPN","KOR","MEX","RUS","SAU","ZAF","TUR","GBR","USA"]
G20_NAMES = {"ARG":"Argentina","AUS":"Australia","BRA":"Brazil","CAN":"Canada",
             "CHN":"China","FRA":"France","DEU":"Germany","IND":"India",
             "IDN":"Indonesia","ITA":"Italy","JPN":"Japan","KOR":"Korea, Rep.",
             "MEX":"Mexico","RUS":"Russian Federation","SAU":"Saudi Arabia",
             "ZAF":"South Africa","TUR":"Turkiye","GBR":"United Kingdom","USA":"United States"}

# ===== RE指数变量映射（最终版） =====
RE_VAR_MAP = {
    # === V (Velocity/Reaction Elasticity) - 财政响应能力 & 监管适应性 ===
    # 财政响应 (V_fiscal)
    'wdi_taxrev':     ('V_fiscal', 'Tax revenue (% GDP)', 'higher', 1.0),
    # 监管适应性 (V_reg)  
    'wbgi_rqe':       ('V_reg',    'WGI Regulatory Quality', 'higher', 1.0),
    'fi_reg':         ('V_reg',    'Fraser Regulation Freedom', 'higher', 1.0),
    'fi_legprop':     ('V_reg',    'Fraser Legal System & Property', 'higher', 1.0),
    'wbgi_rle':       ('V_reg',    'WGI Rule of Law', 'higher', 1.0),
    # 政治稳定与执行能力
    'wbgi_pve':       ('V_pol',    'WGI Political Stability', 'higher', 1.0),
    'wbgi_vae':       ('V_pol',    'WGI Voice & Accountability', 'higher', 1.0),
    
    # === C (Conversion Rate) - 技术转换速度 ===
    'wdi_internet':   ('C_tech',   'Internet users (% pop)', 'higher', 1.0),
    'vdem_academ':    ('C_inst',   'V-Dem Academic Freedom', 'higher', 1.0),
    
    # === F (Friction) - 制度摩擦 & 腐败 ===
    'wbgi_cce':       ('F_corr',   'WGI Control of Corruption', 'lower', 1.0),  # higher = less F
    'icrg_qog':       ('F_bureau', 'ICRG Quality of Government', 'lower', 1.0),  # higher = less F
    'vdem_corr':      ('F_corr',   'V-Dem Political Corruption Index', 'higher_friction', 1.0), # higher = more F
    'wbgi_gee':       ('F_bureau', 'WGI Government Effectiveness', 'lower', 1.0),  # higher = less F
    'fi_sog':         ('F_govsize','Fraser Size of Government', 'higher', 1.0),  # higher = less F
}

# ID columns
ID_COLS = ['cname','ccodealp','year']
data_vars = list(RE_VAR_MAP.keys())
all_cols = ID_COLS + data_vars

print("Loading QoG data with selected variables...")
qog = pd.read_csv(RAW/"qog_std_ts.csv", usecols=all_cols, low_memory=False)
qog = qog[qog['ccodealp'].isin(G20_ISO3)]
qog = qog[(qog.year>=2000)&(qog.year<=2023)].sort_values(['ccodealp','year']).reset_index(drop=True)
qog = qog.rename(columns={'ccodealp':'iso3','cname':'country'})
print(f"QoG G20 filtered: {qog.shape}")
print(f"Years range: {qog.year.min()}-{qog.year.max()}")
print(f"Countries: {sorted(qog.iso3.unique())}")

# ===== 从World Bank API补充R&D和高科技出口 =====
print("\n=== Fetching supplementary WB indicators ===")
def fetch_wb(code, countries=G20_ISO3):
    iso = ";".join(countries)
    url = f"https://api.worldbank.org/v2/country/{iso}/indicator/{code}"
    params = {"format":"json","date":"2000:2023","per_page":20000}
    r = requests.get(url, params=params, timeout=60)
    if r.status_code!=200 or len(r.json())<2 or r.json()[1] is None:
        return None
    rows = [{'iso3':o['countryiso3code'],'year':int(o['date']),'value':o['value']} 
            for o in r.json()[1]]
    df = pd.DataFrame(rows)
    df.columns = ['iso3','year',code.replace('.','_')]
    return df

SUPP = {
    'GB.XPD.RSDV.GD.ZS': 'C_rd',       # R&D % GDP - C: Conversion
    'TX.VAL.TECH.MF.ZS': 'C_hitech',   # High-tech exports - C: Conversion
    'GC.XPN.TOTL.GD.ZS': 'V_govexp',   # Gov expenditure - V: Fiscal
    'NY.GDP.PCAP.CD':    'gdppc',
    'NY.GDP.MKTP.KD.ZG': 'gdpgrowth',
    'SP.POP.TOTL':       'pop',
}
extra_dfs = []
for code, alias in SUPP.items():
    print(f"  Fetching {code} ({alias})...")
    df = fetch_wb(code)
    if df is not None:
        col_old = code.replace('.','_')
        df = df.rename(columns={col_old: alias})
        extra_dfs.append(df)
        n_nn = df[alias].notna().sum()
        print(f"    ✓ {len(df)} rows, non-null={n_nn}")
    time.sleep(0.25)

# Merge
master = qog.copy()
for ed in extra_dfs:
    master = master.merge(ed, on=['iso3','year'], how='left')
print(f"\nMaster shape: {master.shape}")

# Save raw merged
master.to_csv(PROC/"re_master_raw.csv", index=False)
print(f"✓ Saved: {PROC/'re_master_raw.csv'}")

# ===== Coverage report =====
print("\n=== Variable Coverage (% non-null over 19 countries × 24 years = 456 cells) ===")
total = 19*24
coverage = {}
for col in master.columns:
    if col in ['iso3','country','year']: continue
    nn = master[col].notna().sum()
    coverage[col] = (nn, round(100*nn/total,1))
cov_df = pd.DataFrame([(k,v[0],v[1]) for k,v in coverage.items()], 
                       columns=['variable','non_null','pct']).sort_values('pct', ascending=False)
print(cov_df.to_string(index=False))
cov_df.to_csv(PROC/"variable_coverage.csv", index=False)

# Variable mapping export
mapping_rows = []
for var, info in RE_VAR_MAP.items():
    mapping_rows.append({'variable':var,'dimension':info[0],'description':info[1],
                        'direction':info[2],'weight_hint':info[3],'source':'QoG'})
for code, alias in SUPP.items():
    dim = 'C_rd' if 'RSDV' in code else 'C_hitech' if 'TECH' in code else 'V_govexp' if 'XPN' in code else 'control'
    desc = {'GB.XPD.RSDV.GD.ZS':'R&D expenditure (% GDP)',
            'TX.VAL.TECH.MF.ZS':'High-tech exports (% manuf)',
            'GC.XPN.TOTL.GD.ZS':'Gov final consumption (% GDP)',
            'NY.GDP.PCAP.CD':'GDP per capita US$',
            'NY.GDP.MKTP.KD.ZG':'GDP growth %',
            'SP.POP.TOTL':'Population'}[code]
    mapping_rows.append({'variable':alias,'dimension':dim,'description':desc,
                        'direction':'higher','weight_hint':1.0,'source':'WB API'})
pd.DataFrame(mapping_rows).to_csv(PROC/"re_variable_dimension_map.csv", index=False)
print(f"\n✓ Variable-dimension mapping saved")
"""
04_re_index_construction.py
=========================================================================
资源重配置效率（RE）指数构建主程序
Resource Reallocation Efficiency (RE) Index Construction

核心方程 / Core Equation:
    RE = V · C / (1 + F)

方法学创新 / Methodological Innovations:
  [创新1] 双轨权重 (Dual-Track Weighting)：一级维度德尔菲规范权重 + 
          二级指标PCA数据驱动权重
  [创新2] 乘积-比率聚合 (Multiplicative-Ratio Aggregation)：保留RE方程
          的乘积-比率结构 + 对数变换稳健版本
  [创新3] 时变权重 (Dynamic Time-Varying Weights)：5年滚动PCA权重
  [创新4] 领域子指数 (Domain-Specific Sub-Indices)：财政/技术/治理子指数
  [创新5] 不确定性量化 (Uncertainty Quantification)：bootstrap重抽样
          构造RE指数的90%置信区间
=========================================================================
"""

import pandas as pd
import numpy as np
from pathlib import Path
from sklearn.experimental import enable_iterative_imputer  # noqa
from sklearn.impute import IterativeImputer
from sklearn.preprocessing import StandardScaler, MinMaxScaler
from sklearn.decomposition import PCA
import warnings
warnings.filterwarnings('ignore')

np.random.seed(42)

PROC = Path("/home/user/re_index/processed")
OUT  = Path("/home/user/re_index/output")
OUT.mkdir(exist_ok=True)

# ======================================================================
# Step 1: Load and prepare data
# ======================================================================
print("="*72)
print("RE INDEX CONSTRUCTION — 资源重配置效率指数构建")
print("="*72)

df = pd.read_csv(PROC/"re_master_raw.csv")
mapping = pd.read_csv(PROC/"re_variable_dimension_map.csv")
print(f"Loaded master data: {df.shape}")

# Define dimension membership (final theoretical structure)
DIM_V = {  # Velocity / Reaction Elasticity
    'wbgi_rqe':  ('higher', 'Regulatory Quality'),
    'wbgi_pve':  ('higher', 'Political Stability'),
    'wbgi_vae':  ('higher', 'Voice & Accountability'),
    'fi_reg':    ('higher', 'Fraser Regulatory Freedom'),
    'fi_legprop':('higher', 'Legal System Quality'),
    'wbgi_rle':  ('higher', 'Rule of Law'),
    'wdi_taxrev':('higher', 'Tax Revenue % GDP (fiscal capacity)'),
    'V_govexp':  ('higher', 'Gov Expenditure % GDP'),
}
DIM_C = {  # Conversion (technological)
    'wdi_internet':('higher', 'Internet Penetration'),
    'C_rd':       ('higher', 'R&D Expenditure % GDP'),
    'C_hitech':   ('higher', 'High-Tech Exports'),
    'vdem_academ':('higher', 'Academic Freedom'),
}
DIM_F = {  # Friction (bureaucratic + corruption); needs sign reversal for some
    'vdem_corr':  ('friction', 'V-Dem Political Corruption (0-1, higher=more corrupt)'),
    'wbgi_cce':   ('reverse',  'WGI Control of Corruption (reverse for F)'),
    'icrg_qog':   ('reverse',  'ICRG QoG (reverse for F)'),
    'wbgi_gee':   ('reverse',  'WGI Gov Effectiveness (reverse for F)'),
    'fi_sog':     ('reverse',  'Fraser Size of Gov (reverse for F)'),
}

ALL_VARS = list(DIM_V.keys()) + list(DIM_C.keys()) + list(DIM_F.keys())
print(f"Dimensions: V={len(DIM_V)}, C={len(DIM_C)}, F={len(DIM_F)}")
print(f"Total indicators: {len(ALL_VARS)}")

# Subset
work = df[['iso3','country','year']+ALL_VARS].copy()

# ======================================================================
# Step 2: Missing data imputation (Multiple Iterative Imputation)
# ======================================================================
print("\n" + "="*72)
print("Step 2: Missing Data Imputation (Iterative MICE-like)")
print("="*72)
miss_before = work[ALL_VARS].isna().sum().sum()
print(f"Missing cells before imputation: {miss_before} of {len(work)*len(ALL_VARS)} "
      f"({100*miss_before/(len(work)*len(ALL_VARS)):.2f}%)")

# Country-fixed-effects via pivot-wise imputation
imp = IterativeImputer(max_iter=20, random_state=42, sample_posterior=False)
work_imp = work.copy()
work_imp[ALL_VARS] = imp.fit_transform(work[ALL_VARS])
miss_after = work_imp[ALL_VARS].isna().sum().sum()
print(f"Missing cells after imputation:  {miss_after}")
work_imp.to_csv(PROC/"re_data_imputed.csv", index=False)

# ======================================================================
# Step 3: Normalization (Min-Max to [0,1])
# ======================================================================
print("\n" + "="*72)
print("Step 3: Min-Max Normalization to [0,1]")
print("="*72)
norm = work_imp[['iso3','country','year']].copy()
for v in ALL_VARS:
    x = work_imp[v].values
    rng = np.nanmax(x) - np.nanmin(x)
    if rng > 0:
        norm[v] = (x - np.nanmin(x)) / rng
    else:
        norm[v] = 0.5

# Direction handling for F dimension
F_vars = list(DIM_F.keys())
F_norm = norm[F_vars].copy()
for v, (dir_, _) in DIM_F.items():
    if dir_ == 'reverse':
        F_norm[v] = 1.0 - F_norm[v]   # so that higher = more friction
    # 'friction' means already higher = more friction
# Replace in norm
for v in F_vars:
    norm[v+'_F'] = F_norm[v]

print(f"Normalized data shape: {norm.shape}")

# ======================================================================
# Step 4: Dimension scores via PCA (innovation 1: data-driven weights)
# ======================================================================
print("\n" + "="*72)
print("Step 4: PCA-based Sub-dimension Aggregation (Innovation #1)")
print("="*72)

def pca_weight(df_dim, var_list, label):
    """Calculate PCA-based weights for a dimension."""
    X = df_dim[var_list].values
    # Standardize first for PCA
    Xs = (X - X.mean(0)) / (X.std(0) + 1e-9)
    pca = PCA(n_components=1)
    pca.fit(Xs)
    loadings = np.abs(pca.components_[0])
    weights = loadings / loadings.sum()
    explained = pca.explained_variance_ratio_[0]
    print(f"  {label}: PC1 explains {explained*100:.1f}% variance")
    for v, w in zip(var_list, weights):
        print(f"    {v}: weight = {w:.4f}")
    return weights, explained

V_vars = list(DIM_V.keys())
C_vars = list(DIM_C.keys())
F_vars_signed = [v+'_F' for v in DIM_F.keys()]

print("\n[V dimension - Velocity/Reaction Elasticity]")
w_V, ev_V = pca_weight(norm, V_vars, "V")
print("\n[C dimension - Conversion]")
w_C, ev_C = pca_weight(norm, C_vars, "C")
print("\n[F dimension - Institutional Friction]")
w_F, ev_F = pca_weight(norm, F_vars_signed, "F")

# Composite sub-scores
norm['V_score'] = np.dot(norm[V_vars].values, w_V)
norm['C_score'] = np.dot(norm[C_vars].values, w_C)
norm['F_score'] = np.dot(norm[F_vars_signed].values, w_F)

print("\n=== Dimension score statistics ===")
print(norm[['V_score','C_score','F_score']].describe().round(4))

# Equal-weight (Delphi normative) as comparison: each sub-indicator equally weighted
norm['V_score_eq'] = norm[V_vars].mean(axis=1)
norm['C_score_eq'] = norm[C_vars].mean(axis=1)
norm['F_score_eq'] = norm[F_vars_signed].mean(axis=1)

# ======================================================================
# Step 5: RE Index aggregation (Innovation #2: multiplicative-ratio)
# ======================================================================
print("\n" + "="*72)
print("Step 5: RE Index Aggregation — RE = V·C / (1+F)  (Innovation #2)")
print("="*72)

# Rescale V,C,F to (0,1] avoiding 0 to keep formula meaningful
def safe01(x, eps=1e-3):
    return np.clip(x, eps, 1.0)

norm['V_s'] = safe01(norm['V_score'])
norm['C_s'] = safe01(norm['C_score'])
norm['F_s'] = safe01(norm['F_score'])

# Main RE (PCA weights)
norm['RE'] = (norm['V_s'] * norm['C_s']) / (1.0 + norm['F_s'])

# RE_eq (equal weights as comparison)
norm['RE_eq'] = (safe01(norm['V_score_eq']) * safe01(norm['C_score_eq'])) / (1.0 + safe01(norm['F_score_eq']))

# Log-transformed robust version
norm['lnRE'] = np.log(norm['V_s']) + np.log(norm['C_s']) - np.log(1+norm['F_s'])

# Rescale RE to 0-100
re_min = norm['RE'].min(); re_max = norm['RE'].max()
norm['RE_0_100'] = 100*(norm['RE']-re_min)/(re_max-re_min)

print("\n=== RE Index (raw multiplicative) ===")
print(norm[['RE','RE_eq','lnRE','RE_0_100']].describe().round(4))

# ======================================================================
# Step 6: Innovation #3 — Time-varying PCA weights (5-year window)
# ======================================================================
print("\n" + "="*72)
print("Step 6: Dynamic Time-Varying Weights (Innovation #3)")
print("="*72)
years = sorted(norm.year.unique())
tv_weights_V = {}
tv_weights_C = {}
tv_weights_F = {}
for y in years:
    # Use 5-year window centered on y
    sub = norm[(norm.year>=y-2)&(norm.year<=y+2)]
    if len(sub) < 30:
        sub = norm  # fallback
    for vars_, store in [(V_vars,tv_weights_V),(C_vars,tv_weights_C),(F_vars_signed,tv_weights_F)]:
        X = sub[vars_].values
        Xs = (X - X.mean(0))/(X.std(0)+1e-9)
        pca = PCA(n_components=1).fit(Xs)
        ld = np.abs(pca.components_[0])
        store[y] = ld/ld.sum()

# Apply time-varying weights
norm['V_tv'] = 0.0; norm['C_tv']=0.0; norm['F_tv']=0.0
for i, row in norm.iterrows():
    y = row['year']
    norm.at[i,'V_tv'] = np.dot(row[V_vars].values, tv_weights_V[y])
    norm.at[i,'C_tv'] = np.dot(row[C_vars].values, tv_weights_C[y])
    norm.at[i,'F_tv'] = np.dot(row[F_vars_signed].values, tv_weights_F[y])

norm['RE_timevar'] = safe01(norm['V_tv'])*safe01(norm['C_tv'])/(1+safe01(norm['F_tv']))
print(f"Time-varying RE range: [{norm['RE_timevar'].min():.4f}, {norm['RE_timevar'].max():.4f}]")

# ======================================================================
# Step 7: Innovation #4 — Domain-specific sub-indices
# ======================================================================
print("\n" + "="*72)
print("Step 7: Domain-Specific Sub-Indices (Innovation #4)")
print("="*72)

# RE_fiscal: emphasizes fiscal (taxrev, V_govexp, govEffectiveness)
norm['V_fiscal'] = norm[['wdi_taxrev','V_govexp']].mean(axis=1)
norm['F_fiscal'] = (1 - norm['wbgi_gee']).clip(lower=0)
norm['RE_fiscal'] = safe01(norm['V_fiscal'])*safe01(norm['C_score'])/(1+safe01(norm['F_fiscal']))

# RE_tech: emphasizes innovation (R&D, hightech, internet, academic)
norm['C_tech_domain'] = norm[['wdi_internet','C_rd','C_hitech','vdem_academ']].mean(axis=1)
norm['RE_tech'] = safe01(norm['V_score'])*safe01(norm['C_tech_domain'])/(1+safe01(norm['F_score']))

# RE_governance: regulatory-focused
norm['V_gov'] = norm[['wbgi_rqe','wbgi_rle','fi_reg','fi_legprop']].mean(axis=1)
norm['RE_gov'] = safe01(norm['V_gov'])*safe01(norm['C_score'])/(1+safe01(norm['F_score']))

print("Domain RE statistics:")
print(norm[['RE_fiscal','RE_tech','RE_gov']].describe().round(4))

# ======================================================================
# Step 8: Innovation #5 — Bootstrap uncertainty quantification
# ======================================================================
print("\n" + "="*72)
print("Step 8: Bootstrap Uncertainty Quantification (Innovation #5)")
print("="*72)

B = 500
boot_re = np.zeros((B, len(norm)))
for b in range(B):
    # Resample variables with replacement (block-bootstrap on weights)
    # Bootstrap the weights by perturbing them with a Dirichlet noise
    wV_b = np.random.dirichlet(w_V*100 + 1)
    wC_b = np.random.dirichlet(w_C*100 + 1)
    wF_b = np.random.dirichlet(w_F*100 + 1)
    V_b = safe01(np.dot(norm[V_vars].values, wV_b))
    C_b = safe01(np.dot(norm[C_vars].values, wC_b))
    F_b = safe01(np.dot(norm[F_vars_signed].values, wF_b))
    boot_re[b,:] = V_b*C_b/(1+F_b)

norm['RE_mean'] = boot_re.mean(0)
norm['RE_lower90'] = np.percentile(boot_re, 5, axis=0)
norm['RE_upper90'] = np.percentile(boot_re, 95, axis=0)
norm['RE_se'] = boot_re.std(0)
print(f"Bootstrap done: B={B} replications")
print(f"Mean SE: {norm['RE_se'].mean():.4f}")
print(f"Average 90% CI width: {(norm['RE_upper90']-norm['RE_lower90']).mean():.4f}")

# ======================================================================
# Step 9: Save full dataset
# ======================================================================
print("\n" + "="*72)
print("Step 9: Save Final RE Index Dataset")
print("="*72)

# Save final dataset
out_cols = ['iso3','country','year'] + ALL_VARS + [v+'_F' for v in DIM_F.keys()] + \
           ['V_score','C_score','F_score','V_score_eq','C_score_eq','F_score_eq',
            'RE','RE_eq','lnRE','RE_0_100','V_tv','C_tv','F_tv','RE_timevar',
            'RE_fiscal','RE_tech','RE_gov','RE_mean','RE_lower90','RE_upper90','RE_se']
final = norm[out_cols].copy()
final.to_csv(OUT/"G20_RE_Index_2000_2023.csv", index=False)
print(f"✓ Saved: {OUT/'G20_RE_Index_2000_2023.csv'} ({final.shape})")

# Also save weights tables
pd.DataFrame({'variable':V_vars,'weight_PCA':w_V}).to_csv(OUT/"weights_V.csv", index=False)
pd.DataFrame({'variable':C_vars,'weight_PCA':w_C}).to_csv(OUT/"weights_C.csv", index=False)
pd.DataFrame({'variable':F_vars_signed,'weight_PCA':w_F}).to_csv(OUT/"weights_F.csv", index=False)

# Time-varying weights
tv_V_df = pd.DataFrame(tv_weights_V, index=V_vars).T
tv_V_df.index.name = 'year'
tv_V_df.to_csv(OUT/"time_varying_weights_V.csv")
tv_C_df = pd.DataFrame(tv_weights_C, index=C_vars).T
tv_C_df.index.name = 'year'
tv_C_df.to_csv(OUT/"time_varying_weights_C.csv")
tv_F_df = pd.DataFrame(tv_weights_F, index=F_vars_signed).T
tv_F_df.index.name = 'year'
tv_F_df.to_csv(OUT/"time_varying_weights_F.csv")

# Quick summary by country (latest year)
print("\n=== RE Index Top/Bottom (year 2022) ===")
latest = final[final.year==2022].sort_values('RE_0_100', ascending=False)
print(latest[['country','V_score','C_score','F_score','RE','RE_0_100','RE_lower90','RE_upper90']].round(3).to_string(index=False))

# Save PCA explained variance metadata
meta = pd.DataFrame({
    'dimension':['V','C','F'],
    'PC1_explained_var':[ev_V, ev_C, ev_F],
    'n_indicators':[len(V_vars),len(C_vars),len(F_vars_signed)]
})
meta.to_csv(OUT/"pca_metadata.csv", index=False)
print("\n✓ All output files saved to:", OUT)
"""
05_reliability_validity.py
信度与效度检验
Reliability and Validity Testing for RE Index
"""
import pandas as pd
import numpy as np
from pathlib import Path
from scipy import stats
import warnings
warnings.filterwarnings('ignore')

OUT = Path("/home/user/re_index/output")
PROC = Path("/home/user/re_index/processed")

df = pd.read_csv(OUT/"G20_RE_Index_2000_2023.csv")
print("="*72)
print("RELIABILITY AND VALIDITY TESTING")
print("="*72)
print(f"Data: {df.shape}")

# ======================================================================
# 1. Cronbach's alpha for each dimension
# ======================================================================
def cronbach_alpha(items):
    items = items.dropna()
    k = items.shape[1]
    var_sum = items.sum(axis=1).var(ddof=1)
    var_items = items.var(ddof=1).sum()
    return (k/(k-1))*(1 - var_items/var_sum)

V_vars = ['wbgi_rqe','wbgi_pve','wbgi_vae','fi_reg','fi_legprop','wbgi_rle','wdi_taxrev','V_govexp']
C_vars = ['wdi_internet','C_rd','C_hitech','vdem_academ']
F_vars = ['vdem_corr_F','wbgi_cce_F','icrg_qog_F','wbgi_gee_F','fi_sog_F']

# Need normalized data (re-scale 0-1 first since unit differences)
def minmax01(df_, cols):
    out = df_[cols].copy()
    for c in cols:
        x = out[c]
        rng = x.max()-x.min()
        if rng>0: out[c]=(x-x.min())/rng
    return out

V_n = minmax01(df, V_vars)
C_n = minmax01(df, C_vars)
F_n = df[F_vars].copy()  # already 0-1 normalized

print("\n=== Cronbach's α ===")
for label, items in [('V', V_n), ('C', C_n), ('F', F_n)]:
    a = cronbach_alpha(items)
    n = items.shape[1]
    print(f"  {label} (n={n} items): α = {a:.4f}  "
          f"{'(Excellent)' if a>=0.9 else '(Good)' if a>=0.8 else '(Acceptable)' if a>=0.7 else '(Questionable)'}")

# ======================================================================
# 2. Composite Reliability (CR) and AVE
# ======================================================================
def cr_ave(items, label):
    """Composite Reliability and Average Variance Extracted (factor-loading-based)."""
    # Use PC1 loadings as factor loadings
    from sklearn.decomposition import PCA
    X = items.dropna().values
    Xs = (X-X.mean(0))/(X.std(0)+1e-9)
    pca = PCA(n_components=1).fit(Xs)
    loadings = pca.components_[0] * np.sqrt(pca.explained_variance_[0])
    loadings = np.abs(loadings)
    # CR = (Σλ)² / [(Σλ)² + Σ(1-λ²)]
    sum_l = loadings.sum()
    sum_l2 = (loadings**2).sum()
    cr = sum_l**2 / (sum_l**2 + (len(loadings) - sum_l2))
    ave = sum_l2 / len(loadings)
    print(f"  {label}: CR = {cr:.4f}, AVE = {ave:.4f}  "
          f"{'(AVE>=0.5 ✓)' if ave>=0.5 else '(AVE<0.5 ✗)'} "
          f"{'(CR>=0.7 ✓)' if cr>=0.7 else '(CR<0.7 ✗)'}")
    return cr, ave

print("\n=== Composite Reliability (CR) and AVE ===")
results = {}
for label, items in [('V', V_n), ('C', C_n), ('F', F_n)]:
    results[label] = cr_ave(items, label)

# ======================================================================
# 3. Convergent validity: correlation between RE and existing indices
# ======================================================================
print("\n=== Convergent Validity (correlation with established indices) ===")
df_test = df.copy()
df_test['icrg_qog'] = df_test['icrg_qog']  # ICRG QoG
df_test['fi_index'] = pd.read_csv(PROC/"qog_g20_re_vars.csv")['fi_index']

# Use original RE_0_100 vs individual benchmark indicators
benchmarks = {
    'WGI Government Effectiveness': 'wbgi_gee',
    'ICRG Quality of Government':   'icrg_qog',
    'WGI Regulatory Quality':       'wbgi_rqe',
    'WGI Rule of Law':              'wbgi_rle',
}
for name, col in benchmarks.items():
    sub = df.dropna(subset=[col,'RE_0_100'])
    r, p = stats.pearsonr(sub[col], sub['RE_0_100'])
    s, sp = stats.spearmanr(sub[col], sub['RE_0_100'])
    print(f"  RE vs {name}: Pearson r={r:.3f} (p={p:.3g}), Spearman ρ={s:.3f}")

# ======================================================================
# 4. Temporal stability (rolling correlation across years)
# ======================================================================
print("\n=== Temporal Stability (year-on-year rank correlation) ===")
years = sorted(df.year.unique())
ranks = df.pivot(index='iso3', columns='year', values='RE_0_100')
ts = []
for y in years[1:]:
    if y-1 in ranks.columns and y in ranks.columns:
        r = stats.spearmanr(ranks[y-1], ranks[y])[0]
        ts.append((y, r))
ts_df = pd.DataFrame(ts, columns=['year','spearman_with_prev'])
print(f"  Mean year-to-year rank correlation: {ts_df['spearman_with_prev'].mean():.4f}")
print(f"  Min: {ts_df['spearman_with_prev'].min():.4f}, Max: {ts_df['spearman_with_prev'].max():.4f}")
ts_df.to_csv(OUT/"temporal_stability.csv", index=False)

# ======================================================================
# 5. Predictive validity: RE 2019 → COVID-19 response (2020-2021)
# ======================================================================
print("\n=== Predictive Validity: RE_2019 → Crisis Response ===")
# Proxy: GDP growth decline in 2020
re_2019 = df[df.year==2019].set_index('iso3')['RE_0_100']
growth_2020 = df[df.year==2020].set_index('iso3')['gdpgrowth'] if 'gdpgrowth' in df.columns else None

# Reload gdpgrowth from raw
gdp = pd.read_csv(PROC/"re_master_raw.csv")[['iso3','year','gdpgrowth']]
g2020 = gdp[gdp.year==2020].set_index('iso3')['gdpgrowth']
g2021 = gdp[gdp.year==2021].set_index('iso3')['gdpgrowth']
combined = pd.DataFrame({'RE_2019': re_2019, 'GDP_2020': g2020, 'GDP_2021': g2021,
                        'Recovery_index': g2020 + g2021}).dropna()
r1, p1 = stats.pearsonr(combined['RE_2019'], combined['GDP_2020'])
r2, p2 = stats.pearsonr(combined['RE_2019'], combined['Recovery_index'])
print(f"  RE_2019 vs GDP_growth_2020:   r={r1:.3f} (p={p1:.3g})")
print(f"  RE_2019 vs 2020+2021 growth:  r={r2:.3f} (p={p2:.3g})")
combined.to_csv(OUT/"predictive_validity_covid.csv")

# ======================================================================
# 6. Discriminant validity: HTMT ratio
# ======================================================================
print("\n=== Discriminant Validity (Inter-dimension correlations) ===")
dim_corr = df[['V_score','C_score','F_score']].corr()
print(dim_corr.round(3))

# ======================================================================
# 7. Save reliability/validity summary
# ======================================================================
summary = []
for label, (cr, ave) in results.items():
    items = {'V':V_n,'C':C_n,'F':F_n}[label]
    a = cronbach_alpha(items)
    summary.append({
        'dimension': label,
        'n_items':   items.shape[1],
        'cronbach_alpha': round(a,4),
        'composite_reliability': round(cr,4),
        'AVE': round(ave,4),
        'alpha_acceptable': a>=0.7,
        'CR_acceptable': cr>=0.7,
        'AVE_acceptable': ave>=0.5,
    })
sum_df = pd.DataFrame(summary)
sum_df.to_csv(OUT/"reliability_validity_summary.csv", index=False)
print("\n=== SUMMARY ===")
print(sum_df.to_string(index=False))
"""
06_visualizations.py
RE指数可视化套件
"""
import pandas as pd
import numpy as np
import matplotlib.pyplot as plt
import matplotlib as mpl
import seaborn as sns
from pathlib import Path

OUT = Path("/home/user/re_index/output")
FIG = Path("/home/user/re_index/output/figures")
FIG.mkdir(exist_ok=True)

mpl.rcParams['font.family'] = 'DejaVu Sans'
mpl.rcParams['axes.unicode_minus'] = False
sns.set_style('whitegrid')

df = pd.read_csv(OUT/"G20_RE_Index_2000_2023.csv")
weights_V = pd.read_csv(OUT/"weights_V.csv")
weights_C = pd.read_csv(OUT/"weights_C.csv")
weights_F = pd.read_csv(OUT/"weights_F.csv")

COUNTRY_SHORT = {
    "Argentina":"ARG","Australia":"AUS","Brazil":"BRA","Canada":"CAN","China":"CHN",
    "Germany":"DEU","France":"FRA","India":"IND","Indonesia":"IDN","Italy":"ITA",
    "Japan":"JPN","Korea (the Republic of)":"KOR","Mexico":"MEX",
    "Russian Federation (the)":"RUS","Saudi Arabia":"SAU","South Africa":"ZAF",
    "Turkey":"TUR","Türkiye":"TUR",
    "United Kingdom of Great Britain and Northern Ireland (the)":"GBR",
    "United States of America (the)":"USA"
}
df['iso3'] = df['country'].map(COUNTRY_SHORT).fillna(df['iso3'])

# ============================================================
# Figure 1: Heatmap of RE_0_100 across countries × years
# ============================================================
fig, ax = plt.subplots(figsize=(14, 7))
pivot = df.pivot(index='iso3', columns='year', values='RE_0_100')
order = pivot.mean(axis=1).sort_values(ascending=False).index
pivot = pivot.loc[order]
sns.heatmap(pivot, cmap='RdYlGn', center=50, ax=ax, cbar_kws={'label':'RE Index (0-100)'},
            annot=False, linewidths=0.3)
ax.set_title('Figure 1: RE Index Heatmap by Country and Year (2000–2023)', fontsize=14, fontweight='bold')
ax.set_xlabel('Year', fontsize=11); ax.set_ylabel('Country (ISO3)', fontsize=11)
plt.tight_layout()
plt.savefig(FIG/"fig1_heatmap.png", dpi=150, bbox_inches='tight')
plt.close()
print("✓ Figure 1: Heatmap saved")

# ============================================================
# Figure 2: Time-series for major countries
# ============================================================
fig, ax = plt.subplots(figsize=(14, 7))
focus = ['USA','CHN','DEU','JPN','KOR','GBR','FRA','IND','BRA','RUS']
colors = plt.cm.tab10(np.linspace(0,1,len(focus)))
for c, col in zip(focus, colors):
    sub = df[df.iso3==c].sort_values('year')
    if len(sub)>0:
        ax.plot(sub.year, sub.RE_0_100, marker='o', label=c, color=col, linewidth=2, markersize=4)
ax.set_title('Figure 2: RE Index Evolution for Major Economies (2000–2023)', fontsize=14, fontweight='bold')
ax.set_xlabel('Year'); ax.set_ylabel('RE Index (0-100)')
ax.legend(loc='center left', bbox_to_anchor=(1, 0.5), fontsize=10)
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(FIG/"fig2_timeseries.png", dpi=150, bbox_inches='tight')
plt.close()
print("✓ Figure 2: Time series saved")

# ============================================================
# Figure 3: V-C-F decomposition for latest year
# ============================================================
fig, ax = plt.subplots(figsize=(14, 7))
latest = df[df.year==2022].sort_values('RE', ascending=True)
x = np.arange(len(latest))
ax.barh(x-0.25, latest.V_score, height=0.25, label='V (Velocity)', color='#2E86AB')
ax.barh(x,      latest.C_score, height=0.25, label='C (Conversion)', color='#A23B72')
ax.barh(x+0.25, latest.F_score, height=0.25, label='F (Friction)', color='#F18F01')
ax.set_yticks(x); ax.set_yticklabels(latest.iso3)
ax.set_xlabel('Normalized Score (0-1)')
ax.set_title('Figure 3: V-C-F Decomposition by Country (2022)', fontsize=14, fontweight='bold')
ax.legend(loc='lower right')
plt.tight_layout()
plt.savefig(FIG/"fig3_VCF_decomposition.png", dpi=150, bbox_inches='tight')
plt.close()
print("✓ Figure 3: VCF decomposition saved")

# ============================================================
# Figure 4: PCA weights bar chart
# ============================================================
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
for ax, (wdf, title, color) in zip(axes, 
    [(weights_V,'V (Velocity)','#2E86AB'),
     (weights_C,'C (Conversion)','#A23B72'),
     (weights_F,'F (Friction)','#F18F01')]):
    wdf = wdf.sort_values('weight_PCA')
    ax.barh(wdf.variable, wdf.weight_PCA, color=color)
    ax.set_title(f'{title}: PCA Weights')
    ax.set_xlabel('Weight (Σ=1)')
    for i, (v, w) in enumerate(zip(wdf.variable, wdf.weight_PCA)):
        ax.text(w+0.005, i, f'{w:.3f}', va='center', fontsize=8)
plt.suptitle('Figure 4: PCA-Derived Weights for Each Dimension', fontsize=14, fontweight='bold')
plt.tight_layout()
plt.savefig(FIG/"fig4_pca_weights.png", dpi=150, bbox_inches='tight')
plt.close()
print("✓ Figure 4: PCA weights saved")

# ============================================================
# Figure 5: RE with 90% bootstrap CI (2022)
# ============================================================
fig, ax = plt.subplots(figsize=(13, 7))
latest = df[df.year==2022].sort_values('RE')
y = np.arange(len(latest))
errs = [latest.RE - latest.RE_lower90, latest.RE_upper90 - latest.RE]
ax.errorbar(latest.RE, y, xerr=errs, fmt='o', color='darkred', 
            ecolor='gray', capsize=4, markersize=8)
for i, (re, iso) in enumerate(zip(latest.RE, latest.iso3)):
    ax.text(re+0.01, i, iso, va='center', fontsize=9)
ax.set_yticks(y); ax.set_yticklabels([])
ax.set_xlabel('RE Index with 90% Bootstrap CI')
ax.set_title('Figure 5: RE Index with Bootstrap Uncertainty (2022) — Innovation #5',
             fontsize=14, fontweight='bold')
ax.grid(True, alpha=0.3)
plt.tight_layout()
plt.savefig(FIG/"fig5_bootstrap_ci.png", dpi=150, bbox_inches='tight')
plt.close()
print("✓ Figure 5: Bootstrap CI saved")

# ============================================================
# Figure 6: Domain sub-indices comparison
# ============================================================
fig, axes = plt.subplots(1, 3, figsize=(16, 6))
latest = df[df.year==2022].sort_values('RE', ascending=False)
for ax, (col, title, color) in zip(axes,
    [('RE_fiscal','Fiscal Domain RE','#2E86AB'),
     ('RE_tech','Technology Domain RE','#A23B72'),
     ('RE_gov','Governance Domain RE','#F18F01')]):
    sorted_d = latest.sort_values(col, ascending=True).head(19)
    ax.barh(sorted_d.iso3, sorted_d[col], color=color)
    ax.set_title(title)
    ax.set_xlabel('Sub-Index')
plt.suptitle('Figure 6: Domain-Specific RE Sub-Indices (2022) — Innovation #4',
             fontsize=14, fontweight='bold')
plt.tight_layout()
plt.savefig(FIG/"fig6_domain_subindices.png", dpi=150, bbox_inches='tight')
plt.close()
print("✓ Figure 6: Domain sub-indices saved")

# ============================================================
# Figure 7: Inter-dimension correlation matrix
# ============================================================
fig, ax = plt.subplots(figsize=(7, 6))
corr = df[['V_score','C_score','F_score','RE','RE_eq','RE_timevar']].corr()
sns.heatmap(corr, annot=True, fmt='.3f', cmap='RdBu_r', center=0, vmin=-1, vmax=1, ax=ax)
ax.set_title('Figure 7: Correlation Matrix — Dimensions & RE Variants',
             fontsize=13, fontweight='bold')
plt.tight_layout()
plt.savefig(FIG/"fig7_correlation_matrix.png", dpi=150, bbox_inches='tight')
plt.close()
print("✓ Figure 7: Correlation matrix saved")

# ============================================================
# Figure 8: Radar chart for selected countries
# ============================================================
fig = plt.figure(figsize=(14, 10))
focus = ['USA','CHN','DEU','JPN','KOR','BRA','IND','RUS']
n = len(focus)
ncols=4; nrows=2
for i, c in enumerate(focus):
    ax = fig.add_subplot(nrows, ncols, i+1, projection='polar')
    sub = df[(df.iso3==c)&(df.year==2022)]
    if len(sub)==0: continue
    s = sub.iloc[0]
    vals = [s['V_score'], s['C_score'], 1-s['F_score']]  # invert F so higher=better
    cats = ['V\n(Velocity)','C\n(Conversion)','1-F\n(Low Friction)']
    angles = np.linspace(0, 2*np.pi, len(cats), endpoint=False)
    vals_p = list(vals)+[vals[0]]
    angles_p = list(angles)+[angles[0]]
    ax.plot(angles_p, vals_p, 'o-', linewidth=2, color='#A23B72')
    ax.fill(angles_p, vals_p, alpha=0.25, color='#A23B72')
    ax.set_xticks(angles); ax.set_xticklabels(cats, fontsize=9)
    ax.set_ylim(0,1); ax.set_yticks([0.2,0.4,0.6,0.8])
    ax.set_yticklabels(['0.2','0.4','0.6','0.8'], fontsize=7)
    ax.set_title(f'{c} (RE={s.RE_0_100:.1f})', fontsize=11, fontweight='bold', pad=15)
plt.suptitle('Figure 8: V-C-F Radar Profiles (2022)', fontsize=14, fontweight='bold', y=1.0)
plt.tight_layout()
plt.savefig(FIG/"fig8_radar.png", dpi=150, bbox_inches='tight')
plt.close()
print("✓ Figure 8: Radar charts saved")

# ============================================================
# Figure 9: Time-varying weights evolution
# ============================================================
tv_V = pd.read_csv(OUT/"time_varying_weights_V.csv")
tv_C = pd.read_csv(OUT/"time_varying_weights_C.csv")
tv_F = pd.read_csv(OUT/"time_varying_weights_F.csv")
fig, axes = plt.subplots(1, 3, figsize=(18, 5))
for ax, (tv, title, cmap) in zip(axes,
    [(tv_V,'V Dimension Time-Varying Weights','Blues'),
     (tv_C,'C Dimension Time-Varying Weights','Purples'),
     (tv_F,'F Dimension Time-Varying Weights','Oranges')]):
    years = tv['year'].values
    cols = [c for c in tv.columns if c!='year']
    cm = plt.get_cmap(cmap)
    for j, c in enumerate(cols):
        ax.plot(years, tv[c], marker='.', label=c, 
                color=cm(0.3+0.6*j/len(cols)), linewidth=2)
    ax.set_title(title)
    ax.set_xlabel('Year'); ax.set_ylabel('Weight')
    ax.legend(fontsize=7, loc='best')
    ax.grid(True, alpha=0.3)
plt.suptitle('Figure 9: Dynamic Time-Varying Weights (Innovation #3)',
             fontsize=14, fontweight='bold')
plt.tight_layout()
plt.savefig(FIG/"fig9_timevarying_weights.png", dpi=150, bbox_inches='tight')
plt.close()
print("✓ Figure 9: Time-varying weights saved")

# ============================================================
# Figure 10: RE vs benchmark indices (scatter)
# ============================================================
fig, axes = plt.subplots(1, 3, figsize=(16, 5))
bench = [('wbgi_gee','WGI Gov Effectiveness'),
         ('icrg_qog','ICRG QoG'),
         ('wbgi_rqe','WGI Regulatory Quality')]
for ax, (col, name) in zip(axes, bench):
    ax.scatter(df[col], df['RE_0_100'], alpha=0.4, s=20)
    from scipy.stats import pearsonr
    r, p = pearsonr(df[col].dropna(), df.loc[df[col].notna(),'RE_0_100'])
    z = np.polyfit(df[col].dropna(), df.loc[df[col].notna(),'RE_0_100'], 1)
    xs = np.linspace(df[col].min(), df[col].max(), 50)
    ax.plot(xs, z[0]*xs+z[1], 'r--', linewidth=2)
    ax.set_xlabel(name); ax.set_ylabel('RE Index (0-100)')
    ax.set_title(f'RE vs {name}\n(r={r:.3f})')
    ax.grid(True, alpha=0.3)
plt.suptitle('Figure 10: Convergent Validity — RE vs Existing Indices',
             fontsize=14, fontweight='bold')
plt.tight_layout()
plt.savefig(FIG/"fig10_convergent_validity.png", dpi=150, bbox_inches='tight')
plt.close()
print("✓ Figure 10: Convergent validity saved")

print(f"\n=== All figures saved to {FIG} ===")
print("Files:")
for f in sorted(FIG.glob("*.png")):
    print(f"  - {f.name}")
"""
07_export_excel.py
导出多工作表Excel数据集
"""
import pandas as pd
from pathlib import Path

OUT = Path("/home/user/re_index/output")
PROC = Path("/home/user/re_index/processed")

re_data = pd.read_csv(OUT/"G20_RE_Index_2000_2023.csv")
raw_data = pd.read_csv(PROC/"re_master_raw.csv")
imputed = pd.read_csv(PROC/"re_data_imputed.csv")
mapping = pd.read_csv(PROC/"re_variable_dimension_map.csv")
weights_V = pd.read_csv(OUT/"weights_V.csv")
weights_C = pd.read_csv(OUT/"weights_C.csv")
weights_F = pd.read_csv(OUT/"weights_F.csv")
rv_summary = pd.read_csv(OUT/"reliability_validity_summary.csv")
ts_stab = pd.read_csv(OUT/"temporal_stability.csv")
coverage = pd.read_csv(PROC/"variable_coverage.csv")
pca_meta = pd.read_csv(OUT/"pca_metadata.csv")

# Master summary by country (mean over years)
country_summary = re_data.groupby(['iso3','country']).agg(
    n_years=('year','count'),
    RE_mean=('RE_0_100','mean'),
    RE_std=('RE_0_100','std'),
    RE_min=('RE_0_100','min'),
    RE_max=('RE_0_100','max'),
    V_mean=('V_score','mean'),
    C_mean=('C_score','mean'),
    F_mean=('F_score','mean')
).round(3).reset_index()

# Year summary
year_summary = re_data.groupby('year').agg(
    RE_mean=('RE_0_100','mean'),
    RE_std=('RE_0_100','std'),
    V_mean=('V_score','mean'),
    C_mean=('C_score','mean'),
    F_mean=('F_score','mean')
).round(3).reset_index()

# Data dictionary
data_dict = pd.DataFrame([
    {'column':'iso3','description':'ISO 3166-1 alpha-3 country code'},
    {'column':'country','description':'Country full name'},
    {'column':'year','description':'Calendar year (2000–2023)'},
    {'column':'wbgi_*','description':'World Bank Worldwide Governance Indicators (WGI) estimates (-2.5 to +2.5)'},
    {'column':'fi_*','description':'Fraser Institute Economic Freedom of the World components (0-10)'},
    {'column':'vdem_*','description':'V-Dem (Varieties of Democracy) indicators (0-1)'},
    {'column':'icrg_qog','description':'ICRG Quality of Government composite (0-1, higher better)'},
    {'column':'wdi_*','description':'World Bank World Development Indicators'},
    {'column':'C_rd','description':'R&D expenditure as % of GDP'},
    {'column':'C_hitech','description':'High-tech exports as % of manufactured exports'},
    {'column':'V_govexp','description':'General government final consumption expenditure (% of GDP)'},
    {'column':'V_score','description':'Composite V (Velocity/Reaction Elasticity) score [0,1]'},
    {'column':'C_score','description':'Composite C (Conversion) score [0,1]'},
    {'column':'F_score','description':'Composite F (Friction) score [0,1]'},
    {'column':'RE','description':'Resource Reallocation Efficiency Index = V·C/(1+F)'},
    {'column':'RE_eq','description':'Alternative RE using equal-weights (robustness)'},
    {'column':'lnRE','description':'Log-transformed RE = ln(V)+ln(C)-ln(1+F)'},
    {'column':'RE_0_100','description':'RE rescaled to 0-100 (main reporting metric)'},
    {'column':'RE_timevar','description':'RE with time-varying PCA weights (Innovation #3)'},
    {'column':'RE_fiscal','description':'Fiscal-domain RE sub-index (Innovation #4)'},
    {'column':'RE_tech','description':'Technology-domain RE sub-index (Innovation #4)'},
    {'column':'RE_gov','description':'Governance-domain RE sub-index (Innovation #4)'},
    {'column':'RE_mean','description':'Bootstrap mean of RE (B=500)'},
    {'column':'RE_lower90','description':'5th percentile (lower 90% CI) of bootstrap RE'},
    {'column':'RE_upper90','description':'95th percentile (upper 90% CI) of bootstrap RE'},
    {'column':'RE_se','description':'Bootstrap standard error of RE'},
])

# Data source inventory
sources = pd.DataFrame([
    {'source':'World Bank WGI','url':'https://www.worldbank.org/en/publication/worldwide-governance-indicators','variables':'wbgi_gee, wbgi_rqe, wbgi_cce, wbgi_rle, wbgi_vae, wbgi_pve','license':'Open (CC BY 4.0)','access_method':'World Bank API source=3'},
    {'source':'Quality of Government (QoG) Standard TS','url':'https://www.qogdata.pol.gu.se/data/qog_std_ts_jan25.csv','variables':'Wrapper for V-Dem, ICRG, Fraser, etc.','license':'Free academic use','access_method':'Direct CSV download (~62MB)'},
    {'source':'V-Dem v15','url':'https://www.v-dem.net/data/the-v-dem-dataset/','variables':'vdem_corr, vdem_libdem, vdem_academ, etc.','license':'Open via QoG','access_method':'Via QoG (since QoG embeds V-Dem)'},
    {'source':'Fraser Economic Freedom of the World','url':'https://efotw.org/economic-freedom/dataset','variables':'fi_index, fi_reg, fi_legprop, fi_sog','license':'Free for research','access_method':'Via QoG'},
    {'source':'ICRG (via QoG aggregate)','url':'https://www.prsgroup.com/','variables':'icrg_qog (composite)','license':'Aggregate index released by QoG','access_method':'Via QoG'},
    {'source':'World Bank WDI','url':'https://api.worldbank.org/v2/','variables':'GB.XPD.RSDV.GD.ZS, TX.VAL.TECH.MF.ZS, GC.XPN.TOTL.GD.ZS, NY.GDP.*','license':'Open (CC BY 4.0)','access_method':'World Bank API'},
])

# Write Excel
with pd.ExcelWriter(OUT/"G20_RE_Index_2000_2023.xlsx", engine='openpyxl') as writer:
    re_data.to_excel(writer, sheet_name='RE_Index_Master', index=False)
    country_summary.to_excel(writer, sheet_name='Country_Summary', index=False)
    year_summary.to_excel(writer, sheet_name='Year_Summary', index=False)
    raw_data.to_excel(writer, sheet_name='Raw_Data', index=False)
    imputed.to_excel(writer, sheet_name='Imputed_Data', index=False)
    weights_V.to_excel(writer, sheet_name='Weights_V', index=False)
    weights_C.to_excel(writer, sheet_name='Weights_C', index=False)
    weights_F.to_excel(writer, sheet_name='Weights_F', index=False)
    rv_summary.to_excel(writer, sheet_name='Reliability_Validity', index=False)
    ts_stab.to_excel(writer, sheet_name='Temporal_Stability', index=False)
    coverage.to_excel(writer, sheet_name='Coverage_Report', index=False)
    pca_meta.to_excel(writer, sheet_name='PCA_Metadata', index=False)
    mapping.to_excel(writer, sheet_name='Variable_Dimension_Map', index=False)
    data_dict.to_excel(writer, sheet_name='Data_Dictionary', index=False)
    sources.to_excel(writer, sheet_name='Data_Sources', index=False)

print(f"✓ Excel saved: {OUT/'G20_RE_Index_2000_2023.xlsx'}")
import os
print(f"  File size: {os.path.getsize(OUT/'G20_RE_Index_2000_2023.xlsx')/1024:.1f} KB")
print(f"  Sheets: 15 worksheets")