etf_constituent_beta_forecasts.py

"""Forecast constituent betas with estimator and sector-shrinkage comparisons."""

from __future__ import annotations

import argparse
import math
import sys
import time
from pathlib import Path

import numpy as np
import pandas as pd

from beta_analysis import rolling_beta_estimates, rolling_ewma_beta_estimates
from etf_data import (
    DEFAULT_YFINANCE_CODE_DIR,
    daily_returns,
    ensure_saved_history,
    read_equity_holdings,
    sample_date_range,
    yahoo_symbol,
)
from forecast_evaluation import (
    add_common_configuration_flag,
    leave_one_out_group_priors,
)
from volatility_analysis import ewma_parameters


def positive_int(value: str) -> int:
    parsed = int(value)
    if parsed <= 0:
        raise argparse.ArgumentTypeError("must be positive")
    return parsed


def beta_window(value: str) -> int:
    parsed = positive_int(value)
    if parsed < 3:
        raise argparse.ArgumentTypeError("beta windows must be at least 3")
    return parsed


def unit_float(value: str) -> float:
    parsed = float(value)
    if not 0 <= parsed <= 1:
        raise argparse.ArgumentTypeError("must be between 0 and 1")
    return parsed


def parse_args() -> argparse.Namespace:
    parser = argparse.ArgumentParser(
        description=(
            "Compare ordinary and exponentially weighted constituent betas, "
            "with optional sector shrinkage, as forecasts of future beta."
        )
    )
    parser.add_argument("holdings_file", type=Path, help="ETF holdings CSV")
    parser.add_argument(
        "--benchmark", default="IVV", help="Beta benchmark (default: IVV)"
    )
    parser.add_argument(
        "--lookbacks",
        type=beta_window,
        nargs="+",
        default=[21, 63, 126, 252, 504, 756, 1260],
        help="Past-return windows used to forecast beta",
    )
    parser.add_argument(
        "--forecast-horizons",
        "--horizons",
        type=beta_window,
        nargs="+",
        default=[126],
        help="One or more future beta horizons (default: 126)",
    )
    parser.add_argument(
        "--beta-estimators",
        choices=["ols", "ewma"],
        nargs="+",
        default=["ols"],
        help=(
            "Historical beta estimators to compare; future realized beta "
            "always uses ordinary OLS (default: ols)"
        ),
    )
    parser.add_argument(
        "--fixed-stock-weights",
        type=unit_float,
        nargs="*",
        default=[0.25, 0.50, 0.75],
        metavar="WEIGHT",
        help=(
            "Fixed weights on raw stock beta in stock/sector blends; pass the "
            "option with no values to omit fixed blends"
        ),
    )
    parser.add_argument(
        "--min-sector-peers",
        type=positive_int,
        default=3,
        help=(
            "Minimum other stocks for a sector prior; otherwise use all "
            "available stocks (default: 3)"
        ),
    )
    parser.add_argument(
        "--coverage",
        type=unit_float,
        default=0.90,
        help="Required return coverage within every beta window (default: 0.90)",
    )
    parser.add_argument(
        "--history-days",
        type=positive_int,
        default=10000,
        help="Requested return history; all available data are used (default: 10000)",
    )
    parser.add_argument(
        "--annualization-days",
        type=float,
        default=252.0,
        help="Trading days used to annualize hedge volatility (default: 252)",
    )
    parser.add_argument(
        "--max-symbols",
        "--max-sym",
        type=positive_int,
        metavar="N",
        help="Use only the first N holdings rows for a faster test",
    )
    parser.add_argument(
        "--summary-only",
        action="store_true",
        help=(
            "Skip the large detailed forecast CSV; still write overall, sector, "
            "stock, and best-lookback summaries"
        ),
    )
    parser.add_argument("--name", help="Short filename prefix for result CSVs")
    parser.add_argument("--prices-file", type=Path, help="Adjusted-close cache CSV")
    parser.add_argument(
        "--no-download-missing",
        action="store_true",
        help="Fail instead of downloading missing price history",
    )
    parser.add_argument(
        "--output-dir", type=Path, default=Path("output"), help="Output directory"
    )
    parser.add_argument(
        "--code-dir",
        type=Path,
        default=DEFAULT_YFINANCE_CODE_DIR,
        help="Directory containing yfinance_util.py",
    )
    args = parser.parse_args()
    args.benchmark = yahoo_symbol(args.benchmark)
    args.lookbacks = sorted(set(args.lookbacks))
    args.forecast_horizons = sorted(set(args.forecast_horizons))
    args.beta_estimators = list(dict.fromkeys(args.beta_estimators))
    args.fixed_stock_weights = sorted(set(args.fixed_stock_weights))
    if args.coverage <= 0:
        parser.error("--coverage must be greater than zero")
    if args.annualization_days <= 0:
        parser.error("--annualization-days must be positive")
    return args


