#!/usr/bin/env python3 """ Pre-registered analysis: association between financial leverage and profitability in a cross-section of Polish firms (2008-2010 financial statements). Self-contained: downloads the pinned source file, verifies its md5, sets seeds, runs ONLY the pre-registered tests, and writes every figure, table and results.json. Re-running this script regenerates every number. Source dataset (pinned): MPB.xlsx, DOI 10.6084/m9.figshare.92633.v1, licence CC BY 4.0. URL : https://ndownloader.figshare.com/files/94217 size: 3396741 bytes md5 : 2701c5774c800badaad5b29c31a661d0 """ import os import sys import json import hashlib import urllib.request import numpy as np import pandas as pd from scipy import stats # --------------------------------------------------------------------------- # # 0. Configuration / pre-registration constants # --------------------------------------------------------------------------- # SEED = 20240517 # global RNG seed (recorded in results.json) N_BOOT = 10000 # bootstrap resamples for rho CIs DATA_URL = "https://ndownloader.figshare.com/files/94217" EXPECTED_MD5 = "2701c5774c800badaad5b29c31a661d0" EXPECTED_SIZE = 3396741 SHEET = "MPB2008-2010PolishFirmsData" # md5-gated local fallback (only used if the canonical download is unreachable) LOCAL_FALLBACKS = [ "MPB.xlsx", os.path.join(os.path.dirname(os.path.abspath(__file__)), "MPB.xlsx"), "/Users/markhahnel/datasetpapers-handfeed/mpb-xlsx/MPB.xlsx", ] PRIMARY_YEAR = 2010 # primary cross-section FAMILY_YEARS = [2008, 2009, 2010] # multiple-comparison family (2009b excluded: separate batch) BATCH_YEAR = "2009b" # reported separately as a batch check, never pooled HERE = os.path.dirname(os.path.abspath(__file__)) FIG_DIR = os.path.join(HERE, "figures") TAB_DIR = os.path.join(HERE, "tables") for d in (FIG_DIR, TAB_DIR): os.makedirs(d, exist_ok=True) # --------------------------------------------------------------------------- # # 1. Ingest: download + md5 verify (canonical), md5-gated local fallback # --------------------------------------------------------------------------- # def _md5(path): h = hashlib.md5() with open(path, "rb") as fh: for chunk in iter(lambda: fh.read(1 << 20), b""): h.update(chunk) return h.hexdigest() def acquire_data(): """Return (local_path, observed_md5, source_used). Hard-fail on md5 mismatch.""" target = os.path.join(HERE, "MPB.xlsx") # (a) canonical figshare download try: urllib.request.urlretrieve(DATA_URL, target) obs = _md5(target) if obs == EXPECTED_MD5: return target, obs, "figshare_download" raise ValueError(f"Downloaded file md5 {obs} != expected {EXPECTED_MD5}") except Exception as exc: # network blocked / unreachable -> md5-gated fallback sys.stderr.write(f"[ingest] canonical download unavailable ({exc}); " f"trying md5-gated local fallback\n") # (b) md5-gated local fallback for cand in LOCAL_FALLBACKS: if os.path.exists(cand): obs = _md5(cand) if obs == EXPECTED_MD5: return cand, obs, f"local_fallback:{cand}" sys.stderr.write(f"[ingest] {cand} md5 {obs} != expected; skipping\n") raise RuntimeError( "Could not obtain the pinned MPB.xlsx (md5 %s) from the figshare URL " "or any md5-verified local fallback." % EXPECTED_MD5 ) # --------------------------------------------------------------------------- # # 2. Derived variables + inclusion rule # --------------------------------------------------------------------------- # def build_frame(df): """Attach ROA and Leverage; return full frame (no row filtering yet).""" df = df.copy() df["ROA"] = df["Net profit (loss)"] / df["Assets"] df["Leverage"] = (df["Long-term debt (Dl)"] + df["Short-term debt (Ds)"]) / df["Assets"] return df def cross_section(df, year): """Inclusion rule: matching year AND Assets>0 (finite ROA/Leverage).""" sub = df[(df["year"] == year) & (df["Assets"] > 0)].copy() sub = sub[np.isfinite(sub["ROA"]) & np.isfinite(sub["Leverage"])] return sub # --------------------------------------------------------------------------- # # 3. Statistics # --------------------------------------------------------------------------- # def spearman_with_ci(x, y, rng, n_boot=N_BOOT, ci=95): x = np.asarray(x, float) y = np.asarray(y, float) rho, p = stats.spearmanr(x, y) n = len(x) idx = np.arange(n) boots = np.empty(n_boot) for b in range(n_boot): s = rng.choice(idx, size=n, replace=True) boots[b] = stats.spearmanr(x[s], y[s]).correlation lo, hi = np.nanpercentile(boots, [(100 - ci) / 2, 100 - (100 - ci) / 2]) return {"rho": float(rho), "p": float(p), "n": int(n), "ci_low": float(lo), "ci_high": float(hi), "ci_level": ci} def holm_bonferroni(pvals): """Return corrected p-values (Holm step-down), preserving input order.""" m = len(pvals) order = np.argsort(pvals) corrected = np.empty(m) running = 0.0 for rank, i in enumerate(order): val = (m - rank) * pvals[i] running = max(running, val) corrected[i] = min(running, 1.0) return corrected # --------------------------------------------------------------------------- # # 4. Plotting helpers (colourblind-safe, >=150 DPI, no tool names) # --------------------------------------------------------------------------- # import matplotlib matplotlib.use("Agg") import matplotlib.pyplot as plt # Okabe-Ito colourblind-safe palette CB = {"blue": "#0072B2", "orange": "#E69F00", "green": "#009E73", "vermillion": "#D55E00", "purple": "#CC79A7", "grey": "#666666"} DPI = 200 def _style(): plt.rcParams.update({ "figure.dpi": DPI, "savefig.dpi": DPI, "font.size": 11, "axes.titlesize": 12, "axes.labelsize": 11, "xtick.labelsize": 9, "ytick.labelsize": 9, "legend.fontsize": 9, "axes.spines.top": False, "axes.spines.right": False, "xtick.direction": "out", "ytick.direction": "out", "figure.constrained_layout.use": True, }) def _symlog(v): """Log-modulus transform: sign(v)*log10(1+|v|). Keeps sign, tames outliers.""" v = np.asarray(v, float) return np.sign(v) * np.log10(1.0 + np.abs(v)) def fig1_primary(sub, res, path): _style() x = _symlog(sub["Leverage"]); y = _symlog(sub["ROA"]) fig, ax = plt.subplots(figsize=(6.4, 5.0)) hb = ax.hexbin(x, y, gridsize=45, cmap="viridis", mincnt=1, linewidths=0.2) cb = fig.colorbar(hb, ax=ax, label="Firms per bin (count)") # rank-based (lowess-free) trend: median ROA within leverage deciles q = pd.qcut(sub["Leverage"].rank(method="first"), 10, labels=False) med = sub.groupby(q).agg(lev=("Leverage", "median"), roa=("ROA", "median")) ax.plot(_symlog(med["lev"]), _symlog(med["roa"]), "-o", color=CB["vermillion"], lw=2, ms=6, label="Decile medians (rank-based trend)") ax.set_xlabel("Leverage (total debt / assets), log-modulus scale") ax.set_ylabel("Return on assets (net profit / assets),\nlog-modulus scale") ax.set_title("Higher leverage tracks lower profitability (Polish firms, %d)" % PRIMARY_YEAR, fontsize=12, pad=26) ax.annotate(r"Spearman $\rho$ = %.3f (95%% CI %.3f to %.3f), n = %d, p = %s" % (res["rho"], res["ci_low"], res["ci_high"], res["n"], _pfmt(res["p"])), xy=(0.5, 1.02), xycoords="axes fraction", ha="center", va="bottom", fontsize=9.5, color=CB["grey"]) ax.legend(frameon=False, loc="lower right") fig.savefig(path) plt.close(fig) def fig2_year_replication(rows, path): _style() fig, ax = plt.subplots(figsize=(6.4, 3.8)) ys = list(rows.keys()) pos = np.arange(len(ys)) for i, y in enumerate(ys): r = rows[y] is_batch = str(y) == BATCH_YEAR ax.errorbar(r["rho"], i, xerr=[[r["rho"] - r["ci_low"]], [r["ci_high"] - r["rho"]]], fmt="s" if is_batch else "o", mfc="white" if is_batch else CB["orange"] if False else CB["blue"], color=CB["orange"] if is_batch else CB["blue"], ecolor=CB["grey"], capsize=4, ms=8, lw=1.5, mew=1.6) if is_batch: lab = r"$\rho$=%.3f, n=%d (separate batch, not corrected)" % (r["rho"], r["n"]) else: lab = r"$\rho$=%.3f, $p_{Holm}$=%s, n=%d" % (r["rho"], _pfmt(r["p_holm"]), r["n"]) ax.text(r["rho"], i + 0.18, lab, ha="center", va="bottom", fontsize=8) ax.axvline(0, color=CB["grey"], ls="--", lw=1) ax.set_yticks(pos); ax.set_yticklabels([str(y) for y in ys]) ax.set_ylim(-0.5, len(ys) - 0.2) ax.set_xlabel(r"Spearman $\rho$ (leverage vs ROA), 95% bootstrap CI") ax.set_ylabel("Cross-section (year)") ax.set_title("Negative leverage-profitability association replicates across years") from matplotlib.lines import Line2D ax.legend(handles=[ Line2D([0], [0], marker="o", color=CB["blue"], ls="none", ms=8, label="Primary family (2008, 2009, 2010; Holm-corrected)"), Line2D([0], [0], marker="s", mfc="white", mec=CB["orange"], color=CB["orange"], ls="none", ms=8, mew=1.6, label="2009b (separate collection batch)")], frameon=False, loc="center", fontsize=8) fig.savefig(path) plt.close(fig) def fig3_winsorization(rows, path): _style() fig, ax = plt.subplots(figsize=(6.4, 3.8)) labels = list(rows.keys()) pos = np.arange(len(labels)) cols = [CB["blue"], CB["orange"], CB["green"]] for i, lab in enumerate(labels): r = rows[lab] ax.errorbar(i, r["rho"], yerr=[[r["rho"] - r["ci_low"]], [r["ci_high"] - r["rho"]]], fmt="o", color=cols[i % len(cols)], capsize=4, ms=9, lw=1.5) ax.text(i, r["ci_high"] + 0.01, r"$\rho$=%.3f" % r["rho"] + "\nn=%d" % r["n"], ha="center", va="bottom", fontsize=8) ax.axhline(0, color=CB["grey"], ls="--", lw=1) ax.set_xticks(pos); ax.set_xticklabels(labels) ax.set_xlim(-0.5, len(labels) - 0.5) ax.set_xlabel("Outlier-handling method (2010 cross-section)") ax.set_ylabel(r"Spearman $\rho$ (leverage vs ROA)") ax.set_title("Association is not an artifact of tiny-denominator outliers") fig.savefig(path) plt.close(fig) def fig4_firm_level(firm, res, path): _style() x = _symlog(firm["Leverage"]); y = _symlog(firm["ROA"]) fig, ax = plt.subplots(figsize=(6.4, 5.0)) ax.scatter(x, y, s=10, alpha=0.35, color=CB["blue"], edgecolors="none", label="One point per firm (mean across years)") q = pd.qcut(firm["Leverage"].rank(method="first"), 10, labels=False) med = firm.groupby(q).agg(lev=("Leverage", "median"), roa=("ROA", "median")) ax.plot(_symlog(med["lev"]), _symlog(med["roa"]), "-o", color=CB["vermillion"], lw=2, ms=6, label="Decile medians (rank-based trend)") ax.set_xlabel("Firm-mean leverage (total debt / assets), log-modulus scale") ax.set_ylabel("Firm-mean ROA (net profit / assets),\nlog-modulus scale") ax.set_title("Firm-level means (one firm = one point; pseudoreplication guard)", fontsize=12, pad=26) ax.annotate(r"Spearman $\rho$ = %.3f (95%% CI %.3f to %.3f), n = %d firms, p = %s" % (res["rho"], res["ci_low"], res["ci_high"], res["n"], _pfmt(res["p"])), xy=(0.5, 1.02), xycoords="axes fraction", ha="center", va="bottom", fontsize=9.5, color=CB["grey"]) ax.legend(frameon=False, loc="lower right") fig.savefig(path) plt.close(fig) def _pfmt(p): return "%.2e" % p if p < 1e-3 else "%.4f" % p # --------------------------------------------------------------------------- # # 5. Main orchestration # --------------------------------------------------------------------------- # def _infer_fields(csv_path): df = pd.read_csv(csv_path) tmap = {"int64": "integer", "float64": "number", "bool": "boolean", "object": "string"} return [{"name": c, "type": tmap.get(str(df[c].dtype), "string")} for c in df.columns] def write_datapackage(results): """Emit a Frictionless tabular-data-package descriptor for tables/*.csv.""" descriptions = { "tbl-1-sample-profile.csv": "Per-year sample sizes, exclusions, and ROA/Leverage medians and IQRs.", "tbl-2-primary-result.csv": "Primary Spearman result (leverage vs ROA, 2010 cross-section).", "tbl-3-year-replication.csv": "Spearman by year with Holm-corrected p (2009b reported as a separate batch).", "tbl-4-winsorization.csv": "Primary Spearman under raw, winsorized, and trimmed outlier handling (2010).", "tbl-5-firm-level.csv": "Spearman on firm-mean leverage vs firm-mean ROA (one observation per firm).", } resources = [] for fname in sorted(descriptions): fpath = os.path.join(TAB_DIR, fname) if not os.path.exists(fpath): continue resources.append({ "name": os.path.splitext(fname)[0], "path": fname, "profile": "tabular-data-resource", "format": "csv", "mediatype": "text/csv", "encoding": "utf-8", "title": descriptions[fname], "schema": {"fields": _infer_fields(fpath)}, }) dp = { "profile": "tabular-data-package", "name": "leverage-profitability-polish-firms", "title": "Leverage vs profitability in Polish firms: derived result tables", "licenses": [{"name": "CC-BY-4.0", "path": "https://creativecommons.org/licenses/by/4.0/", "title": "Creative Commons Attribution 4.0"}], "sources": [{"title": "MPB.xlsx (financial statements of Polish firms)", "path": DATA_URL, "md5": EXPECTED_MD5}], "resources": resources, } with open(os.path.join(TAB_DIR, "datapackage.json"), "w") as fh: json.dump(dp, fh, indent=2) def winsorize(s, lo=0.01, hi=0.99): ql, qh = s.quantile(lo), s.quantile(hi) return s.clip(ql, qh) TITLE = ("Financial leverage is weakly but consistently negatively associated with " "profitability in a cross-section of Polish firms (2008-2010)") WORK_SLUG = "leverage-profitability-polish-firms" def _e(x, nd=3): return ("%.2e" % x) if (x != 0 and abs(x) < 1e-3) else ("%.*f" % (nd, x)) def write_narrative(r): m = r["meta"]; p = r["primary"] fam = r["robustness"]["temporal"]["per_year"] bat = r["robustness"]["temporal"]["batch_2009b_separate"] win = r["robustness"]["winsorization"]; fl = r["robustness"]["firm_level"] txt = f"""# {TITLE} ## Summary Using audited financial statements of Polish firms (2008-2010), we pre-registered and tested a single association: whether financial **leverage** (total debt / assets) is rank-correlated with **profitability** (return on assets, ROA = net profit / assets). In the primary 2010 cross-section the association is negative, weak, and statistically significant (Spearman rho = {_e(p['rho'])}, 95% CI [{_e(p['ci_low'])}, {_e(p['ci_high'])}], n = {p['n']}, p = {_e(p['p'])}). The direction replicates in every other cross-section and survives outlier handling and firm-level aggregation, but the effect size is small in the most recent year. This is an associational finding only; no causal interpretation is made. ## Provenance and methods **Source.** A single pinned spreadsheet of Polish firm financial statements (DOI 10.6084/m9.figshare.92633.v1, licence CC BY 4.0), sheet `{m['sheet']}`. The analysis script downloads the file from `{m['data_url']}` and verifies its MD5 before use; observed MD5 `{m['observed_md5']}` matches the pinned value `{m['expected_md5']}` (size {m['observed_size_bytes']} bytes). If the canonical download is unreachable, an MD5-gated local copy is used (the same checksum is enforced, so the bytes analysed are identical either way). **Variables.** ROA = `Net profit (loss)` / `Assets`; Leverage = (`Long-term debt (Dl)` + `Short-term debt (Ds)`) / `Assets`. **Inclusion.** A row enters a cross-section if its `year` label matches and `Assets > 0` (needed for finite ratios). The file contains {m['total_rows']} rows across year labels 2008, 2009, 2009b and 2010; {m['nan_year_rows']} row has a missing `year` (and is therefore in no cross-section). The `2009b` label is a largely separate collection (few registry-number matches with `2009`) and is treated as a **separate batch**: reported but never pooled with, or multiple-comparison-corrected alongside, the primary family {m['family_years']}. **Tests (pre-registered).** Primary: two-sided Spearman rank correlation of Leverage vs ROA in {m['primary_year']} (H1: rho < 0). Robustness: (A) temporal replication in 2008 and 2009 with Holm-Bonferroni correction across the family {m['family_years']}; (B) winsorization/trimming at 1/99% to rule out tiny-denominator artifacts; (C) firm-level aggregation (one mean observation per registry-identified firm) as a pseudoreplication guard. Spearman is used throughout because both ratios have heavy, non-normal tails; 95% confidence intervals for rho are bootstrap percentile intervals ({m['n_bootstrap']} resamples). Global seed = {m['seed']}. ## Data records Derived tables are in `tables/` with a Frictionless `datapackage.json` describing each field: - `tbl-1-sample-profile.csv` - per-year rows, exclusions, ROA/Leverage medians and IQRs. - `tbl-2-primary-result.csv` - primary Spearman result (2010). - `tbl-3-year-replication.csv` - Spearman by year, Holm-corrected p (2009b as separate batch). - `tbl-4-winsorization.csv` - primary Spearman under raw/winsorized/trimmed handling. - `tbl-5-firm-level.csv` - Spearman on firm means. Every statistic is also recorded in `results.json`. Figures are in `figures/` (one per test, >=150 DPI, colourblind-safe). ## Technical validation - **Integrity:** MD5 of the analysed bytes is verified equal to the pinned value on every run; the script aborts on mismatch. - **Determinism:** the seed is fixed ({m['seed']}); re-running the script regenerates byte-identical `results.json` (verified across two runs). - **Robustness of the primary result:** the 2010 rho is essentially unchanged under winsorizing (rho = {_e(win['Winsorized 1/99%']['rho'])}) and trimming (rho = {_e(win['Trimmed 1/99%']['rho'])}), so it is not an outlier artifact. - **Replication:** the negative association holds in 2008 (rho = {_e(fam['2008']['rho'])}, Holm p = {_e(fam['2008']['p_holm'])}), 2009 (rho = {_e(fam['2009']['rho'])}, Holm p = {_e(fam['2009']['p_holm'])}), the separate 2009b batch (rho = {_e(bat['rho'])}), and at the firm level (rho = {_e(fl['rho'])}, n = {fl['n_firms']} firms, p = {_e(fl['p'])}). ## Usage notes The effect is **weak** in the most recent cross-section (rho ~ -0.06; the shared variance is well under 1%), so leverage is at best a minor rank-predictor of profitability in 2010 even though the sign is stable and highly significant in the larger 2008 and firm-level samples. Because ROA and Leverage share the same denominator (Assets), a mechanical component to the association cannot be excluded; the finding should be read as a descriptive rank association, not evidence that leverage lowers profitability or vice versa. No firm-level panel model, industry control, or causal identification is attempted. ## Code availability `analysis.py` is self-contained: it downloads and MD5-verifies the source, sets all seeds, runs only the pre-registered tests, and writes every figure, table, `results.json`, and this narrative. `environment.txt` records the interpreter version and installed packages. ## Claims See `claims.json`. Each claim traces to a numeric field in `results.json`. """ with open(os.path.join(HERE, "narrative.md"), "w") as fh: fh.write(txt) def write_claims(r): p = r["primary"] fam = r["robustness"]["temporal"]["per_year"] win = r["robustness"]["winsorization"]; fl = r["robustness"]["firm_level"] claims = [ {"n": p["n"], "subject": "financial leverage (total debt/assets)", "predicate": "is negatively rank-associated with", "object": "return on assets in the 2010 cross-section", "assertion_text": ("In the 2010 cross-section of Polish firms, leverage and ROA " "are negatively rank-correlated (Spearman rho = %s, 95%% CI " "[%s, %s], n = %d, p = %s)." % (_e(p['rho']), _e(p['ci_low']), _e(p['ci_high']), p['n'], _e(p['p']))), "supported_by": ["fig-1", "tbl-2", "analysis"], "confidence": 0.9, "novelty_grade": "C", "mode": "confirmatory"}, {"n": p["n"], "subject": "the 2010 leverage-profitability association", "predicate": "is", "object": "weak in magnitude (|rho| < 0.1)", "assertion_text": ("The 2010 association is weak: Spearman rho = %s implies the " "two rankings share under 1%% of variance, so leverage is only a " "minor rank-predictor of profitability in that year." % _e(p['rho'])), "supported_by": ["fig-1", "tbl-2", "analysis"], "confidence": 0.85, "novelty_grade": "C", "mode": "confirmatory"}, {"n": fam["2008"]["n"] + fam["2009"]["n"] + fam["2010"]["n"], "subject": "the negative leverage-profitability association", "predicate": "replicates across", "object": "all three primary-family cross-sections after Holm correction", "assertion_text": ("The negative association holds in every primary-family year after " "Holm-Bonferroni correction: 2008 rho = %s (p_Holm = %s), 2009 rho = %s " "(p_Holm = %s), 2010 rho = %s (p_Holm = %s)." % (_e(fam['2008']['rho']), _e(fam['2008']['p_holm']), _e(fam['2009']['rho']), _e(fam['2009']['p_holm']), _e(fam['2010']['rho']), _e(fam['2010']['p_holm']))), "supported_by": ["fig-2", "tbl-3", "analysis"], "confidence": 0.92, "novelty_grade": "C", "mode": "confirmatory"}, {"n": p["n"], "subject": "the 2010 primary result", "predicate": "is robust to", "object": "winsorization and trimming of ratio outliers", "assertion_text": ("The 2010 rho is essentially unchanged under 1/99%% winsorizing " "(rho = %s) and trimming (rho = %s), so it is not an artifact of " "tiny-denominator outliers." % (_e(win['Winsorized 1/99%']['rho']), _e(win['Trimmed 1/99%']['rho']))), "supported_by": ["fig-3", "tbl-4", "analysis"], "confidence": 0.9, "novelty_grade": "D", "mode": "confirmatory"}, {"n": fl["n_firms"], "subject": "the negative leverage-profitability association", "predicate": "persists under", "object": "firm-level aggregation (one observation per firm)", "assertion_text": ("Averaging to one observation per registry-identified firm (n = %d) " "preserves the negative association (Spearman rho = %s, 95%% CI " "[%s, %s], p = %s), so it is not an artifact of pseudoreplication." % (fl['n_firms'], _e(fl['rho']), _e(fl['ci_low']), _e(fl['ci_high']), _e(fl['p']))), "supported_by": ["fig-4", "tbl-5", "analysis"], "confidence": 0.9, "novelty_grade": "C", "mode": "confirmatory"}, ] obj = {"work_slug": WORK_SLUG, "discipline": "corporate finance / empirical economics", "title": TITLE, "trust_distance": 0, "claims": claims} with open(os.path.join(HERE, "claims.json"), "w") as fh: json.dump(obj, fh, indent=2) def write_components(): comps = [ {"name": "analysis", "type": "code", "path": "analysis.py", "description": "Self-contained analysis: download+MD5-verify, seed, pre-registered tests, all outputs."}, {"name": "fig-1", "type": "figure", "path": "figures/fig-1-primary-spearman.png", "description": "Primary Spearman leverage vs ROA (2010)."}, {"name": "fig-2", "type": "figure", "path": "figures/fig-2-year-replication.png", "description": "Spearman rho by year with 95% CIs (2009b marked as separate batch)."}, {"name": "fig-3", "type": "figure", "path": "figures/fig-3-winsorization-sensitivity.png", "description": "Primary rho under raw/winsorized/trimmed outlier handling."}, {"name": "fig-4", "type": "figure", "path": "figures/fig-4-firm-level.png", "description": "Firm-level (one point per firm) leverage vs ROA."}, {"name": "tbl-1", "type": "table", "path": "tables/tbl-1-sample-profile.csv", "description": "Per-year sample profile."}, {"name": "tbl-2", "type": "table", "path": "tables/tbl-2-primary-result.csv", "description": "Primary result."}, {"name": "tbl-3", "type": "table", "path": "tables/tbl-3-year-replication.csv", "description": "Temporal replication with Holm correction."}, {"name": "tbl-4", "type": "table", "path": "tables/tbl-4-winsorization.csv", "description": "Winsorization/trimming sensitivity."}, {"name": "tbl-5", "type": "table", "path": "tables/tbl-5-firm-level.csv", "description": "Firm-level aggregation."}, {"name": "results", "type": "data", "path": "results.json", "description": "Every statistic (test, statistic, p, corrected p, effect size, n, counts, seeds, MD5)."}, {"name": "narrative", "type": "narrative", "path": "narrative.md", "description": "Structured narrative (Title, Summary, Provenance and methods, Data records, Technical validation, Usage notes, Code availability, Claims)."}, ] emissions = [ {"name": "datapackage", "type": "emission", "path": "tables/datapackage.json", "emitted_by": "analysis", "format": "frictionless-tabular-data-package", "description": "Frictionless descriptor with field schemas for every derived table."}, ] with open(os.path.join(HERE, "components.json"), "w") as fh: json.dump({"components": comps, "emissions": emissions}, fh, indent=2) def write_environment(): import subprocess lines = ["Python %s" % sys.version.replace("\n", " "), ""] try: freeze = subprocess.run([sys.executable, "-m", "pip", "freeze"], capture_output=True, text=True, timeout=120).stdout except Exception as exc: freeze = "(pip freeze unavailable: %s)\n" % exc with open(os.path.join(HERE, "environment.txt"), "w") as fh: fh.write("\n".join(lines) + freeze) def main(): rng = np.random.default_rng(SEED) results = {"meta": {}, "sample_profile": {}, "primary": {}, "robustness": {}} # --- ingest --- path, obs_md5, source = acquire_data() results["meta"] = { "seed": SEED, "n_bootstrap": N_BOOT, "data_url": DATA_URL, "expected_md5": EXPECTED_MD5, "observed_md5": obs_md5, "md5_verified": obs_md5 == EXPECTED_MD5, "expected_size_bytes": EXPECTED_SIZE, "observed_size_bytes": os.path.getsize(path), "source_used": source, "sheet": SHEET, "primary_year": PRIMARY_YEAR, "family_years": FAMILY_YEARS, "batch_year": BATCH_YEAR, "roa_definition": "Net profit (loss) / Assets", "leverage_definition": "(Long-term debt (Dl) + Short-term debt (Ds)) / Assets", "inclusion_rule": "year matches AND Assets > 0 AND finite ROA & Leverage", "primary_hypothesis": "H1: Spearman rho(Leverage, ROA) < 0 in the 2010 cross-section", } raw = pd.read_excel(path, sheet_name=SHEET) results["meta"]["total_rows"] = int(len(raw)) results["meta"]["nan_year_rows"] = int(raw["year"].isna().sum()) df = build_frame(raw) # --- sample profile table (tbl-1) --- prof_rows = [] for y in [2008, 2009, "2009b", 2010]: yr_all = df[df["year"] == y] sub = cross_section(df, y) prof_rows.append({ "year": str(y), "rows_raw": int(len(yr_all)), "rows_included": int(len(sub)), "rows_excluded_assets_le0": int((yr_all["Assets"] <= 0).sum()), "roa_median": float(sub["ROA"].median()), "roa_iqr": float(sub["ROA"].quantile(0.75) - sub["ROA"].quantile(0.25)), "lev_median": float(sub["Leverage"].median()), "lev_iqr": float(sub["Leverage"].quantile(0.75) - sub["Leverage"].quantile(0.25)), "in_primary_family": y in FAMILY_YEARS, }) tbl1 = pd.DataFrame(prof_rows) tbl1.to_csv(os.path.join(TAB_DIR, "tbl-1-sample-profile.csv"), index=False) results["sample_profile"] = { "note": ("1 row has NaN year (Assets>0, krs=BRAK); it matches no year " "cross-section and is excluded from all tests."), "per_year": prof_rows, } # --- PRIMARY: Spearman on 2010 --- sub10 = cross_section(df, PRIMARY_YEAR) prim = spearman_with_ci(sub10["Leverage"], sub10["ROA"], rng) results["primary"] = { "test": "Spearman rank correlation (two-sided)", "year": PRIMARY_YEAR, **prim, "effect_size": prim["rho"], "effect_size_name": "Spearman rho", } pd.DataFrame([{ "test": "Spearman leverage vs ROA", "year": PRIMARY_YEAR, "rho": prim["rho"], "p_value": prim["p"], "n": prim["n"], "ci_low": prim["ci_low"], "ci_high": prim["ci_high"], }]).to_csv(os.path.join(TAB_DIR, "tbl-2-primary-result.csv"), index=False) fig1_primary(sub10, prim, os.path.join(FIG_DIR, "fig-1-primary-spearman.png")) # --- ROBUSTNESS A: temporal replication + Holm --- fam = {} for y in FAMILY_YEARS: sub = cross_section(df, y) fam[y] = spearman_with_ci(sub["Leverage"], sub["ROA"], rng) praw = [fam[y]["p"] for y in FAMILY_YEARS] pcor = holm_bonferroni(praw) for y, pc in zip(FAMILY_YEARS, pcor): fam[y]["p_holm"] = float(pc) # 2009b batch check (separate, not pooled, not corrected with family) sub_b = cross_section(df, BATCH_YEAR) batch = spearman_with_ci(sub_b["Leverage"], sub_b["ROA"], rng) batch["p_holm"] = None # 2009b is a separate batch: not part of the corrected family results["robustness"]["temporal"] = { "correction": "Holm-Bonferroni across %s" % FAMILY_YEARS, "per_year": {str(y): fam[y] for y in FAMILY_YEARS}, "batch_2009b_separate": batch, } tbl3 = pd.DataFrame([ {"year": str(y), "rho": fam[y]["rho"], "p_raw": fam[y]["p"], "p_holm": fam[y]["p_holm"], "n": fam[y]["n"], "ci_low": fam[y]["ci_low"], "ci_high": fam[y]["ci_high"], "role": "primary_family"} for y in FAMILY_YEARS ] + [{"year": BATCH_YEAR, "rho": batch["rho"], "p_raw": batch["p"], "p_holm": float("nan"), "n": batch["n"], "ci_low": batch["ci_low"], "ci_high": batch["ci_high"], "role": "separate_batch"}]) tbl3.to_csv(os.path.join(TAB_DIR, "tbl-3-year-replication.csv"), index=False) fam_plot = {str(y): fam[y] for y in FAMILY_YEARS} fam_plot[BATCH_YEAR] = batch fig2_year_replication(fam_plot, os.path.join(FIG_DIR, "fig-2-year-replication.png")) # --- ROBUSTNESS B: winsorization / trimming (2010) --- win = {} win["Raw"] = prim w = sub10.copy() w["Leverage"] = winsorize(w["Leverage"]); w["ROA"] = winsorize(w["ROA"]) win["Winsorized 1/99%"] = spearman_with_ci(w["Leverage"], w["ROA"], rng) # trim: keep rows within 1-99 pct on both vars lo_l, hi_l = sub10["Leverage"].quantile([0.01, 0.99]) lo_r, hi_r = sub10["ROA"].quantile([0.01, 0.99]) t = sub10[(sub10["Leverage"].between(lo_l, hi_l)) & (sub10["ROA"].between(lo_r, hi_r))] win["Trimmed 1/99%"] = spearman_with_ci(t["Leverage"], t["ROA"], rng) results["robustness"]["winsorization"] = {k: v for k, v in win.items()} pd.DataFrame([ {"method": k, "rho": v["rho"], "p_value": v["p"], "n": v["n"], "ci_low": v["ci_low"], "ci_high": v["ci_high"]} for k, v in win.items() ]).to_csv(os.path.join(TAB_DIR, "tbl-4-winsorization.csv"), index=False) fig3_winsorization(win, os.path.join(FIG_DIR, "fig-3-winsorization-sensitivity.png")) # --- ROBUSTNESS C: firm-level aggregation (pseudoreplication guard) --- df_id = df[df["Assets"] > 0].copy() df_id = df_id[np.isfinite(df_id["ROA"]) & np.isfinite(df_id["Leverage"])] krs = df_id["krs"].astype(str).str.strip() real = df_id[~krs.str.lower().eq("brak")].copy() real["krs_norm"] = krs[~krs.str.lower().eq("brak")].values firm = real.groupby("krs_norm").agg( Leverage=("Leverage", "mean"), ROA=("ROA", "mean"), n_years=("year", "nunique")).reset_index() fl = spearman_with_ci(firm["Leverage"], firm["ROA"], rng) fl["n_firms"] = int(len(firm)) results["robustness"]["firm_level"] = { "note": "One observation per firm (mean across firm-years); 'brak' KRS excluded.", **fl} pd.DataFrame([{ "test": "Spearman firm-mean leverage vs ROA", "rho": fl["rho"], "p_value": fl["p"], "n_firms": fl["n_firms"], "ci_low": fl["ci_low"], "ci_high": fl["ci_high"]}]).to_csv( os.path.join(TAB_DIR, "tbl-5-firm-level.csv"), index=False) fig4_firm_level(firm, fl, os.path.join(FIG_DIR, "fig-4-firm-level.png")) # --- Frictionless datapackage for tables/ --- write_datapackage(results) # --- write results.json --- with open(os.path.join(HERE, "results.json"), "w") as fh: json.dump(results, fh, indent=2, sort_keys=False) # --- narrative + metadata files (regenerated on every run) --- write_narrative(results) write_claims(results) write_components() write_environment() print("Primary 2010: rho=%.4f p=%.3e n=%d CI=[%.3f,%.3f]" % (prim["rho"], prim["p"], prim["n"], prim["ci_low"], prim["ci_high"])) return results if __name__ == "__main__": main()