diff --git a/src/models/causal_forest.py b/src/models/causal_forest.py
new file mode 100644
index 0000000..a7f2a17
--- /dev/null
+++ b/src/models/causal_forest.py
@@ -0,0 +1,1340 @@
+"""Causal forest for transfer effect estimation in Fantacalcio.
+
+Estimates the causal effect of roster changes on team performance using
+Generalized Random Forest (Athey et al., 2019). Separates a player's own
+contribution from confounders: replacement effects, synergies, redundancies,
+and concurrent roster changes.
+
+Two backends:
+- econml (preferred): econml.grf.CausalForest
+- Pure sklearn fallback: honest causal forest with treatment-effect splits
+"""
+
+import logging
+from dataclasses import dataclass
+from typing import Dict, List, Optional, Tuple
+
+import numpy as np
+import pandas as pd
+from sklearn.ensemble import RandomForestRegressor
+from sklearn.preprocessing import StandardScaler
+
+from .base_model import BaseModel
+
+logger = logging.getLogger(__name__)
+
+_VALID_ROLES = {"P", "D", "C", "A"}
+
+# ---------------------------------------------------------------------------
+# Feature extraction helpers
+# ---------------------------------------------------------------------------
+
+_ROLE_KEYWORDS: Dict[str, List[str]] = {
+    "P": ["portiere", "goalkeeper", "gk", "P_"],
+    "D": ["difensore", "defender", "df", "D_"],
+    "C": ["centrocampista", "midfielder", "mid", "mf", "C_"],
+    "A": ["attaccante", "forward", "striker", "fw", "A_"],
+}
+
+
+def _detect_role(col: str) -> Optional[str]:
+    lower = col.lower()
+    for role, keywords in _ROLE_KEYWORDS.items():
+        for kw in keywords:
+            if kw.lower() in lower:
+                return role
+    return None
+
+
+def _guess_role_column(df: pd.DataFrame) -> Optional[str]:
+    for col in df.columns:
+        if col.lower() in ("role", "ruolo", "position", "pos"):
+            return col
+    return None
+
+
+def _safe_float_series(col: pd.Series) -> np.ndarray:
+    return pd.to_numeric(col, errors="coerce").fillna(0).values
+
+
+# ---------------------------------------------------------------------------
+# Roster feature builder
+# ---------------------------------------------------------------------------
+
+def _build_roster_features(roster_df: pd.DataFrame) -> pd.Series:
+    """Compute roster-level summary features from a squad dataframe.
+
+    Parameters
+    ----------
+    roster_df : pd.DataFrame
+        One row per player currently on the roster. Expected columns:
+        role (or auto-detected), points, minutes, age, team, cost.
+
+    Returns
+    -------
+    pd.Series
+        Single-row series of roster-level features.
+    """
+    numeric_cols = roster_df.select_dtypes(include=[np.number]).columns.tolist()
+    role_col = _guess_role_column(roster_df)
+    n = len(roster_df)
+
+    features: Dict[str, float] = {}
+
+    # Per-role counts
+    if role_col is not None:
+        role_counts = roster_df[role_col].value_counts()
+        for r in _VALID_ROLES:
+            features[f"count_{r}"] = float(role_counts.get(r, 0))
+    else:
+        logger.warning("No role column found in roster_df; per-role features skipped")
+        for r in _VALID_ROLES:
+            features[f"count_{r}"] = 0.0
+
+    # Per-role average points
+    if role_col is not None and "points" in numeric_cols:
+        for r in _VALID_ROLES:
+            mask = roster_df[role_col] == r
+            if mask.sum() > 0:
+                features[f"avg_points_{r}"] = float(roster_df.loc[mask, "points"].mean())
+            else:
+                features[f"avg_points_{r}"] = 0.0
+
+    # Per-role average minutes
+    if role_col is not None and "minutes" in numeric_cols:
+        for r in _VALID_ROLES:
+            mask = roster_df[role_col] == r
+            if mask.sum() > 0:
+                features[f"avg_minutes_{r}"] = float(roster_df.loc[mask, "minutes"].mean())
+            else:
+                features[f"avg_minutes_{r}"] = 0.0
+
+    # Total points
+    if "points" in numeric_cols:
+        features["total_points"] = float(roster_df["points"].sum())
+        features["points_std"] = float(roster_df["points"].std()) if n > 1 else 0.0
+    else:
+        features["total_points"] = 0.0
+        features["points_std"] = 0.0
+
+    # Age
+    if "age" in numeric_cols:
+        features["age_avg"] = float(roster_df["age"].mean()) if n > 0 else 0.0
+        features["age_std"] = float(roster_df["age"].std()) if n > 1 else 0.0
+    else:
+        features["age_avg"] = 0.0
+        features["age_std"] = 0.0
+
+    # Formation entropy (diversity of roles)
+    if role_col is not None:
+        total_players = n if n > 0 else 1
+        entropy = 0.0
+        for r in _VALID_ROLES:
+            p = features.get(f"count_{r}", 0.0) / total_players
+            if p > 0:
+                entropy -= p * np.log(p)
+        features["formation_entropy"] = entropy
+    else:
+        features["formation_entropy"] = 0.0
+
+    # Budget / slots remaining (if available)
+    features["budget_remaining"] = float(roster_df.get("budget", 0)) if isinstance(roster_df.get("budget", (int, float)), (int, float, np.integer, np.floating)) else 0.0
+    # Try columns for budget
+    if "budget" in roster_df.columns:
+        features["budget_remaining"] = float(roster_df["budget"].iloc[0]) if len(roster_df) > 0 else 0.0
+
+    if "slots_remaining" in roster_df.columns:
+        features["slots_remaining"] = float(roster_df["slots_remaining"].iloc[0]) if len(roster_df) > 0 else 0.0
+    else:
+        features["slots_remaining"] = 0.0
+
+    # Interaction level (shared-team pairs)
+    if "team" in roster_df.columns:
+        teams = roster_df["team"].values
+        interaction = 0
+        for i in range(len(teams)):
+            for j in range(i + 1, len(teams)):
+                if teams[i] == teams[j]:
+                    interaction += 1
+        max_pairs = n * (n - 1) / 2 if n > 1 else 1
+        features["interaction_level"] = interaction / max(max_pairs, 1)
+    else:
+        features["interaction_level"] = 0.0
+
+    # Role scarcity: count of unique positions within each role
+    if role_col is not None and "position" in roster_df.columns:
+        for r in _VALID_ROLES:
+            mask = roster_df[role_col] == r
+            if mask.sum() > 0:
+                features[f"role_diversity_{r}"] = float(roster_df.loc[mask, "position"].nunique())
+            else:
+                features[f"role_diversity_{r}"] = 0.0
+
+    return pd.Series(features)
+
+
+# ---------------------------------------------------------------------------
+# Player / treatment feature builder
+# ---------------------------------------------------------------------------
+
+def _build_player_features(player_row: pd.Series, player_pool_df: Optional[pd.DataFrame] = None) -> pd.Series:
+    """Compute treatment features for a player being added or removed.
+
+    Parameters
+    ----------
+    player_row : pd.Series
+        One player's data. Expected keys: projected_points, role, team,
+        age, minutes_last_season, projected_points (or points), name.
+    player_pool_df : pd.DataFrame, optional
+        Full player pool for computing scarcity / VORP.
+
+    Returns
+    -------
+    pd.Series
+        Treatment feature vector.
+    """
+    feats: Dict[str, float] = {}
+
+    # Core player stats
+    feats["player_projected_points"] = float(player_row.get("projected_points", player_row.get("points", 0))) if any(k in player_row.index for k in ("projected_points", "points")) else 0.0
+    feats["player_age"] = float(player_row.get("age", 0)) if "age" in player_row.index else 0.0
+    feats["player_minutes_last_season"] = float(player_row.get("minutes_last_season", player_row.get("minutes", 0))) if any(k in player_row.index for k in ("minutes_last_season", "minutes")) else 0.0
+
+    # Role one-hot
+    role = str(player_row.get("role", player_row.get("ruolo", player_row.get("position", ""))))
+    for r in _VALID_ROLES:
+        feats[f"role_is_{r}"] = 1.0 if r == role else 0.0
+
+    feats["player_role"] = role
+
+    # Role scarcity: how many comparable players exist in the pool
+    if player_pool_df is not None and not player_pool_df.empty:
+        pool_role_col = _guess_role_column(player_pool_df)
+        query_role = role
+        if pool_role_col is not None and query_role:
+            same_role = player_pool_df[player_pool_df[pool_role_col] == query_role]
+            feats["role_scarcity"] = len(same_role)
+        else:
+            feats["role_scarcity"] = float(len(player_pool_df))
+
+        # Value over replacement
+        if "projected_points" in player_pool_df.columns or "points" in player_pool_df.columns:
+            pts_col = "projected_points" if "projected_points" in player_pool_df.columns else "points"
+            pool_pts = _safe_float_series(player_pool_df[pts_col])
+            if pool_role_col is not None and query_role:
+                same_role_mask = player_pool_df[pool_role_col] == query_role
+                if same_role_mask.sum() > 0:
+                    replacement_pts = np.quantile(pool_pts[same_role_mask], 0.25)
+                else:
+                    replacement_pts = np.quantile(pool_pts, 0.25)
+            else:
+                replacement_pts = np.quantile(pool_pts, 0.25)
+            feats["value_over_replacement"] = feats.get("player_projected_points", 0.0) - replacement_pts
+        else:
+            feats["value_over_replacement"] = 0.0
+    else:
+        feats["role_scarcity"] = 0.0
+        feats["value_over_replacement"] = 0.0
+
+    return pd.Series(feats)
+
+
+def _build_treatment_features(
+    players_df: pd.DataFrame,
+    player_pool_df: Optional[pd.DataFrame] = None,
+) -> pd.DataFrame:
+    """Build treatment feature matrix from a dataframe of transferred players.
+
+    Parameters
+    ----------
+    players_df : pd.DataFrame
+        One row per transferred player.
+    player_pool_df : pd.DataFrame, optional
+        Full pool for scarcity/VORP computation.
+
+    Returns
+    -------
+    pd.DataFrame
+        Treatment feature matrix (n_players x n_features).
+    """
+    records = []
+    for _, row in players_df.iterrows():
+        records.append(_build_player_features(row, player_pool_df))
+    return pd.DataFrame(records)
+
+
+# ---------------------------------------------------------------------------
+# Sklearn fallback: Honest Causal Forest
+# ---------------------------------------------------------------------------
+
+class _HonestCausalTree:
+    """Single honest causal tree.
+
+    Splits the tree-building half to maximize treatment effect heterogeneity
+    (causal criterion), then estimates treatment effects on the estimation half.
+    """
+
+    def __init__(
+        self,
+        min_samples_leaf: int = 5,
+        max_depth: int = 5,
+        random_state: Optional[int] = None,
+    ):
+        self.min_samples_leaf = min_samples_leaf
+        self.max_depth = max_depth
+        self.random_state = random_state
+        self._rng = np.random.RandomState(random_state)
+        self._split_col: Optional[int] = None
+        self._split_val: Optional[float] = None
+        self._left: Optional["_HonestCausalTree"] = None
+        self._right: Optional["_HonestCausalTree"] = None
+        self._tau: Optional[float] = None
+        self._is_leaf: bool = True
+
+    def _causal_criterion(
+        self, X: np.ndarray, T: np.ndarray, Y: np.ndarray, col: int, val: float
+    ) -> float:
+        """Compute the causal splitting criterion.
+
+        Split to maximize n_L * n_R / n^2 * (τ_L - τ_R)^2
+        where τ is the estimated treatment effect in each child node.
+        """
+        n = len(X)
+        if n < 2:
+            return -1.0
+
+        mask_left = X[:, col] <= val
+        mask_right = ~mask_left
+        n_left = mask_left.sum()
+        n_right = mask_right.sum()
+
+        if n_left < self.min_samples_leaf or n_right < self.min_samples_leaf:
+            return -1.0
+
+        tau_left = _estimate_leaf_effect(T[mask_left], Y[mask_left])
+        tau_right = _estimate_leaf_effect(T[mask_right], Y[mask_right])
+
+        if tau_left is None or tau_right is None:
+            return -1.0
+
+        weight = (n_left * n_right) / (n * n)
+        return weight * (tau_left - tau_right) ** 2
+
+    def fit(self, X: np.ndarray, T: np.ndarray, Y: np.ndarray, depth: int = 0):
+        """Grow the tree on the tree-building half."""
+        n = len(X)
+        if n < 2 * self.min_samples_leaf or depth >= self.max_depth:
+            self._tau = _estimate_leaf_effect(T, Y) or 0.0
+            self._is_leaf = True
+            return
+
+        best_score = -1.0
+        best_col = -1
+        best_val = 0.0
+
+        n_features = X.shape[1]
+        feature_subset = min(n_features, max(1, int(np.sqrt(n_features))))
+        candidate_cols = self._rng.choice(n_features, size=feature_subset, replace=False)
+
+        for col in candidate_cols:
+            col_vals = X[:, col]
+            unique_vals = np.unique(col_vals)
+            if len(unique_vals) <= 1:
+                continue
+            # Try a random subset of splits
+            n_trials = min(10, len(unique_vals) - 1)
+            trial_vals = self._rng.choice(unique_vals[:-1], size=n_trials, replace=False)
+            for val in trial_vals:
+                score = self._causal_criterion(X, T, Y, col, val)
+                if score > best_score:
+                    best_score = score
+                    best_col = col
+                    best_val = val
+
+        if best_score <= 0:
+            self._tau = _estimate_leaf_effect(T, Y) or 0.0
+            self._is_leaf = True
+            return
+
+        self._is_leaf = False
+        self._split_col = best_col
+        self._split_val = best_val
+
+        mask_left = X[:, best_col] <= best_val
+        mask_right = ~mask_left
+
+        self._left = _HonestCausalTree(
+            min_samples_leaf=self.min_samples_leaf,
+            max_depth=self.max_depth,
+            random_state=self._rng.randint(0, 2 ** 31),
+        )
+        self._right = _HonestCausalTree(
+            min_samples_leaf=self.min_samples_leaf,
+            max_depth=self.max_depth,
+            random_state=self._rng.randint(0, 2 ** 31),
+        )
+        self._left.fit(X[mask_left], T[mask_left], Y[mask_left], depth + 1)
+        self._right.fit(X[mask_right], T[mask_right], Y[mask_right], depth + 1)
+
+    def estimate(self, X: np.ndarray, T: np.ndarray, Y: np.ndarray):
+        """Assign estimation-half observations to leaves and compute τ per leaf."""
+        if self._is_leaf:
+            self._tau = _estimate_leaf_effect(T, Y) or 0.0
+            return
+
+        mask_left = X[:, self._split_col] <= self._split_val
+        mask_right = ~mask_left
+
+        if mask_left.sum() > 0:
+            self._left.estimate(X[mask_left], T[mask_left], Y[mask_left])
+        else:
+            self._left._tau = 0.0
+
+        if mask_right.sum() > 0:
+            self._right.estimate(X[mask_right], T[mask_right], Y[mask_right])
+        else:
+            self._right._tau = 0.0
+
+    def predict(self, X: np.ndarray) -> np.ndarray:
+        """Predict τ for each row in X."""
+        n = X.shape[0]
+        if n == 0:
+            return np.array([])
+        if self._is_leaf:
+            return np.full(n, self._tau or 0.0)
+
+        preds = np.empty(n)
+        mask_left = X[:, self._split_col] <= self._split_val
+        mask_right = ~mask_left
+
+        if mask_left.sum() > 0 and self._left is not None:
+            preds[mask_left] = self._left.predict(X[mask_left])
+        elif self._left is not None:
+            pass
+
+        if mask_right.sum() > 0 and self._right is not None:
+            preds[mask_right] = self._right.predict(X[mask_right])
+
+        return preds
+
+
+def _estimate_leaf_effect(T: np.ndarray, Y: np.ndarray) -> Optional[float]:
+    """Estimate τ = E[Y|T=1] - E[Y|T=0] for a leaf node.
+
+    If treatment is continuous, use covariance-based estimator:
+    τ = Cov(Y, T) / Var(T)
+
+    If treatment is binary, use difference in means.
+    """
+    if len(T) < 2:
+        return None
+    unique_t = np.unique(T)
+    if len(unique_t) <= 1:
+        return None
+
+    # Check if binary treatment
+    if set(unique_t).issubset({0.0, 1.0}) and len(unique_t) == 2:
+        t1_mask = T >= 0.5
+        t0_mask = ~t1_mask
+        n1 = t1_mask.sum()
+        n0 = t0_mask.sum()
+        if n1 < 1 or n0 < 1:
+            return None
+        return float(Y[t1_mask].mean() - Y[t0_mask].mean())
+
+    # Continuous treatment: use OLS slope
+    T_centered = T - T.mean()
+    denom = np.dot(T_centered, T_centered)
+    if denom < 1e-12:
+        return None
+    Y_centered = Y - Y.mean()
+    beta = np.dot(T_centered, Y_centered) / denom
+    return float(beta)
+
+
+class _SklearnCausalForest:
+    """Pure sklearn honest causal forest.
+
+    Grows multiple _HonestCausalTree on bootstrapped splits of the data,
+    using honest estimation (split on one half, estimate on the other).
+    """
+
+    def __init__(
+        self,
+        n_estimators: int = 100,
+        min_samples_leaf: int = 5,
+        max_depth: int = 5,
+        random_state: int = 42,
+    ):
+        self.n_estimators = n_estimators
+        self.min_samples_leaf = min_samples_leaf
+        self.max_depth = max_depth
+        self.random_state = random_state
+        self._rng = np.random.RandomState(random_state)
+        self._trees: List[_HonestCausalTree] = []
+        self._fitted = False
+
+    def fit(self, X: np.ndarray, T: np.ndarray, Y: np.ndarray):
+        X = np.asarray(X, dtype=np.float64)
+        T = np.asarray(T, dtype=np.float64).ravel()
+        Y = np.asarray(Y, dtype=np.float64).ravel()
+
+        n = len(X)
+        self._trees = []
+
+        for i in range(self.n_estimators):
+            # Bootstrap sample
+            idx = self._rng.randint(0, n, size=n)
+            X_boot = X[idx]
+            T_boot = T[idx]
+            Y_boot = Y[idx]
+
+            # Honest split: 50% for tree building, 50% for estimation
+            n_boot = len(X_boot)
+            n_build = n_boot // 2
+            if n_build < 2 * self.min_samples_leaf:
+                continue
+
+            idx_build = self._rng.choice(n_boot, size=n_build, replace=False)
+            idx_est = np.setdiff1d(np.arange(n_boot), idx_build)
+            if len(idx_est) < 2 * self.min_samples_leaf:
+                continue
+
+            tree = _HonestCausalTree(
+                min_samples_leaf=self.min_samples_leaf,
+                max_depth=self.max_depth,
+                random_state=self._rng.randint(0, 2 ** 31),
+            )
+            # Split on build half
+            tree.fit(X_boot[idx_build], T_boot[idx_build], Y_boot[idx_build])
+            # Estimate leaf effects on estimation half
+            tree.estimate(X_boot[idx_est], T_boot[idx_est], Y_boot[idx_est])
+
+            self._trees.append(tree)
+
+        self._fitted = True
+        logger.info(
+            f"SklearnCausalForest fitted: {len(self._trees)}/{self.n_estimators} trees "
+            f"(min_samples_leaf={self.min_samples_leaf}, max_depth={self.max_depth})"
+        )
+
+    def predict(self, X: np.ndarray) -> Tuple[np.ndarray, np.ndarray]:
+        """Return (point predictions, variance across trees)."""
+        if not self._fitted:
+            raise RuntimeError("Model not fitted.")
+        X = np.asarray(X, dtype=np.float64)
+        n_trees = len(self._trees)
+        if n_trees == 0:
+            return np.zeros(X.shape[0]), np.zeros(X.shape[0])
+
+        all_preds = np.zeros((n_trees, X.shape[0]))
+        for i, tree in enumerate(self._trees):
+            all_preds[i] = tree.predict(X)
+
+        means = all_preds.mean(axis=0)
+        vars_ = all_preds.var(axis=0, ddof=1) if n_trees > 1 else np.zeros_like(means)
+        return means, vars_
+
+
+# ---------------------------------------------------------------------------
+# EconML dependency check
+# ---------------------------------------------------------------------------
+
+def _has_econml() -> bool:
+    try:
+        import econml.grf  # noqa: F401
+        return True
+    except ImportError:
+        return False
+
+
+# ---------------------------------------------------------------------------
+# TransferCausalModel
+# ---------------------------------------------------------------------------
+
+class TransferCausalModel(BaseModel):
+    """Causal forest for estimating transfer effects on team performance.
+
+    Estimates the causal effect of adding or removing a player from a Fantacalcio
+    roster, separating the player's own contribution from confounders
+    (replacement effects, synergies, concurrent roster moves).
+
+    Parameters
+    ----------
+    model_dir : str
+        Directory for persisting trained models.
+    n_estimators : int
+        Number of trees in the causal forest.
+    min_samples_leaf : int
+        Minimum samples per leaf.
+    max_depth : int
+        Maximum tree depth.
+    random_state : int
+        Random seed for reproducibility.
+    """
+
+    def __init__(
+        self,
+        model_dir: str = "models_trained",
+        n_estimators: int = 100,
+        min_samples_leaf: int = 5,
+        max_depth: int = 5,
+        random_state: int = 42,
+    ):
+        super().__init__(model_dir)
+        self.n_estimators = n_estimators
+        self.min_samples_leaf = min_samples_leaf
+        self.max_depth = max_depth
+        self.random_state = random_state
+
+        self._model = None
+        self._use_econml = False
+        self._scaler = StandardScaler()
+        self._feature_names: List[str] = []
+        self._treatment_names: List[str] = []
+        self._fitted = False
+
+    # ------------------------------------------------------------------
+    # Fit
+    # ------------------------------------------------------------------
+
+    def fit(
+        self,
+        X: pd.DataFrame,
+        treatment: pd.DataFrame,
+        outcome: pd.Series,
+        **kwargs,
+    ):
+        """Fit the causal forest.
+
+        Parameters
+        ----------
+        X : pd.DataFrame
+            Pre-transfer roster features (one row per transfer event).
+        treatment : pd.DataFrame
+            Player features for the transferred player. Rows must align
+            with X. Can include binary indicators or continuous quality scores.
+        outcome : pd.Series
+            Post-transfer team score change (can be positive or negative).
+        """
+        if len(X) == 0:
+            raise ValueError("X cannot be empty")
+        if len(X) != len(treatment):
+            raise ValueError(
+                f"X and treatment must have same length: {len(X)} vs {len(treatment)}"
+            )
+        if len(X) != len(outcome):
+            raise ValueError(
+                f"X and outcome must have same length: {len(X)} vs {len(outcome)}"
+            )
+
+        self._feature_names = list(X.columns)
+        self._treatment_names = list(treatment.columns)
+
+        X_num = X.select_dtypes(include=[np.number]).fillna(0)
+        T_num = treatment.select_dtypes(include=[np.number]).fillna(0)
+        Y_num = np.asarray(outcome, dtype=np.float64).ravel()
+
+        X_scaled = self._scaler.fit_transform(X_num)
+
+        # Build a combined feature matrix: roster features + treatment features
+        W = np.hstack([X_scaled, np.asarray(T_num, dtype=np.float64)])
+
+        if _has_econml():
+            try:
+                self._fit_econml(W, T_num, Y_num)
+                self._use_econml = True
+                self._fitted = True
+                return self
+            except Exception as exc:
+                logger.warning(
+                    f"econml failed ({exc}); falling back to sklearn causal forest"
+                )
+
+        self._use_econml = False
+        self._fit_sklearn(W, T_num, Y_num)
+        self._fitted = True
+        return self
+
+    def _fit_econml(self, W: np.ndarray, T: pd.DataFrame, Y: np.ndarray):
+        from econml.grf import CausalForest
+
+        T_num = np.asarray(T, dtype=np.float64)
+        # Use the first treatment column as the primary treatment
+        # If the dataframe has multiple columns, combine them
+        if T_num.ndim == 2 and T_num.shape[1] > 1:
+            logger.info(
+                f"Multiple treatment columns ({T_num.shape[1]}). "
+                "Using first column as primary treatment."
+            )
+            T_primary = T_num[:, 0]
+        elif T_num.ndim == 2:
+            T_primary = T_num[:, 0]
+        else:
+            T_primary = T_num.ravel()
+
+        self._model = CausalForest(
+            n_estimators=self.n_estimators,
+            min_samples_leaf=self.min_samples_leaf,
+            max_depth=self.max_depth,
+            random_state=self.random_state,
+        )
+        self._model.fit(W, T_primary, Y)
+        logger.info(
+            f"TransferCausalModel fitted via econml (n={len(W)}, "
+            f"features={W.shape[1]})"
+        )
+
+    def _fit_sklearn(self, W: np.ndarray, T: pd.DataFrame, Y: np.ndarray):
+        T_num = np.asarray(T, dtype=np.float64)
+        if T_num.ndim == 2 and T_num.shape[1] > 1:
+            T_primary = T_num[:, 0]
+        elif T_num.ndim == 2:
+            T_primary = T_num[:, 0]
+        else:
+            T_primary = T_num.ravel()
+
+        self._model = _SklearnCausalForest(
+            n_estimators=self.n_estimators,
+            min_samples_leaf=self.min_samples_leaf,
+            max_depth=self.max_depth,
+            random_state=self.random_state,
+        )
+        self._model.fit(W, T_primary, Y)
+        logger.info(
+            f"TransferCausalModel fitted via sklearn causal forest "
+            f"(n={len(W)}, features={W.shape[1]})"
+        )
+
+    # ------------------------------------------------------------------
+    # Predict effect
+    # ------------------------------------------------------------------
+
+    def predict(self, X: pd.DataFrame) -> np.ndarray:
+        """BaseModel interface — returns point prediction of treatment effect."""
+        if not self._fitted:
+            raise RuntimeError("Model not fitted. Call fit() first.")
+        ate, cate, _, _ = self.predict_effect(X)
+        return cate
+
+    def predict_effect(
+        self,
+        X: pd.DataFrame,
+        treatment: Optional[pd.DataFrame] = None,
+    ) -> Tuple[float, np.ndarray, np.ndarray, np.ndarray]:
+        """Predict average and conditional treatment effects.
+
+        Parameters
+        ----------
+        X : pd.DataFrame
+            Roster feature matrix.
+        treatment : pd.DataFrame, optional
+            Treatment feature matrix. If None, effects are predicted at the
+            observed treatment levels from training.
+
+        Returns
+        -------
+        ate : float
+            Average treatment effect across all observations.
+        cate : np.ndarray
+            Conditional average treatment effect per observation.
+        ci_lower : np.ndarray
+            Lower 95% confidence bound per observation.
+        ci_upper : np.ndarray
+            Upper 95% confidence bound per observation.
+        """
+        if not self._fitted:
+            raise RuntimeError("Model not fitted. Call fit() first.")
+
+        X_num = X.select_dtypes(include=[np.number]).fillna(0)
+        X_scaled = self._scaler.transform(X_num)
+
+        if treatment is not None:
+            T_num = treatment.select_dtypes(include=[np.number]).fillna(0)
+            W = np.hstack([X_scaled, np.asarray(T_num, dtype=np.float64)])
+        else:
+            # Use zeros for treatment features — predicts baseline
+            n_treat = len(self._treatment_names) if self._treatment_names else 1
+            W = np.hstack([X_scaled, np.zeros((X_scaled.shape[0], n_treat))])
+
+        if self._use_econml and self._model is not None:
+            cate = self._model.predict(W)
+            cate = np.asarray(cate, dtype=np.float64).ravel()
+            # EconML provides confidence intervals via bootstrap
+            ate = float(cate.mean())
+            # Approximate CI via jackknife/variance
+            var_est = float(cate.var(ddof=1)) if len(cate) > 1 else 0.0
+            se = np.sqrt(max(var_est, 1e-10))
+            if len(cate) <= 1:
+                se_arr = np.full_like(cate, se)
+            else:
+                se_arr = _jackknife_se(cate)
+            ci_lower = cate - 1.96 * se_arr
+            ci_upper = cate + 1.96 * se_arr
+            return ate, cate, ci_lower, ci_upper
+
+        # sklearn fallback
+        if self._model is not None:
+            cate, var_est = self._model.predict(W)
+            cate = np.asarray(cate, dtype=np.float64).ravel()
+            var_est = np.asarray(var_est, dtype=np.float64).ravel()
+            ate = float(cate.mean())
+            se = np.sqrt(np.maximum(var_est, 1e-10))
+            ci_lower = cate - 1.96 * se
+            ci_upper = cate + 1.96 * se
+            return ate, cate, ci_lower, ci_upper
+
+        raise RuntimeError("No underlying model available.")
+
+    # ------------------------------------------------------------------
+    # Individual transfer prediction
+    # ------------------------------------------------------------------
+
+    def predict_individual_effect(
+        self,
+        current_roster_df: pd.DataFrame,
+        candidate_player: pd.Series,
+        role_to_replace: Optional[str] = None,
+        player_pool_df: Optional[pd.DataFrame] = None,
+    ) -> Dict:
+        """Predict the net effect of adding a candidate player to the roster.
+
+        Parameters
+        ----------
+        current_roster_df : pd.DataFrame
+            Current roster (one row per player).
+        candidate_player : pd.Series
+            Candidate player data.
+        role_to_replace : str, optional
+            If provided, simulate removing a player of this role first.
+        player_pool_df : pd.DataFrame, optional
+            Full player pool for scarcity/VORP computation.
+
+        Returns
+        -------
+        dict with keys:
+            effect, ci_lower, ci_upper, recommendation, confidence,
+            synergy_bonus, redundancy_penalty, base_projection
+        """
+        if not self._fitted:
+            raise RuntimeError("Model not fitted. Call fit() first.")
+
+        # Build roster features for the current state
+        roster_feats = _build_roster_features(current_roster_df)
+        X_df = pd.DataFrame([roster_feats])
+
+        # Build treatment features for the candidate
+        player_feats = _build_player_features(candidate_player, player_pool_df)
+        T_df = pd.DataFrame([player_feats])
+
+        # Align with training feature names
+        for col in self._feature_names:
+            if col not in X_df.columns:
+                X_df[col] = 0.0
+        X_df = X_df[self._feature_names]
+
+        for col in self._treatment_names:
+            if col not in T_df.columns:
+                T_df[col] = 0.0
+        T_df = T_df[self._treatment_names]
+
+        _, cate, ci_lower, ci_upper = self.predict_effect(X_df, T_df)
+        effect = float(cate[0])
+        ci_lo = float(ci_lower[0])
+        ci_hi = float(ci_upper[0])
+
+        # Base projection
+        base_proj = float(player_feats.get("player_projected_points", 0.0))
+
+        # Synergy / redundancy analysis
+        role_col = _guess_role_column(current_roster_df)
+        candidate_role = player_feats.get("player_role", "")
+        synergy_bonus = 0.0
+        redundancy_penalty = 0.0
+
+        if role_col is not None and candidate_role:
+            role_count = int((current_roster_df[role_col] == candidate_role).sum())
+            role_points = _safe_float_series(
+                current_roster_df.loc[current_roster_df[role_col] == candidate_role, "points"]
+                if "points" in current_roster_df.columns
+                else pd.Series([0.0])
+            ).mean()
+
+            # Synergy bonus: positive effect from interaction with teammates
+            if "interaction_level" in roster_feats.index:
+                synergy_bonus = float(roster_feats["interaction_level"]) * 0.5
+
+            # Redundancy penalty: already strong in this role
+            if role_count >= 4:
+                redundancy_penalty = min(float(role_points), 2.0)
+
+        # Confidence based on CI width
+        ci_width = ci_hi - ci_lo
+        if ci_width < 1.0:
+            confidence = "HIGH"
+        elif ci_width < 3.0:
+            confidence = "MEDIUM"
+        else:
+            confidence = "LOW"
+
+        # Recommendation
+        if effect > ci_width * 0.5 and effect > 0:
+            recommendation = "BUY"
+        elif effect < -ci_width * 0.5 and effect < 0:
+            recommendation = "SELL"
+        else:
+            recommendation = "AVOID"
+
+        return {
+            "effect": round(effect, 3),
+            "ci_lower": round(ci_lo, 3),
+            "ci_upper": round(ci_hi, 3),
+            "recommendation": recommendation,
+            "confidence": confidence,
+            "synergy_bonus": round(synergy_bonus, 3),
+            "redundancy_penalty": round(redundancy_penalty, 3),
+            "base_projection": round(base_proj, 2),
+        }
+
+    # ------------------------------------------------------------------
+    # Rank transfers
+    # ------------------------------------------------------------------
+
+    def rank_transfers(
+        self,
+        current_roster_df: pd.DataFrame,
+        candidate_pool_df: pd.DataFrame,
+        n_recommendations: int = 10,
+    ) -> pd.DataFrame:
+        """Rank all possible transfers by predicted causal effect.
+
+        Parameters
+        ----------
+        current_roster_df : pd.DataFrame
+            Current roster.
+        candidate_pool_df : pd.DataFrame
+            Pool of candidate players.
+        n_recommendations : int
+            Number of top recommendations to return.
+
+        Returns
+        -------
+        pd.DataFrame
+            Ranked transfer recommendations with columns:
+            player_name, role, current_team, base_projection, causal_effect,
+            ci_lower, ci_upper, synergy_bonus, redundancy_penalty,
+            recommendation, confidence
+        """
+        if not self._fitted:
+            raise RuntimeError("Model not fitted. Call fit() first.")
+
+        roster_feats = _build_roster_features(current_roster_df)
+        X_df = pd.DataFrame([roster_feats])
+
+        role_col = _guess_role_column(current_roster_df)
+
+        results = []
+        for idx, row in candidate_pool_df.iterrows():
+            player_feats = _build_player_features(row, candidate_pool_df)
+            T_df = pd.DataFrame([player_feats])
+
+            X_aligned = X_df.copy()
+            for col in self._feature_names:
+                if col not in X_aligned.columns:
+                    X_aligned[col] = 0.0
+            X_aligned = X_aligned[self._feature_names]
+
+            T_aligned = T_df.copy()
+            for col in self._treatment_names:
+                if col not in T_aligned.columns:
+                    T_aligned[col] = 0.0
+            T_aligned = T_aligned[self._treatment_names]
+
+            _, cate, ci_lower, ci_upper = self.predict_effect(X_aligned, T_aligned)
+            effect = float(cate[0])
+            ci_lo = float(ci_lower[0])
+            ci_hi = float(ci_upper[0])
+
+            base_proj = float(player_feats.get("player_projected_points", 0.0))
+            candidate_role = player_feats.get("player_role", "")
+            synergy_bonus = 0.0
+            redundancy_penalty = 0.0
+
+            if role_col is not None and candidate_role:
+                role_count = int((current_roster_df[role_col] == candidate_role).sum())
+                role_points = _safe_float_series(
+                    current_roster_df.loc[current_roster_df[role_col] == candidate_role, "points"]
+                    if "points" in current_roster_df.columns
+                    else pd.Series([0.0])
+                ).mean()
+                if role_count >= 4:
+                    redundancy_penalty = min(float(role_points), 2.0)
+                if "interaction_level" in roster_feats.index:
+                    synergy_bonus = float(roster_feats["interaction_level"]) * 0.5
+
+            ci_width = ci_hi - ci_lo
+            if ci_width < 1.0:
+                confidence = "HIGH"
+            elif ci_width < 3.0:
+                confidence = "MEDIUM"
+            else:
+                confidence = "LOW"
+
+            if effect > ci_width * 0.5 and effect > 0:
+                recommendation = "BUY"
+            elif effect < -ci_width * 0.5 and effect < 0:
+                recommendation = "SELL"
+            else:
+                recommendation = "AVOID"
+
+            results.append({
+                "player_name": row.get("name", row.get("player_name", str(idx))),
+                "role": candidate_role,
+                "current_team": row.get("team", row.get("squadra", "")),
+                "base_projection": round(base_proj, 2),
+                "causal_effect": round(effect, 3),
+                "ci_lower": round(ci_lo, 3),
+                "ci_upper": round(ci_hi, 3),
+                "synergy_bonus": round(synergy_bonus, 3),
+                "redundancy_penalty": round(redundancy_penalty, 3),
+                "recommendation": recommendation,
+                "confidence": confidence,
+            })
+
+        df = pd.DataFrame(results)
+        df = df.sort_values("causal_effect", ascending=False)
+        df = df.head(n_recommendations).reset_index(drop=True)
+
+        logger.info(
+            f"Ranked {len(results)} candidates; returning top {len(df)}"
+        )
+        return df
+
+    # ------------------------------------------------------------------
+    # Confounder analysis
+    # ------------------------------------------------------------------
+
+    def analyze_confounders(
+        self,
+        X: pd.DataFrame,
+        treatment: pd.DataFrame,
+        outcome: pd.Series,
+    ) -> List[Dict]:
+        """Identify roster features that confound the transfer effect.
+
+        A confounder correlates with both the treatment assignment and
+        the outcome — meaning it could bias a naive estimate.
+
+        Parameters
+        ----------
+        X : pd.DataFrame
+            Roster features.
+        treatment : pd.DataFrame
+            Treatment features.
+        outcome : pd.Series
+            Outcome values.
+
+        Returns
+        -------
+        list of dict
+            Each dict has keys: variable, corr_with_treatment, corr_with_outcome,
+            confounder (bool), partial_corr
+        """
+        X_num = X.select_dtypes(include=[np.number]).fillna(0)
+        T_num = np.asarray(treatment.select_dtypes(include=[np.number]).fillna(0), dtype=np.float64)
+        Y_num = np.asarray(outcome, dtype=np.float64).ravel()
+
+        # Determine treatment variable (first treatment column, or binary indicator)
+        if T_num.ndim == 2 and T_num.shape[1] > 0:
+            T_primary = T_num[:, 0]
+        else:
+            T_primary = T_num.ravel()
+
+        results = []
+
+        for col in X_num.columns:
+            x_vals = _safe_float_series(X_num[col])
+
+            # Correlation with treatment
+            corr_T = _safe_corr(x_vals, T_primary)
+
+            # Correlation with outcome
+            corr_Y = _safe_corr(x_vals, Y_num)
+
+            # Partial correlation: corr(X, Y | T)
+            partial = _partial_corr(x_vals, Y_num, T_primary)
+
+            # Confounder if correlates with both
+            is_conf = abs(corr_T) > 0.05 and abs(corr_Y) > 0.05
+
+            results.append({
+                "variable": col,
+                "corr_with_treatment": round(corr_T, 4),
+                "corr_with_outcome": round(corr_Y, 4),
+                "confounder": is_conf,
+                "partial_corr": round(partial, 4),
+            })
+
+        results.sort(key=lambda d: abs(d["partial_corr"]), reverse=True)
+
+        n_conf = sum(1 for r in results if r["confounder"])
+        logger.info(
+            f"Confounder analysis: {n_conf}/{len(results)} variables identified "
+            "as potential confounders"
+        )
+        return results
+
+    # ------------------------------------------------------------------
+    # Subgroup effects
+    # ------------------------------------------------------------------
+
+    def subgroup_effects(
+        self,
+        X: pd.DataFrame,
+        treatment: pd.DataFrame,
+        outcome: pd.Series,
+        subgroups: Dict[str, np.ndarray],
+    ) -> pd.DataFrame:
+        """Estimate treatment effect by subgroup.
+
+        Parameters
+        ----------
+        X : pd.DataFrame
+            Roster features.
+        treatment : pd.DataFrame
+            Treatment features.
+        outcome : pd.Series
+            Outcome values.
+        subgroups : dict
+            Mapping of subgroup name → boolean mask array.
+
+        Returns
+        -------
+        pd.DataFrame
+            Columns: subgroup, n, ate, ci_lower, ci_upper
+        """
+        X_num = X.select_dtypes(include=[np.number]).fillna(0)
+        T_num = treatment.select_dtypes(include=[np.number]).fillna(0)
+        Y_num = np.asarray(outcome, dtype=np.float64).ravel()
+
+        results = []
+        for name, mask in subgroups.items():
+            mask = np.asarray(mask, dtype=bool)
+            n = int(mask.sum())
+            if n < 5:
+                logger.warning(
+                    f"Subgroup '{name}' has only {n} samples; skipping"
+                )
+                continue
+
+            X_sub = X_num.loc[mask] if isinstance(X_num, pd.DataFrame) else X_num[mask]
+            T_sub = T_num.loc[mask] if isinstance(T_num, pd.DataFrame) else T_num[mask]
+            Y_sub = Y_num[mask]
+
+            if isinstance(X_sub, pd.DataFrame):
+                X_sub = X_sub.values
+            if isinstance(T_sub, pd.DataFrame):
+                T_sub = T_sub.values
+
+            # Compute ATE within subgroup using difference-in-means
+            if T_sub.ndim == 2 and T_sub.shape[1] > 0:
+                T_flat = T_sub[:, 0].ravel()
+            else:
+                T_flat = T_sub.ravel()
+
+            unique_t = np.unique(T_flat)
+            if set(unique_t).issubset({0.0, 1.0}) and len(unique_t) >= 2:
+                t1 = Y_sub[T_flat >= 0.5]
+                t0 = Y_sub[T_flat < 0.5]
+                if len(t1) > 0 and len(t0) > 0:
+                    ate = float(t1.mean() - t0.mean())
+                    se = np.sqrt(
+                        t1.var(ddof=1) / max(len(t1), 1) +
+                        t0.var(ddof=1) / max(len(t0), 1)
+                    ) if len(t1) > 1 and len(t0) > 1 else 1.0
+                else:
+                    ate = 0.0
+                    se = 1.0
+            else:
+                # Continuous treatment: use covariance
+                T_c = T_flat - T_flat.mean()
+                Y_c = Y_sub - Y_sub.mean()
+                denom = np.dot(T_c, T_c)
+                if denom > 1e-12:
+                    ate = float(np.dot(T_c, Y_c) / denom)
+                else:
+                    ate = 0.0
+                residual = Y_c - ate * T_c
+                se = (
+                    np.sqrt(np.dot(residual, residual) / (max(n - 2, 1)) / max(denom, 1e-10))
+                    if n > 2
+                    else 1.0
+                )
+
+            ci_lo = ate - 1.96 * se
+            ci_hi = ate + 1.96 * se
+
+            results.append({
+                "subgroup": name,
+                "n": n,
+                "ate": round(ate, 4),
+                "ci_lower": round(ci_lo, 4),
+                "ci_upper": round(ci_hi, 4),
+                "significant": abs(ate) > 1.96 * se,
+            })
+
+        df = pd.DataFrame(results)
+        df = df.sort_values("ate", ascending=False).reset_index(drop=True)
+        logger.info(
+            f"Subgroup effects computed: {len(results)} subgroups"
+        )
+        return df
+
+
+# ---------------------------------------------------------------------------
+# Helper functions
+# ---------------------------------------------------------------------------
+
+def _safe_corr(a: np.ndarray, b: np.ndarray) -> float:
+    """Compute correlation, returning 0.0 on degenerate input."""
+    a = np.asarray(a, dtype=np.float64).ravel()
+    b = np.asarray(b, dtype=np.float64).ravel()
+    if len(a) < 2 or len(b) < 2:
+        return 0.0
+    a_std = a.std()
+    b_std = b.std()
+    if a_std < 1e-12 or b_std < 1e-12:
+        return 0.0
+    return float(np.corrcoef(a, b)[0, 1])
+
+
+def _partial_corr(x: np.ndarray, y: np.ndarray, z: np.ndarray) -> float:
+    """Compute partial correlation corr(x, y | z)."""
+    x = np.asarray(x, dtype=np.float64).ravel()
+    y = np.asarray(y, dtype=np.float64).ravel()
+    z = np.asarray(z, dtype=np.float64).ravel()
+    r_xy = _safe_corr(x, y)
+    r_xz = _safe_corr(x, z)
+    r_yz = _safe_corr(y, z)
+    denom = np.sqrt((1 - r_xz ** 2) * (1 - r_yz ** 2))
+    if denom < 1e-12:
+        return 0.0
+    return float((r_xy - r_xz * r_yz) / denom)
+
+
+def _jackknife_se(values: np.ndarray) -> np.ndarray:
+    """Jackknife estimate of standard error."""
+    n = len(values)
+    if n <= 1:
+        return np.full_like(values, 1.0)
+    total = values.sum()
+    inf = (total - values) / (n - 1)
+    se_global = inf.std(ddof=1) * np.sqrt(n - 1)
+    return np.full_like(values, max(se_global, 1e-6))
+
+
+# ---------------------------------------------------------------------------
+# AuctionEffectAnalyzer
+# ---------------------------------------------------------------------------
+
+@dataclass
+class AuctionEffectAnalyzer:
+    """Integrates causal transfer effects with auction bidding strategy.
+
+    Adjusts fantasy-auction bids based on the predicted causal effect
+    of adding a player to the current roster composition.
+
+    Parameters
+    ----------
+    causal_model : TransferCausalModel
+        Fitted causal forest model.
+    current_roster_df : pd.DataFrame
+        Current team roster.
+    """
+
+    causal_model: TransferCausalModel
+    roster: pd.DataFrame
+
+    def recommend_bid_adjustment(
+        self,
+        player: pd.Series,
+        base_bid: float,
+        player_pool_df: pd.DataFrame,
+        max_total_budget: float = 500.0,
+    ) -> Dict:
+        """Recommend an adjusted bid based on causal effect on team composition.
+
+        Parameters
+        ----------
+        player : pd.Series
+            Player data for the auction target.
+        base_bid : float
+            Baseline bid (e.g., from projected points model).
+        player_pool_df : pd.DataFrame
+            Full pool of available players (for scarcity analysis).
+        max_total_budget : float
+            League budget cap.
+
+        Returns
+        -------
+        dict with keys: adjusted_bid, adjustment, adjustment_reason,
+            causal_effect, confidence
+        """
+        if not self.causal_model._fitted:
+            raise RuntimeError("Causal model must be fitted before use.")
+
+        result = self.causal_model.predict_individual_effect(
+            self.roster, player, player_pool_df=player_pool_df
+        )
+
+        effect = result["effect"]
+        ci_width = result["ci_upper"] - result["ci_lower"]
+        confidence = result["confidence"]
+        synergy = result["synergy_bonus"]
+        redundancy = result["redundancy_penalty"]
+
+        # Rule-based adjustment
+        reasons: List[str] = []
+
+        # Base adjustment proportional to effect, bounded
+        adjustment_pct = np.clip(effect / 10.0, -0.3, 0.3)
+
+        # Confidence modifier
+        if confidence == "LOW":
+            adjustment_pct *= 0.25
+            reasons.append("Low confidence — adjustment dampened")
+        elif confidence == "MEDIUM":
+            adjustment_pct *= 0.6
+
+        if synergy > 0.5:
+            adj_synergy = min(synergy / 10.0, 0.15)
+            adjustment_pct += adj_synergy
+            reasons.append(f"Positive synergy bonus (+{adj_synergy:.1%})")
+
+        if redundancy > 1.0:
+            adj_red = min(redundancy / 10.0, 0.2)
+            adjustment_pct -= adj_red
+            reasons.append(f"Redundancy penalty (-{adj_red:.1%})")
+
+        # Budget constraint
+        current_spent = float(self.roster.get("cost", pd.Series([0])).sum()) if "cost" in self.roster.columns else 0.0
+        budget_factor = 1.0 - (current_spent / max_total_budget)
+        if budget_factor < 0.3:
+            adjustment_pct -= 0.1
+            reasons.append("Budget nearly exhausted — downward adjustment")
+
+        adjusted_bid = base_bid * (1.0 + adjustment_pct)
+        adjusted_bid = max(adjusted_bid, 1.0)
+
+        if not reasons:
+            reasons.append(
+                f"Neutral adjustment — net effect {effect:+.2f} points/game"
+            )
+
+        return {
+            "adjusted_bid": round(adjusted_bid, 1),
+            "adjustment": f"{adjustment_pct:+.1%}",
+            "adjustment_reason": "; ".join(reasons),
+            "causal_effect": effect,
+            "confidence": confidence,
+        }