def fixed_method(weight: float) -> str:
    text = f"{weight:.6g}".replace(".", "p")
    return f"fixed_stock_weight_{text}"


def method_definitions(
    fixed_weights: list[float],
) -> list[tuple[str, str, str, str | None]]:
    methods = [
        ("beta_one", "beta_one_forecast", "beta_one_residual_volatility", None),
        ("raw", "raw_beta", "raw_residual_volatility", None),
        ("sector", "peer_mean_beta", "sector_residual_volatility", None),
        (
            "empirical_shrinkage",
            "empirical_beta",
            "empirical_shrinkage_residual_volatility",
            "empirical_stock_weight",
        ),
    ]
    for weight in fixed_weights:
        method = fixed_method(weight)
        methods.append(
            (
                method,
                f"{method}_beta",
                f"{method}_residual_volatility",
                None,
            )
        )
    return methods


def load_holdings(args: argparse.Namespace) -> pd.DataFrame:
    holdings = read_equity_holdings(args.holdings_file)
    if args.max_symbols is not None:
        holdings = holdings.head(args.max_symbols)
    holdings = holdings[holdings["yahoo_symbol"] != args.benchmark].copy()
    if holdings.empty:
        raise ValueError("holdings file has no constituents other than the benchmark")
    holdings = holdings.rename(columns={"ticker": "display_ticker"})
    return holdings.set_index("yahoo_symbol", drop=False)


def leave_one_out_priors(
    raw_betas: pd.Series,
    standard_errors: pd.Series,
    sectors: pd.Series,
    min_sector_peers: int,
) -> pd.DataFrame:
    priors = leave_one_out_group_priors(
        raw_betas,
        standard_errors.pow(2),
        sectors,
        min_sector_peers,
    )
    priors = priors.rename(
        columns={
            "peer_mean": "peer_mean_beta",
            "peer_observed_variance": "peer_observed_beta_variance",
            "peer_mean_estimation_variance": "peer_mean_sampling_variance",
            "estimated_within_group_variance": (
                "estimated_within_group_beta_variance"
            ),
            "empirical_own_weight": "empirical_stock_weight",
            "empirical_estimate": "empirical_beta",
        }
    )
    priors.loc[priors["prior_source"].eq("group"), "prior_source"] = "sector"
    return priors


def residual_volatilities(
    future_assets: pd.DataFrame,
    future_benchmark: pd.Series,
    forecast_betas: pd.Series,
    annualization_days: float,
) -> pd.Series:
    residuals = future_assets.sub(
        future_benchmark.to_numpy()[:, None] * forecast_betas.to_numpy()[None, :],
        axis=None,
    )
    return residuals.std(axis=0, ddof=1) * math.sqrt(annualization_days)


def forecast_schedule(
    observations: int, max_lookback: int, horizons: list[int]
) -> dict[int, list[int]]:
    return {
        horizon: list(range(max_lookback - 1, observations - horizon, horizon))
        for horizon in horizons
    }


def calculate_forecasts(
    returns: pd.DataFrame,
    holdings: pd.DataFrame,
    benchmark: str,
    lookbacks: list[int],
    horizons: list[int],
    beta_estimators: list[str],
    fixed_weights: list[float],
    coverage: float,
    min_sector_peers: int,
    annualization_days: float,
) -> pd.DataFrame:
    assets = holdings.index.tolist()
    benchmark_dates = returns.index[returns[benchmark].notna()]
    sample = returns.reindex(benchmark_dates)[[*assets, benchmark]]
    schedule = forecast_schedule(len(sample), max(lookbacks), horizons)
    if not any(schedule.values()):
        raise ValueError(
            "insufficient benchmark history for the longest lookback plus a "
            "forecast horizon"
        )

    needed_positions = {
        position
        for horizon, origins in schedule.items()
        for origin in origins
        for position in (origin, origin + horizon)
    }
    needed_dates = sample.index[sorted(needed_positions)]
    ols_estimates: dict[
        int, tuple[pd.DataFrame, pd.DataFrame, pd.DataFrame]
    ] = {}
    ols_windows = set(horizons)
    if "ols" in beta_estimators:
        ols_windows.update(lookbacks)
    for window in sorted(ols_windows):
        required = max(3, math.ceil(window * coverage))
        betas, standard_errors, observations = rolling_beta_estimates(
            sample, assets, benchmark, window, required
        )
        ols_estimates[window] = (
            betas.reindex(needed_dates),
            standard_errors.reindex(needed_dates),
            observations.reindex(needed_dates),
        )
    historical_estimates = {
        ("ordinary_ols", lookback): ols_estimates[lookback]
        for lookback in lookbacks
        if "ols" in beta_estimators
    }
    if "ewma" in beta_estimators:
        for lookback in lookbacks:
            required = max(3, math.ceil(lookback * coverage))
            betas, standard_errors, observations = (
                rolling_ewma_beta_estimates(
                    sample, assets, benchmark, lookback, required
                )
            )
            historical_estimates[("ewma", lookback)] = (
                betas.reindex(needed_dates),
                standard_errors.reindex(needed_dates),
                observations.reindex(needed_dates),
            )

    frames = []
    sectors = holdings["sector"]
    annualizer = math.sqrt(annualization_days)
    for horizon, origins in schedule.items():
        future_betas, future_standard_errors, future_observations = ols_estimates[
            horizon
        ]
        for origin in origins:
            origin_date = sample.index[origin]
            future_end_position = origin + horizon
            future_end_date = sample.index[future_end_position]
            future_start_date = sample.index[origin + 1]
            future = sample.iloc[origin + 1 : future_end_position + 1]
            future_assets = future[assets]
            future_benchmark = future[benchmark]
            for estimator in beta_estimators:
                estimator_name = (
                    "ordinary_ols" if estimator == "ols" else "ewma"
                )
                for lookback in lookbacks:
                    raw_betas, raw_standard_errors, raw_observations = (
                        historical_estimates[(estimator_name, lookback)]
                    )
                    raw = raw_betas.loc[origin_date]
                    raw_se = raw_standard_errors.loc[origin_date]
                    priors = leave_one_out_priors(
                        raw, raw_se, sectors, min_sector_peers
                    )
                    actual = future_betas.loc[future_end_date]
                    eligible = (
                        raw.notna()
                        & raw_se.notna()
                        & actual.notna()
                        & priors["peer_mean_beta"].notna()
                        & priors["empirical_beta"].notna()
                    )
                    symbols = eligible.index[eligible].tolist()
                    if not symbols:
                        continue
                    predictions: dict[str, pd.Series] = {
                        "beta_one": pd.Series(1.0, index=symbols),
                        "raw": raw[symbols],
                        "sector": priors.loc[symbols, "peer_mean_beta"],
                        "empirical_shrinkage": priors.loc[
                            symbols, "empirical_beta"
                        ],
                    }
                    for weight in fixed_weights:
                        predictions[fixed_method(weight)] = (
                            weight * raw[symbols]
                            + (1.0 - weight)
                            * priors.loc[symbols, "peer_mean_beta"]
                        )

                    result = pd.DataFrame(index=symbols)
                    result.index.name = "symbol"
                    result["display_ticker"] = holdings.loc[
                        symbols, "display_ticker"
                    ]
                    result["sector"] = holdings.loc[symbols, "sector"]
                    result["holdings_weight"] = holdings.loc[
                        symbols, "holdings_weight"
                    ]
                    result["benchmark"] = benchmark
                    result["beta_estimator"] = estimator_name
                    if estimator == "ewma":
                        alpha, decay, half_life = ewma_parameters(lookback)
                        estimation_start = sample.index[0]
                        result["estimator_parameter"] = decay
                        result["ewma_alpha"] = alpha
                        result["ewma_decay"] = decay
                        result["ewma_half_life"] = half_life
                    else:
                        estimation_start = sample.index[origin - lookback + 1]
                        result["estimator_parameter"] = np.nan
                        result["ewma_alpha"] = np.nan
                        result["ewma_decay"] = np.nan
                        result["ewma_half_life"] = np.nan
                    result["lookback_days"] = lookback
                    result["forecast_horizon_days"] = horizon
                    result["annualization_days"] = annualization_days
                    result["forecast_origin"] = origin_date.date().isoformat()
                    result["estimation_start_date"] = (
                        estimation_start.date().isoformat()
                    )
                    result["estimation_end_date"] = origin_date.date().isoformat()
                    result["future_start_date"] = (
                        future_start_date.date().isoformat()
                    )
                    result["future_end_date"] = future_end_date.date().isoformat()
                    result["estimation_observations"] = raw_observations.loc[
                        origin_date, symbols
                    ]
                    result["future_observations"] = future_observations.loc[
                        future_end_date, symbols
                    ]
                    result["raw_beta"] = raw[symbols]
                    result["raw_beta_standard_error"] = raw_se[symbols]
                    result["future_beta"] = actual[symbols]
                    result["future_beta_standard_error"] = (
                        future_standard_errors.loc[future_end_date, symbols]
                    )
                    for column in priors.columns:
                        result[column] = priors.loc[symbols, column]

                    for method, forecasts in predictions.items():
                        if method == "raw":
                            forecast_column = "raw_beta"
                        elif method == "sector":
                            forecast_column = "peer_mean_beta"
                        elif method == "empirical_shrinkage":
                            forecast_column = "empirical_beta"
                        elif method == "beta_one":
                            forecast_column = "beta_one_forecast"
                            result[forecast_column] = 1.0
                        else:
                            forecast_column = f"{method}_beta"
                            result[forecast_column] = forecasts
                        error_column = f"{method}_error"
                        result[error_column] = (
                            result[forecast_column] - result["future_beta"]
                        )
                        result[f"{method}_squared_error"] = result[
                            error_column
                        ].pow(2)
                        residual_column = f"{method}_residual_volatility"
                        result[residual_column] = residual_volatilities(
                            future_assets[symbols],
                            future_benchmark,
                            forecasts,
                            annualization_days,
                        )
                    oracle_residuals = future_assets[symbols].sub(
                        future_benchmark.to_numpy()[:, None]
                        * actual[symbols].to_numpy()[None, :],
                        axis=None,
                    )
                    result["oracle_residual_volatility"] = (
                        oracle_residuals.std(axis=0, ddof=1) * annualizer
                    )
                    frames.append(result.reset_index())
    if not frames:
        raise ValueError(
            "no forecasts have enough stock, benchmark, and peer-sector data"
        )
    return pd.concat(frames, ignore_index=True)


def summarize_method(
    group: pd.DataFrame,
    method: str,
    forecast_column: str,
    residual_column: str,
    weight_column: str | None,
) -> dict[str, object]:
    valid = group.dropna(subset=[forecast_column, "future_beta", residual_column])
    errors = valid[forecast_column] - valid["future_beta"]
    squared_errors = errors.pow(2)
    raw_squared_errors = (
        valid["raw_beta"] - valid["future_beta"]
    ).pow(2)
    if weight_column is not None:
        mean_stock_weight = float(valid[weight_column].mean())
    elif method == "raw":
        mean_stock_weight = 1.0
    elif method == "sector":
        mean_stock_weight = 0.0
    elif method.startswith("fixed_stock_weight_"):
        mean_stock_weight = float(method.rsplit("_", 1)[1].replace("p", "."))
    else:
        mean_stock_weight = float("nan")
    forecast_correlation = (
        valid[forecast_column].corr(valid["future_beta"])
        if valid[forecast_column].nunique() > 1
        and valid["future_beta"].nunique() > 1
        else float("nan")
    )
    return {
        "method": method,
        "observations": len(valid),
        "first_forecast_origin": valid["forecast_origin"].min(),
        "last_forecast_origin": valid["forecast_origin"].max(),
        "first_future_start_date": valid["future_start_date"].min(),
        "last_future_end_date": valid["future_end_date"].max(),
        "mean_forecast_beta": float(valid[forecast_column].mean()),
        "mean_future_beta": float(valid["future_beta"].mean()),
        "bias": float(errors.mean()),
        "mae": float(errors.abs().mean()),
        "rmse": math.sqrt(float(squared_errors.mean())),
        "forecast_future_correlation": float(forecast_correlation),
        "beats_raw_rate": (
            float((squared_errors < raw_squared_errors).mean())
            if method != "raw"
            else float("nan")
        ),
        "mean_stock_weight": mean_stock_weight,
        "mean_residual_annualized_volatility": float(
            valid[residual_column].mean()
        ),
        "mean_oracle_residual_annualized_volatility": float(
            valid["oracle_residual_volatility"].mean()
        ),
    }


def accuracy_table(
    forecasts: pd.DataFrame,
    group_columns: list[str],
    methods: list[tuple[str, str, str, str | None]],
) -> pd.DataFrame:
    rows = []
    for keys, group in forecasts.groupby(group_columns, sort=True, dropna=False):
        if not isinstance(keys, tuple):
            keys = (keys,)
        base = dict(zip(group_columns, keys))
        for method, forecast_column, residual_column, weight_column in methods:
            rows.append(
                {
                    **base,
                    "benchmark": group["benchmark"].iloc[0],
                    "beta_estimator": group["beta_estimator"].iloc[0],
                    "estimator_parameter": group["estimator_parameter"].iloc[0],
                    "annualization_days": group["annualization_days"].iloc[0],
                    **summarize_method(
                        group,
                        method,
                        forecast_column,
                        residual_column,
                        weight_column,
                    ),
                }
            )
    return pd.DataFrame(rows)


def best_lookbacks(overall: pd.DataFrame) -> pd.DataFrame:
    rows = []
    beta_one = overall[overall["method"].eq("beta_one")]
    for horizon, group in beta_one.groupby("forecast_horizon_days", sort=True):
        best = group.iloc[0]
        rows.append(
            {
                "benchmark": best["benchmark"],
                "beta_estimator": "not_applicable",
                "estimator_parameter": np.nan,
                "forecast_horizon_days": horizon,
                "method": "beta_one",
                "best_lookback_days": np.nan,
                "rmse": best["rmse"],
                "mae": best["mae"],
                "bias": best["bias"],
                "forecast_future_correlation": best[
                    "forecast_future_correlation"
                ],
                "beats_raw_rate": np.nan,
                "mean_stock_weight": best["mean_stock_weight"],
                "mean_residual_annualized_volatility": best[
                    "mean_residual_annualized_volatility"
                ],
                "observations": best["observations"],
                "first_future_start_date": best["first_future_start_date"],
                "last_future_end_date": best["last_future_end_date"],
            }
        )
    compared = overall[~overall["method"].eq("beta_one")]
    for (horizon, estimator, method), group in compared.groupby(
        ["forecast_horizon_days", "beta_estimator", "method"], sort=True
    ):
        best = group.loc[group["rmse"].idxmin()]
        rows.append(
            {
                "benchmark": best["benchmark"],
                "beta_estimator": estimator,
                "estimator_parameter": best["estimator_parameter"],
                "forecast_horizon_days": horizon,
                "method": method,
                "best_lookback_days": int(best["lookback_days"]),
                "rmse": best["rmse"],
                "mae": best["mae"],
                "bias": best["bias"],
                "forecast_future_correlation": best[
                    "forecast_future_correlation"
                ],
                "beats_raw_rate": best["beats_raw_rate"],
                "mean_stock_weight": best["mean_stock_weight"],
                "mean_residual_annualized_volatility": best[
                    "mean_residual_annualized_volatility"
                ],
                "observations": best["observations"],
                "first_future_start_date": best["first_future_start_date"],
                "last_future_end_date": best["last_future_end_date"],
            }
        )
    return pd.DataFrame(rows)


def print_best(best: pd.DataFrame) -> None:
    for horizon in sorted(best["forecast_horizon_days"].unique()):
        print(f"\nBest lookback by method for {horizon}-day future beta:")
        view = best[best["forecast_horizon_days"] == horizon][
            [
                "method",
                "beta_estimator",
                "best_lookback_days",
                "rmse",
                "mae",
                "bias",
                "forecast_future_correlation",
                "beats_raw_rate",
                "mean_stock_weight",
                "mean_residual_annualized_volatility",
                "observations",
            ]
        ].copy()
        for column in [
            "rmse",
            "mae",
            "bias",
            "forecast_future_correlation",
            "mean_stock_weight",
        ]:
            view[column] = view[column].map(
                lambda value: f"{value:.3f}" if pd.notna(value) else ""
            )
        view["beats_raw_rate"] = view["beats_raw_rate"].map(
            lambda value: f"{100 * value:.2f}%" if pd.notna(value) else ""
        )
        view["mean_residual_annualized_volatility"] = view[
            "mean_residual_annualized_volatility"
        ].map(lambda value: f"{100 * value:.2f}%")
        view["best_lookback_days"] = view["best_lookback_days"].map(
            lambda value: str(int(value)) if pd.notna(value) else "--"
        )
        print(view.to_string(index=False))


def main() -> int:
    overall_start = time.perf_counter()
    args = parse_args()
    holdings = load_holdings(args)
    args.output_dir.mkdir(parents=True, exist_ok=True)
    stem = args.name or args.holdings_file.stem
    prices_path = args.prices_file or args.output_dir / f"{stem}_adjusted_close.csv"
    symbols = holdings.index.tolist()
    data_symbols = [*symbols, args.benchmark]

    data_start = time.perf_counter()
    prices, downloaded, unavailable, download_elapsed = ensure_saved_history(
        prices_path,
        data_symbols,
        args.history_days,
        args.code_dir,
        args.no_download_missing,
        critical_symbols=[args.benchmark],
    )
    selected_prices = prices.reindex(columns=data_symbols)
    returns = daily_returns(selected_prices)
    data_elapsed = time.perf_counter() - data_start

    calculation_start = time.perf_counter()
    forecasts = calculate_forecasts(
        returns,
        holdings,
        args.benchmark,
        args.lookbacks,
        args.forecast_horizons,
        args.beta_estimators,
        args.fixed_stock_weights,
        args.coverage,
        args.min_sector_peers,
        args.annualization_days,
    )
    estimator_names = [
        "ordinary_ols" if estimator == "ols" else "ewma"
        for estimator in args.beta_estimators
    ]
    forecasts = add_common_configuration_flag(
        forecasts,
        [
            (estimator, lookback)
            for estimator in estimator_names
            for lookback in args.lookbacks
        ],
        ["beta_estimator", "lookback_days"],
        ["symbol", "forecast_origin", "future_start_date", "future_end_date"],
        flag_column="common_estimator_lookback_sample",
    )
    forecasts["common_lookback_sample"] = forecasts[
        "common_estimator_lookback_sample"
    ]
    common_forecasts = forecasts[
        forecasts["common_estimator_lookback_sample"]
    ].copy()
    if common_forecasts.empty:
        raise ValueError(
            "no forecast observations are available for every requested beta "
            "estimator and lookback"
        )
    methods = method_definitions(args.fixed_stock_weights)
    overall_accuracy = accuracy_table(
        common_forecasts,
        ["beta_estimator", "forecast_horizon_days", "lookback_days"],
        methods,
    )
    sector_accuracy = accuracy_table(
        common_forecasts,
        [
            "sector",
            "beta_estimator",
            "forecast_horizon_days",
            "lookback_days",
        ],
        methods,
    )
    stock_accuracy = accuracy_table(
        common_forecasts,
        [
            "symbol",
            "sector",
            "beta_estimator",
            "forecast_horizon_days",
            "lookback_days",
        ],
        methods,
    )
    best = best_lookbacks(overall_accuracy)
    paths = {
        "forecasts": args.output_dir / f"{stem}_constituent_beta_forecasts.csv",
        "overall": args.output_dir / f"{stem}_beta_forecast_accuracy.csv",
        "sector": args.output_dir / f"{stem}_sector_beta_forecast_accuracy.csv",
        "stock": args.output_dir / f"{stem}_stock_beta_forecast_accuracy.csv",
        "best": args.output_dir / f"{stem}_beta_forecast_best_lookbacks.csv",
    }
    written_paths = []
    if not args.summary_only:
        forecasts.to_csv(paths["forecasts"], index=False)
        written_paths.append(paths["forecasts"])
    overall_accuracy.to_csv(paths["overall"], index=False)
    written_paths.append(paths["overall"])
    sector_accuracy.to_csv(paths["sector"], index=False)
    written_paths.append(paths["sector"])
    stock_accuracy.to_csv(paths["stock"], index=False)
    written_paths.append(paths["stock"])
    best.to_csv(paths["best"], index=False)
    written_paths.append(paths["best"])
    calculation_elapsed = time.perf_counter() - calculation_start

    price_range = sample_date_range(selected_prices[args.benchmark].dropna())
    print(f"Using prices: {prices_path}")
    print(
        f"Benchmark price range: {price_range[0]} to {price_range[1]}; "
        f"requested up to {args.history_days} return days"
    )
    print(
        f"Holdings: {len(holdings)} stocks in "
        f"{holdings['sector'].nunique()} sectors; benchmark: {args.benchmark}"
    )
    print("Historical beta estimators: " + ", ".join(estimator_names))
    excluded_rows = len(forecasts) - len(common_forecasts)
    print(
        f"Common estimator/lookback evaluation rows: {len(common_forecasts)}; "
        f"excluded from summaries: {excluded_rows}"
    )
    available_horizons = set(common_forecasts["forecast_horizon_days"].unique())
    missing_horizons = [
        horizon
        for horizon in args.forecast_horizons
        if horizon not in available_horizons
    ]
    if missing_horizons:
        print(
            "No completed forecasts for horizons: "
            + ", ".join(str(horizon) for horizon in missing_horizons)
        )
    if downloaded:
        print("Downloaded symbols: " + ", ".join(downloaded))
    if unavailable:
        print("Unavailable symbols excluded: " + ", ".join(unavailable))
    print_best(best)
    if args.summary_only:
        print(
            "Skipped detailed forecast CSV (--summary-only): "
            f"{paths['forecasts']} (any existing file was not updated)"
        )
    for path in written_paths:
        print(f"Wrote {path}")
    print(f"Data elapsed: {data_elapsed:.3f} seconds")
    if downloaded:
        print(f"Download elapsed: {download_elapsed:.3f} seconds")
    print(f"Calculations elapsed: {calculation_elapsed:.3f} seconds")
    print(f"Overall elapsed: {time.perf_counter() - overall_start:.3f} seconds")
    return 0


if __name__ == "__main__":
    try:
        raise SystemExit(main())
    except (FileNotFoundError, ImportError, ValueError, RuntimeError) as exc:
        print(f"error: {exc}", file=sys.stderr)
        raise SystemExit(1)