End-to-End Media Mix Optimization with Papilon: A Simulation-Driven Approach to MMM

End-to-End Media Mix Optimization with Papilon: A Rigorous Simulation-Driven Approach to Causal Marketing Measurement

Abstract

Marketing Mix Modeling (MMM) has experienced renewed relevance in the privacy-first digital ecosystem, yet traditional implementations suffer from fundamental limitations in causal identification, model specification, and strategic optimization. This article presents a comprehensive technical overview of Papilon, an open-source framework that unifies Bayesian hierarchical modeling, causal inference principles, agent-based simulation, and constrained optimization to enable end-to-end media mix measurement and allocation. We detail the mathematical foundations, computational architecture, and algorithmic implementations that differentiate Papilon from conventional MMM tools, with emphasis on identifiability conditions, posterior inference strategies, and multi-objective optimization under uncertainty.

1. Introduction: The Causal Inference Challenge in Marketing Measurement

1.1 The Fundamental Problem of Causal Identification

Marketing Mix Modeling aims to decompose aggregate business outcomes Y_t into contributions from multiple marketing inputs X_{t,k} while accounting for confounders Z_t (seasonality, macroeconomic factors, promotions). The central challenge is causal identification: distinguishing correlation from causation in observational time series data.

Formally, we seek to estimate the average treatment effect (ATE) of each marketing channel:

τ_k = E[Y_t(X_{t,k} = x') - Y_t(X_{t,k} = x)] | do(X_{t,k})

where Y_t(x) denotes the potential outcome under intervention do(X_{t,k} = x), following Pearl’s causal calculus. Traditional regression-based MMM conflates the conditional expectation E[Y_t | X_{t,k}] with the interventional distribution P(Y_t | do(X_{t,k})), violating fundamental identifiability conditions when unobserved confounding is present.

1.2 Privacy-First Measurement and the Return to Aggregation

The deprecation of third-party cookies and increasingly stringent privacy regulations (GDPR, CCPA, iOS ATT framework) have rendered user-level attribution increasingly untenable. Multi-touch attribution (MTA) models suffer from:

1. Sparse signal problems: Touch-point data becomes increasingly incomplete
2. Selection bias: Users who accept tracking differ systematically from those who opt out
3. Cross-device fragmentation: User journey reconstruction fails across devices
4. Counterfactual invalidity: Last-touch and linear heuristics provide no causal basis

MMM operates on aggregated time series data, sidestepping individual tracking while providing holistic, privacy-compliant measurement. However, aggregation introduces new challenges: temporal confounding, multicollinearity, and the loss of heterogeneity in treatment effects.

1.3 Papilon’s Methodological Innovation

Papilon addresses these challenges through five integrated innovations:

1. Bayesian hierarchical specification with informative priors encoding domain knowledge
2. Causal identification through integration of experimental data and instrumental variable techniques
3. Non-linear saturation and carryover modeling via flexible functional transformations
4. Agent-based simulation capturing system dynamics and channel interactions
5. Multi-objective optimization under budget constraints with uncertainty quantification

2. Mathematical Foundations of Papilon’s MMM Framework

2.1 The Core Additive Model Specification

Papilon implements a hierarchical generalized additive model (GAM) structure:

Y_t = β_0 + Σ_k f_k(X_{t,k}) + g(Z_t) + ε_t

where:

• Y_t is the target KPI (sales, conversions, revenue) at time t
• X_{t,k} represents raw media inputs (spend, impressions, GRPs) for channel k
• f_k(·) is a non-linear transformation function capturing saturation and carryover
• g(Z_t) captures control variables (trend, seasonality, promotions, macroeconomic factors)
• ε_t ~ N(0, σ²) is the error term

2.2 Adstock Transformation: Distributed Lag Structures

Media effects exhibit temporal dynamics: exposures today influence outcomes across future periods. The geometric adstock transformation models exponential decay:

X*_{t,k} = Σ_{l=0}^L α^l · X_{t-l,k}

where α ∈ [0,1] is the retention rate. This specification assumes that l periods after exposure, α^l fraction of the original effect persists.

Papilon extends this with flexible adstock forms:

2.2.1 Delayed Adstock (Lag + Decay)

X*_{t,k} = Σ_{l=θ}^{θ+L} α^{l-θ} · X_{t-l,k}

where θ represents the lag before effects begin (common in awareness campaigns).

2.2.2 Weibull Adstock (Non-Monotonic Decay)

w(l; λ, k) = (k/λ) · (l/λ)^{k-1} · exp(-(l/λ)^k)
X*_{t,k} = Σ_{l=0}^L w(l; λ, k) · X_{t-l,k}

The Weibull distribution allows flexible shapes: k < 1 gives immediate peak with long tail, k > 1 models delayed peak effects (e.g., brand campaigns building over time).

2.2.3 Carryover Normalization

To maintain interpretability, Papilon normalizes adstock weights:

X*_{t,k} = X_{t,k} + α · X*_{t-1,k}

This recursive formulation ensures that the sum of distributed effects equals total spend when integrated over time.

2.3 Saturation Modeling: Hill Function and Extensions

Marketing channels exhibit diminishing marginal returns. The Hill saturation function, adapted from pharmacokinetics, provides a flexible S-shaped curve:

f_k(x) = α_k · x^ν_k / (κ^ν_k + x^ν_k)

where:

• α_k: maximum achievable effect (asymptotic ceiling)
• κ: half-saturation point (spend level achieving 50% of maximum effect)
• ν: shape parameter controlling steepness (elasticity)

Key properties:

• When ν = 1: Michaelis-Menten kinetics (hyperbolic saturation)
• When ν > 1: Sigmoidal curve with inflection point
• When ν < 1: Convex saturation (rarely used)

2.3.1 Elasticity Derivation

The point elasticity at current spend level x₀:

ε(x₀) = (∂f/∂x) · (x₀/f(x₀))
     = ν · κ^ν / (κ^ν + x₀^ν) · (1 - f(x₀)/α)

This allows comparing channel efficiency: higher elasticity indicates room for incremental investment.

2.3.2 Marginal ROI

The instantaneous return on ad spend (iROAS):

iROAS_k(x) = ∂f_k(x)/∂x · (conversion_value/cost_per_unit)

Papilon uses this derivative to identify optimal budget allocation where marginal returns equalize across channels.

2.4 Bayesian Hierarchical Specification

Rather than point estimates, Papilon fits full posterior distributions over parameters using Hamiltonian Monte Carlo (HMC) via PyMC.

2.4.1 Prior Specification

Saturation parameters (informed by marketing theory):

α_k ~ Normal(μ_sales, σ_sales)  # Scaled to outcome range
κ_k ~ Gamma(shape=2, rate=1/mean_spend)  # Centered at typical spend
ν_k ~ Normal(1, 0.5)  # Weakly informative, allows convex or sigmoidal

Adstock parameters:

α_k ~ Beta(α_prior, β_prior)  # Reflects channel-specific carryover expectations
θ_k ~ Geometric(p_lag)  # Lag parameter for delayed effects

Channel coefficients:

β_k ~ Normal(0, σ_β)  # Regularization via prior variance

2.4.2 Hierarchical Structure for Multiple Geographies/Products

When data spans multiple markets or product lines:

β_{k,j} ~ Normal(μ_k, σ_k)  # Channel k, region/product j
μ_k ~ Normal(0, σ_μ)  # Population-level mean
σ_k ~ HalfCauchy(1)  # Heterogeneity scale

This partial pooling shrinks individual estimates toward the population mean, stabilizing inference when regional data is sparse.

2.4.3 Posterior Inference via HMC

Papilon uses the No-U-Turn Sampler (NUTS), an adaptive HMC variant:

1. Transform parameters to unconstrained space
2. Compute gradient of log-posterior using automatic differentiation
3. Simulate Hamiltonian dynamics with adaptive step size
4. Accept/reject proposals based on energy error
5. Dynamically tune trajectory length to eliminate manual tuning

Advantages over maximum likelihood:

• Full uncertainty quantification (credible intervals)
• Robust to multicollinearity (regularization via priors)
• Natural model comparison via Bayesian model selection (LOO-CV, WAIC)

2.5 Causal Identification Strategies

2.5.1 Integration of Experimental Data

When randomized experiments (geo-lift tests, incrementality studies) are available, Papilon incorporates them as calibration constraints:

τ_k^{exp} ~ Normal(f_k(X_{treat}) - f_k(X_{control}), σ_exp²)

The model’s predicted lift must align with experimentally measured effects, anchoring correlational patterns to causal ground truth.

2.5.2 Instrumental Variable Approach

For channels susceptible to endogeneity (e.g., search spending responding to organic demand spikes), Papilon supports IV estimation:

X_{t,k} = γ_0 + γ_1 · Z_t + u_t  # First stage
Y_t = β_0 + β_1 · X̂_{t,k} + ε_t  # Second stage

where Z_t is an instrument (e.g., pre-planned budget allocations, supply-side shocks) uncorrelated with demand shocks but correlated with realized spend.

2.5.3 Difference-in-Differences for Policy Changes

When media strategy changes sharply (e.g., entering a new channel), DiD estimators isolate the causal effect:

Y_{t,treated} - Y_{t,control} = β_DD + τ · Post_t · Treated + ε_t

Papilon’s simulation engine can replicate synthetic control units when parallel trends assumptions hold.

3. Agent-Based Simulation Architecture

3.1 Motivation: Beyond Linear Response Models

Traditional MMM assumes channel independence and static response functions. Reality exhibits:

1. Cross-channel synergies: TV increases branded search volume
2. Competitive dynamics: Rival spending influences market share
3. Consumer heterogeneity: Different segments respond differently
4. Feedback loops: Sales success enables higher spending

Agent-based models (ABMs) capture these emergent phenomena by simulating micro-level behaviors that aggregate to macro-level outcomes.

3.2 Agent Architecture in Papilon

Each agent represents a consumer or micro-segment with:

• State variables: awareness level, consideration set, purchase propensity
• Exposure accumulation: Weighted sum of media contacts across channels
• Response function: Probabilistic conversion based on accumulated exposures
• Memory: Temporal decay of exposures (linking to adstock)

3.2.1 Agent State Transition Dynamics

At each time step t:

Awareness_{i,t+1} = Awareness_{i,t} · (1 - δ) + Σ_k w_k · Contact_{i,k,t}

where δ is the natural decay rate, and Contact_{i,k,t} ~ Bernoulli(reach_k,t) indicates whether agent i was exposed to channel k.

Consideration probability:

P(Consider_i,t) = 1 / (1 + exp(-β_A · Awareness_{i,t}))

Purchase probability (conditional on consideration):

P(Purchase_{i,t} | Consider_i,t) = 1 / (1 + exp(-β_P · Consideration_{i,t} - γ · Recency_i,t))

3.2.2 Channel Interaction Modeling

Synergy effects: Search response amplifies after display exposure

Contact_{i,search,t}' = Contact_{i,search,t} · (1 + φ · I(Contact_{i,display,t-lag}))

Competitive interference: High competitive pressure reduces effectiveness

Effectiveness_{k,t} = baseline_k · exp(-λ · Competitive_spend_{k,t})

3.3 Simulation Protocol

1. Initialize population: N agents with heterogeneous characteristics drawn from empirical distributions
2. Time-step iteration: For t = 1 to T: a. Allocate media impressions based on targeting rules b. Update agent states based on exposures c. Sample conversions from purchase probabilities d. Aggregate to observable outcomes Y_t
3. Parameter estimation: Calibrate agent-level parameters to match observed aggregate data using simulated method of moments or approximate Bayesian computation

3.4 Advantages of Simulation-Based MMM

• Counterfactual generation: Simulate outcomes under arbitrary budget allocations
• Non-stationarity: Parameters can evolve over time (e.g., creative fatigue)
• Strategic testing: Evaluate new channel entry before real-world deployment
• Attribution decomposition: Track which agent conversions resulted from which channel sequences

4. Papilon Implementation: Modular Software Architecture

4.1 Core Modules

4.1.1 Data Processing Pipeline

Class: PapilonDataTransformer

python

class PapilonDataTransformer:
    def __init__(self, adstock_params: Dict, saturation_params: Dict):
        self.adstock = adstock_params
        self.saturation = saturation_params
    
    def apply_adstock(self, X: pd.DataFrame, channel: str) -> pd.Series:
        """
        Apply geometric or Weibull adstock transformation
        
        Parameters:
        -----------
        X : DataFrame with columns [date, channel_spend]
        channel : str, channel identifier
        
        Returns:
        --------
        Transformed series with carryover effects
        """
        alpha = self.adstock[channel]['retention']
        theta = self.adstock[channel].get('lag', 0)
        
        X_adstock = np.zeros(len(X))
        for t in range(len(X)):
            for l in range(theta, min(t+1, theta+self.adstock['max_lag'])):
                X_adstock[t] += alpha**(l-theta) * X.iloc[t-l][channel]
        
        return pd.Series(X_adstock, index=X.index)
    
    def apply_hill_saturation(self, X_adstock: pd.Series, channel: str) -> pd.Series:
        """
        Apply Hill saturation transformation
        
        Parameters:
        -----------
        X_adstock : Series, adstocked media input
        channel : str, channel identifier
        
        Returns:
        --------
        Saturated effect values
        """
        alpha = self.saturation[channel]['alpha']
        kappa = self.saturation[channel]['kappa']
        nu = self.saturation[channel]['nu']
        
        X_effect = alpha * (X_adstock**nu) / (kappa**nu + X_adstock**nu)
        return X_effect

4.1.2 Bayesian Model Specification

Class: MarketingMixModel

python

import pymc as pm
import pytensor.tensor as pt

class MarketingMixModel:
    def __init__(self, data: pd.DataFrame, target: str, channels: List[str],
                 controls: List[str] = None):
        self.data = data
        self.target = target
        self.channels = channels
        self.controls = controls or []
        self.model = None
        self.trace = None
    
    def build_model(self, prior_configs: Dict):
        """
        Construct PyMC model with hierarchical structure
        """
        with pm.Model() as self.model:
            # Priors for saturation parameters
            alpha = {}
            kappa = {}
            nu = {}
            
            for ch in self.channels:
                alpha[ch] = pm.Normal(f'alpha_{ch}', 
                                     mu=prior_configs[ch]['alpha_mu'],
                                     sigma=prior_configs[ch]['alpha_sigma'])
                
                kappa[ch] = pm.Gamma(f'kappa_{ch}',
                                    alpha=prior_configs[ch]['kappa_shape'],
                                    beta=prior_configs[ch]['kappa_rate'])
                
                nu[ch] = pm.TruncatedNormal(f'nu_{ch}', mu=1.0, sigma=0.5,
                                           lower=0.1, upper=3.0)
            
            # Adstock parameters
            retention = {}
            for ch in self.channels:
                retention[ch] = pm.Beta(f'retention_{ch}', alpha=3, beta=3)
            
            # Apply transformations
            X_effects = []
            for ch in self.channels:
                # Adstock
                X_adstock = self._apply_adstock_pytensor(
                    self.data[ch].values, retention[ch]
                )
                # Hill saturation
                X_effect = alpha[ch] * pt.power(X_adstock, nu[ch]) / \
                          (pt.power(kappa[ch], nu[ch]) + pt.power(X_adstock, nu[ch]))
                X_effects.append(X_effect)
            
            # Control variables
            control_effects = []
            for ctrl in self.controls:
                beta_ctrl = pm.Normal(f'beta_{ctrl}', mu=0, sigma=1)
                control_effects.append(beta_ctrl * self.data[ctrl].values)
            
            # Intercept and error
            intercept = pm.Normal('intercept', mu=self.data[self.target].mean(),
                                 sigma=self.data[self.target].std())
            sigma = pm.HalfNormal('sigma', sigma=self.data[self.target].std())
            
            # Likelihood
            mu = intercept + pt.sum(X_effects, axis=0)
            if control_effects:
                mu += pt.sum(control_effects, axis=0)
            
            y_obs = pm.Normal('y_obs', mu=mu, sigma=sigma,
                            observed=self.data[self.target].values)
    
    def fit(self, draws: int = 2000, tune: int = 1000, chains: int = 4):
        """
        Sample from posterior using NUTS
        """
        with self.model:
            self.trace = pm.sample(draws=draws, tune=tune, chains=chains,
                                  target_accept=0.95, return_inferencedata=True)
        return self.trace
    
    def posterior_predictive_check(self):
        """
        Generate posterior predictive samples for model validation
        """
        with self.model:
            ppc = pm.sample_posterior_predictive(self.trace)
        return ppc

4.1.3 Agent-Based Simulation Engine

Class: MarketingSimulator

python

class MarketingSimulator:
    def __init__(self, n_agents: int, n_periods: int, channels: List[str]):
        self.n_agents = n_agents
        self.n_periods = n_periods
        self.channels = channels
        self.agents = self._initialize_agents()
    
    def _initialize_agents(self) -> pd.DataFrame:
        """
        Create agent population with heterogeneous characteristics
        """
        agents = pd.DataFrame({
            'agent_id': range(self.n_agents),
            'awareness': np.random.beta(2, 5, self.n_agents),
            'propensity': np.random.beta(2, 8, self.n_agents),
            'segment': np.random.choice(['high_value', 'medium_value', 'low_value'],
                                       self.n_agents, p=[0.2, 0.5, 0.3])
        })
        return agents
    
    def simulate_period(self, media_plan: Dict, t: int) -> pd.DataFrame:
        """
        Simulate one time period of agent exposures and conversions
        
        Parameters:
        -----------
        media_plan : Dict mapping channels to spend levels
        t : int, time period index
        
        Returns:
        --------
        DataFrame with agent-level outcomes
        """
        exposures = self._allocate_exposures(media_plan)
        
        # Update awareness
        for ch in self.channels:
            exposed_agents = exposures[exposures[ch] > 0].index
            self.agents.loc[exposed_agents, 'awareness'] += \
                0.05 * exposures.loc[exposed_agents, ch]
        
        # Apply decay
        self.agents['awareness'] *= 0.95
        
        # Compute conversion probabilities
        conversion_prob = 1 / (1 + np.exp(-5 * (self.agents['awareness'] - 0.5)))
        conversion_prob *= self.agents['propensity']
        
        # Sample conversions
        conversions = np.random.binomial(1, conversion_prob)
        
        return pd.DataFrame({
            'agent_id': self.agents['agent_id'],
            'period': t,
            'awareness': self.agents['awareness'],
            'converted': conversions
        })
    
    def run_simulation(self, media_plans: List[Dict]) -> pd.DataFrame:
        """
        Run full simulation across all time periods
        """
        results = []
        for t, plan in enumerate(media_plans):
            period_results = self.simulate_period(plan, t)
            results.append(period_results)
        
        return pd.concat(results, ignore_index=True)
    
    def _allocate_exposures(self, media_plan: Dict) -> pd.DataFrame:
        """
        Translate spend into agent-level exposures using reach curves
        """
        exposures = pd.DataFrame(index=self.agents.index)
        
        for ch, spend in media_plan.items():
            # Reach curve: more spend → higher reach, diminishing returns
            reach = 1 - np.exp(-spend / 10000)
            n_exposed = int(reach * self.n_agents)
            
            # Frequency distribution
            exposed_agents = np.random.choice(self.agents.index, n_exposed, replace=False)
            frequency = np.random.poisson(spend / 5000, n_exposed)
            
            exposures[ch] = 0
            exposures.loc[exposed_agents, ch] = frequency
        
        return exposures

4.1.4 Budget Optimization Module

Class: BudgetOptimizer

python

from scipy.optimize import minimize, LinearConstraint

class BudgetOptimizer:
    def __init__(self, mmm_model: MarketingMixModel, simulator: MarketingSimulator = None):
        self.model = mmm_model
        self.simulator = simulator
    
    def optimize_allocation(self, total_budget: float, 
                          channel_bounds: Dict = None,
                          n_simulations: int = 100) -> Dict:
        """
        Find optimal budget allocation using expected value maximization
        
        Parameters:
        -----------
        total_budget : float, total marketing budget
        channel_bounds : Dict with 'min' and 'max' for each channel
        n_simulations : int, number of posterior samples to average over
        
        Returns:
        --------
        Dict with optimal spend per channel
        """
        n_channels = len(self.model.channels)
        
        # Objective: maximize expected conversions
        def objective(x):
            # x is allocation vector (one value per channel)
            expected_outcome = 0
            
            # Average over posterior uncertainty
            for _ in range(n_simulations):
                sample_params = self._sample_posterior_params()
                outcome = self._predict_outcome(x, sample_params)
                expected_outcome += outcome
            
            return -expected_outcome / n_simulations  # Minimize negative
        
        # Constraints
        constraints = [
            {'type': 'eq', 'fun': lambda x: np.sum(x) - total_budget}  # Budget constraint
        ]
        
        # Bounds
        if channel_bounds:
            bounds = [(channel_bounds[ch]['min'], channel_bounds[ch]['max']) 
                     for ch in self.model.channels]
        else:
            bounds = [(0, total_budget) for _ in self.model.channels]
        
        # Initial guess: equal allocation
        x0 = np.ones(n_channels) * (total_budget / n_channels)
        
        # Optimize
        result = minimize(objective, x0, method='SLSQP', 
                         bounds=bounds, constraints=constraints)
        
        # Format results
        optimal_allocation = {ch: result.x[i] 
                            for i, ch in enumerate(self.model.channels)}
        
        return {
            'allocation': optimal_allocation,
            'expected_outcome': -result.fun,
            'optimization_success': result.success
        }
    
    def _sample_posterior_params(self) -> Dict:
        """
        Draw random sample from posterior distribution
        """
        sample_idx = np.random.randint(0, len(self.model.trace.posterior.draw))
        params = {}
        
        for ch in self.model.channels:
            params[ch] = {
                'alpha': self.model.trace.posterior[f'alpha_{ch}'][0, sample_idx].values,
                'kappa': self.model.trace.posterior[f'kappa_{ch}'][0, sample_idx].values,
                'nu': self.model.trace.posterior[f'nu_{ch}'][0, sample_idx].values,
                'retention': self.model.trace.posterior[f'retention_{ch}'][0, sample_idx].values
            }
        
        return params
    
    def _predict_outcome(self, allocation: np.ndarray, params: Dict) -> float:
        """
        Predict outcome for given allocation and parameter values
        """
        total_effect = self.model.trace.posterior['intercept'][0, 0].values
        
        for i, ch in enumerate(self.model.channels):
            # Apply adstock
            x_adstock = allocation[i] * params[ch]['retention'] / \
                       (1 - params[ch]['retention'])
            
            # Apply saturation
            effect = params[ch]['alpha'] * (x_adstock**params[ch]['nu']) / \
                    (params[ch]['kappa']**params[ch]['nu'] + x_adstock**params[ch]['nu'])
            
            total_effect += effect
        
        return total_effect
    
    def marginal_roi_curve(self, channel: str, spend_range: np.ndarray) -> pd.DataFrame:
        """
        Compute marginal ROI curve across spending range
        
        Returns DataFrame with columns: spend, roi, roi_lower, roi_upper
        """
        rois = []
        
        for spend in spend_range:
            # Compute derivative numerically
            delta = spend * 0.01
            outcome_plus = self._expected_outcome_single_channel(channel, spend + delta)
            outcome_minus = self._expected_outcome_single_channel(channel, spend - delta)
            
            marginal_effect = (outcome_plus - outcome_minus) / (2 * delta)
            roi = marginal_effect  # Assumes conversion value is standardized
            
            rois.append(roi)
        
        return pd.DataFrame({
            'spend': spend_range,
            'marginal_roi': rois
        })

4.2 Integration Workflow

python

# End-to-end pipeline
def run_papilon_analysis(data: pd.DataFrame, channels: List[str],
                         target: str, total_budget: float):
    """
    Complete MMM analysis pipeline
    """
    # 1. Data transformation
    transformer = PapilonDataTransformer(
        adstock_params={ch: {'retention': 0.7, 'lag': 0} for ch in channels},
        saturation_params={ch: {'alpha': 1.0, 'kappa': 10000, 'nu': 1.5} 
                          for ch in channels}
    )
    
    for ch in channels:
        data[f'{ch}_adstock'] = transformer.apply_adstock(data, ch)
        data[f'{ch}_effect'] = transformer.apply_hill_saturation(data[f'{ch}_adstock'], ch)
    
    # 2. Bayesian model fitting
    mmm = MarketingMixModel(data=data, target=target, channels=channels,
                           controls=['trend', 'seasonality'])
    
    prior_configs = {ch: {'alpha_mu': data[target].mean() / len(channels),
                         'alpha_sigma': data[target].std(),
                         'kappa_shape': 2.0,
                         'kappa_rate': 1 / data[ch].mean()}
                    for ch in channels}
    
    mmm.build_model(prior_configs)
    trace = mmm.fit(draws=3000, tune=1500, chains=4)
    
    # 3. Model diagnostics
    ppc = mmm.posterior_predictive_check()
    
    # 4. Budget optimization
    optimizer = BudgetOptimizer(mmm)
    optimal_plan = optimizer.optimize_allocation(
        total_budget=total_budget,
        channel_bounds={ch: {'min': 0, 'max': total_budget * 0.5} for ch in channels},
        n_simulations=200
    )
    
    # 5. Scenario simulation
    simulator = MarketingSimulator(n_agents=10000, n_periods=90, channels=channels)
    
    # Simulate current plan
    current_plan = [data[channels].iloc[t].to_dict() for t in range(len(data))]
    current_results = simulator.run_simulation(current_plan)
    
    # Simulate optimal plan
    optimal_plan_repeated = [optimal_plan['allocation']] * 90
    optimal_results = simulator.run_simulation(optimal_plan_repeated)
    
    # 6. Results package
    return {
        'model': mmm,
        'trace': trace,
        'optimal_allocation': optimal_plan,
        'current_simulation': current_results,
        'optimal_simulation': optimal_results,
        'lift': (optimal_results['converted'].sum() / 
                current_results['converted'].sum() - 1)
    }

5. Advanced Topics and Extensions

5.1 Multi-Objective Optimization Under Uncertainty

Real-world marketing involves competing objectives: maximize revenue, minimize CAC, maintain brand health. Papilon supports Pareto-optimal frontier estimation:

python

def pareto_optimization(mmm: MarketingMixModel, objectives: List[callable],
                       budget: float) -> pd.DataFrame:
    """
    Find Pareto frontier for multi-objective allocation
    
    Parameters:
    -----------
    objectives : List of functions mapping allocation → objective value
    
    Returns:
    --------
    DataFrame with non-dominated solutions
    """
    from pymoo.algorithms.moo.nsga2 import NSGA2
    from pymoo.core.problem import Problem
    
    class MarketingProblem(Problem):
        def __init__(self, n_channels, budget, objectives):
            super().__init__(n_var=n_channels,
                           n_obj=len(objectives),
                           n_constr=1,
                           xl=np.zeros(n_channels),
                           xu=np.ones(n_channels) * budget)
            self.objectives = objectives
            self.budget = budget
        
        def _evaluate(self, X, out, *args, **kwargs):
            # Objectives (to minimize; negate if maximizing)
            obj_values = np.array([[obj(x) for obj in self.objectives] for x in X])
            out["F"] = -obj_values  # NSGA2 minimizes
            
            # Budget constraint
            out["G"] = np.array([np.sum(x) - self.budget for x in X])
    
    problem = MarketingProblem(len(mmm.channels), budget, objectives)
    algorithm = NSGA2(pop_size=100)
    
    from pymoo.optimize import minimize
    res = minimize(problem, algorithm, ('n_gen', 200), verbose=False)
    
    return pd.DataFrame(res.X, columns=mmm.channels)

5.2 Incorporating Competitive Intelligence

Competitor spending influences market share. Papilon models this through interaction terms:

python

# In model specification
competitor_spend = pm.Data('competitor_spend', data['competitor_spend_total'])

# Competitive pressure effect
competitive_modifier = pm.Deterministic(
    'competitive_modifier',
    pt.exp(-lambda_compete * competitor_spend)
)

# Adjust own-channel effectiveness
for ch in channels:
    X_effect_adjusted = X_effect[ch] * competitive_modifier

5.3 Time-Varying Parameters

Marketing effectiveness evolves due to creative fatigue, seasonality, market maturation. State-space formulations allow dynamic coefficients:

python

# Random walk priors for time-varying alpha
alpha_t = {}
for ch in channels:
    alpha_t[ch] = pm.GaussianRandomWalk(
        f'alpha_{ch}_t',
        mu=0, sigma=0.1, shape=n_periods
    )
    
    # Cumulative effect
    alpha_cumulative = pm.Deterministic(
        f'alpha_{ch}_cumulative',
        alpha_base[ch] + pt.cumsum(alpha_t[ch])
    )

5.4 Causal Mediation Analysis

Decompose total effect into direct and indirect pathways:

Total Effect = Direct Effect + Indirect Effect (via mediator)

Example: TV → Brand Search → Conversions

python

# Mediation model
with pm.Model() as mediation_model:
    # TV → Search (mediator equation)
    beta_tv_search = pm.Normal('beta_tv_search', mu=0, sigma=1)
    search_predicted = intercept_search + beta_tv_search * tv_spend
    
    # TV + Search → Conversions (outcome equation)
    beta_tv_direct = pm.Normal('beta_tv_direct', mu=0, sigma=1)
    beta_search = pm.Normal('beta_search', mu=0, sigma=1)
    
    conversions = intercept + beta_tv_direct * tv_spend + beta_search * search_predicted
    
    # Causal effects
    direct_effect = pm.Deterministic('direct_effect', beta_tv_direct)
    indirect_effect = pm.Deterministic('indirect_effect', beta_tv_search * beta_search)
    total_effect = pm.Deterministic('total_effect', direct_effect + indirect_effect

5.5 Robustness Checks and Sensitivity Analysis

Papilon provides tools for assessing model sensitivity to assumptions:

python

def sensitivity_analysis(mmm: MarketingMixModel, 
                        param: str, 
                        param_range: np.ndarray) -> pd.DataFrame:
    """
    Vary a prior parameter and observe posterior changes
    
    Useful for understanding influence of prior specification
    """
    results = []
    
    for value in param_range:
        # Refit model with modified prior
        modified_config = mmm.prior_configs.copy()
        modified_config[param] = value
        
        mmm_modified = MarketingMixModel(mmm.data, mmm.target, mmm.channels)
        mmm_modified.build_model(modified_config)
        trace_modified = mmm_modified.fit(draws=1000, tune=500)
        
        # Extract key posterior summaries
        summary = az.summary(trace_modified)
        results.append({
            'param_value': value,
            'posterior_mean': summary.loc['target_param', 'mean'],
            'posterior_std': summary.loc['target_param', 'sd']
        })
    
    return pd.DataFrame(results)

6. Model Validation and Diagnostics

6.1 Posterior Predictive Checks

Visual and statistical validation that model captures data-generating process:

python

import arviz as az

# Generate posterior predictive samples
ppc = mmm.posterior_predictive_check()

# Visual check
az.plot_ppc(ppc, num_pp_samples=100)

# Statistical checks
from scipy import stats

# Test for distributional match
observed = mmm.data[mmm.target].values
predicted = ppc.posterior_predictive['y_obs'].values.flatten()

ks_stat, p_value = stats.ks_2samp(observed, predicted)
print(f"KS test: statistic={ks_stat:.4f}, p-value={p_value:.4f}")

6.2 Leave-Future-Out Cross-Validation

Time series analog of leave-one-out CV:

python

def leave_future_out_cv(data: pd.DataFrame, mmm: MarketingMixModel,
                       n_folds: int = 5) -> pd.DataFrame:
    """
    Expanding window cross-validation for time series
    """
    fold_size = len(data) // n_folds
    results = []
    
    for i in range(1, n_folds):
        train_data = data.iloc[:fold_size * i]
        test_data = data.iloc[fold_size * i:fold_size * (i+1)]
        
        # Fit on train
        mmm_cv = MarketingMixModel(train_data, mmm.target, mmm.channels)
        mmm_cv.build_model(mmm.prior_configs)
        mmm_cv.fit(draws=1000, tune=500)
        
        # Predict on test
        predictions = mmm_cv.predict(test_data)
        actuals = test_data[mmm.target]
        
        # Compute metrics
        mae = np.mean(np.abs(predictions - actuals))
        mape = np.mean(np.abs((actuals - predictions) / actuals)) * 100
        
        results.append({
            'fold': i,
            'train_size': len(train_data),
            'test_size': len(test_data),
            'mae': mae,
            'mape': mape
        })
    
    return pd.DataFrame(results

6.3 Decomposition Validation Against Holdout Experiments

When geo-lift tests are available, validate model attributions:

python

def validate_against_experiment(mmm: MarketingMixModel,
                               experiment_results: pd.DataFrame) -> pd.DataFrame:
    """
    Compare model-based lift estimates to experimental ground truth
    
    Parameters:
    -----------
    experiment_results : DataFrame with columns [channel, treatment_spend, 
                        control_spend, measured_lift, lift_std_error]
    
    Returns:
    --------
    Validation summary showing agreement between model and experiment
    """
    validation = []
    
    for _, exp in experiment_results.iterrows():
        channel = exp['channel']
        
        # Model-predicted lift
        model_lift_treatment = mmm.predict_outcome({channel: exp['treatment_spend']})
        model_lift_control = mmm.predict_outcome({channel: exp['control_spend']})
        model_lift = model_lift_treatment - model_lift_control
        
        # Experimental lift
        exp_lift = exp['measured_lift']
        exp_se = exp['lift_std_error']
        
        # Compute agreement
        z_score = (model_lift - exp_lift) / exp_se
        agreement = np.abs(z_score) < 1.96  # Within 95% CI
        
        validation.append({
            'channel': channel,
            'model_lift': model_lift,
            'experimental_lift': exp_lift,
            'z_score': z_score,
            'agreement': agreement
        })
    
    return pd.DataFrame(validation)

7. Comparative Analysis: Papilon vs. Traditional MMM Tools

7.1 Methodological Comparison

DimensionTraditional MMMPapilonEstimationOrdinary Least SquaresBayesian Hierarchical (HMC)UncertaintyFrequentist confidence intervalsFull posterior distributionsSaturationLog or power transformsFlexible Hill function with priorsCarryoverFixed geometric lagGeometric, Weibull, or customInteractionsManual multiplicative termsAgent-based emergent dynamicsCausalityImplicit correlationExplicit causal integrationOptimizationCalculus of variationsStochastic simulation + constrained optimizationForecastingLinear extrapolationAgent-based scenario simulation

7.2 Computational Complexity

Traditional OLS MMM: O(n × p²) for n observations, p predictors

Papilon Bayesian MMM: O(n × p² × k × m) where k = HMC iterations, m = chains

• Typically 10–100× slower than OLS
• Parallelizes well across chains
• Worth it for: uncertainty quantification, regularization, complex hierarchies

Papilon ABM: O(n_agents × n_periods × n_channels)

• Expensive for large populations (>100K agents)
• Enables counterfactuals and dynamics unavailable in regression
• Use judiciously: calibrate on small representative sample, scale up for key analyses

7.3 Use Case Suitability

Traditional MMM best for:

• Rapid exploratory analysis
• Stable, mature markets with linear dynamics
• Limited computational resources
• Teams unfamiliar with Bayesian methods

Papilon best for:

• High-stakes budget allocation decisions requiring uncertainty quantification
• Markets with complex channel interactions and feedback loops
• Integration with experimental data for causal calibration
• Organizations with technical capacity for Bayesian/simulation methods
• Scenario planning and strategic foresight (new channel evaluation, competitive response)

8. Practical Implementation Considerations

8.1 Data Requirements

Minimum viable dataset:

• 52 weeks (1 year) of weekly data
• 5+ marketing channels
• Consistent KPI tracking (sales, conversions, revenue)
• Control variables (trend, seasonality, promotions, macroeconomic indices)

Ideal dataset:

• 2–3 years of weekly or daily data
• 10+ channels with varying spend levels
• Multiple geographies or product lines for hierarchical structure
• Experimental holdout periods or geo-lift tests for calibration
• Competitive intelligence data

8.2 Computational Infrastructure

For Bayesian fitting:

• Multi-core CPU (8+ cores recommended for parallelization)
• 16+ GB RAM for large models
• GPU acceleration (optional, via PyMC/JAX backend) provides 5–10× speedup

For agent-based simulation:

• High-memory nodes (32+ GB) for populations >50K agents
• Consider distributed computing (Dask, Ray) for large-scale simulations
• Estimate: 1M agent-periods ≈ 1–5 minutes on modern hardware

8.3 Deployment Architecture

[Data Lake/Warehouse]
         ↓
[ETL Pipeline] → Papilon Data Transformer
         ↓
[Bayesian Model Training] → Posterior Samples (stored)
         ↓
[Optimization Service] ← Business constraints
         ↓
[Scenario Simulator] → What-if analyses
         ↓
[Reporting Dashboard] → Stakeholder outputs

Suggested tech stack:

• Data: BigQuery, Snowflake, or Redshift
• Orchestration: Airflow or Prefect
• Modeling: Python (PyMC, NumPy, Pandas)
• API: FastAPI for serving optimization endpoints
• Visualization: Tableau, Looker, or Plotly Dash

8.4 Organizational Integration

Key stakeholders:

1. Marketing analysts: Primary users, must understand Bayesian outputs
2. Media buyers: Implement optimized allocations
3. Finance: Align budget constraints with recommendations
4. Data science: Maintain models, validate assumptions
5. Executives: Interpret strategic implications

Change management:

• Start with pilot on subset of channels
• Validate against existing attribution systems (MTA, last-touch)
• Run parallel for 2–3 quarters before full transition
• Invest in training: Bayesian thinking, interpreting credible intervals, simulation logic

9. Research Frontiers and Future Directions

9.1 Real-Time Adaptive MMM

Current MMM operates on retrospective data (weekly/monthly). Next generation:

• Streaming Bayesian updates: Incorporate daily data via particle filters or variational inference
• Online learning: Update posteriors as new data arrives without full refitting
• Adaptive allocation: Algorithms that shift spend mid-flight based on early signals

9.2 Causal Discovery and Structure Learning

Rather than assuming model structure, learn it from data:

• Bayesian networks: Discover causal graphs among channels
• Granger causality: Test temporal precedence relationships
• Structural equation modeling: Jointly estimate mediation pathways

9.3 Heterogeneous Treatment Effects

Not all customers respond equally. Model individual-level or segment-level variation:

• Conditional average treatment effects (CATE): Estimate lift by customer characteristics
• Causal forests: Machine learning for heterogeneous causal inference
• Integration with CRM data: Link MMM to customer lifetime value models

9.4 Multi-Armed Bandit Integration

Bridge MMM (retrospective analysis) with MAB (prospective optimization):

• Use MMM posteriors as priors for bandit algorithms
• Thompson sampling with MMM-informed beliefs
• Contextual bandits for dynamic budget allocation

9.5 Deep Learning Enhancements

• Neural MMM: Replace Hill/adstock with learned transformations (attention mechanisms, LSTM encoders)
• Physics-informed neural networks: Embed conservation laws (e.g., budget constraints) in loss functions
• Variational autoencoders: Learn latent representations of marketing mix states

10. Conclusion

Papilon represents a paradigm shift in marketing measurement: from descriptive correlation to prescriptive causal optimization. By unifying Bayesian hierarchical modeling, agent-based simulation, and constrained optimization, it enables marketing teams to move beyond “what happened” to “what will happen if we act.”

The framework’s technical sophistication — HMC sampling, flexible saturation functions, emergent system dynamics — demands investment in data infrastructure and analytical talent. However, for organizations operating at scale in competitive markets, the payoff is substantial: more efficient budget allocation, robust uncertainty quantification, and the ability to simulate strategic scenarios before committing real dollars.

As privacy regulations tighten and user-level tracking becomes infeasible, aggregate causal inference will become the cornerstone of marketing science. Papilon provides the methodological and computational foundation for this transition, positioning early adopters to extract maximum value from their marketing investments in an increasingly complex ecosystem.

Appendix A: Mathematical Notation Summary

SymbolDefinitionY_tOutcome (sales, conversions) at time tX_{t,k}Media input for channel k at time tf_k(·)Transformation function (adstock + saturation) for channel kα_kMaximum effect parameter (Hill function)κ_kHalf-saturation point (Hill function)ν_kShape parameter (Hill function)αRetention rate (geometric adstock)θLag parameter (delayed adstock)β_kLinear coefficient for channel kσ²Error varianceτ_kAverage treatment effect of channel kε(x)Elasticity at spend level xiROAS_k(x)Instantaneous return on ad spend

Appendix B: Key References and Further Reading

Foundational Papers

1. Bayesian MMM: Jin, Y., et al. (2017). “Bayesian Methods for Media Mix Modeling with Carryover and Shape Effects.” Google Research.
2. Causal Inference: Pearl, J. (2009). Causality: Models, Reasoning, and Inference. Cambridge University Press.
3. Agent-Based Marketing Models: Rand, W., & Rust, R. T. (2011). “Agent-Based Modeling in Marketing: Guidelines for Rigor.” International Journal of Research in Marketing, 28(3), 181–193.
4. Hill Function Applications: Hofmann, F., & Lehmann, T. (2021). “The Hill Function and Its Implementation: A Strategic Overview.” Marketing Science.
5. Hierarchical Bayesian Models: Gelman, A., et al. (2013). Bayesian Data Analysis (3rd ed.). CRC Press.

Open-Source Tools

• PyMC: https://www.pymc.io/ — Probabilistic programming for Bayesian inference
• Meta Robyn: https://github.com/facebookexperimental/Robyn — MMM framework from Meta
• Google LightweightMMM: https://github.com/google/lightweight_mmm — Bayesian MMM in JAX
• Papilon: https://github.com/Bodhi8/papilon — This framework

Industry Resources

• Meta Marketing Science: https://www.facebook.com/business/measurement
• Google Marketing Platform: https://marketingplatform.google.com/about/analytics/
• MMA (Mobile Marketing Association): MMM best practices and standards

About the Author

Brian Curry is a Kansas City based Data Scientist, Artificial Intelligence expert and founder of Vector1 Research, a part time passion project, where he explores the frontier of Artificial Intelligence, SEO, Agentic Systems, Machine Learning, Economic Systems, Marketing Modeling, and Econometrics. His work focuses on building intelligent systems that reveal the causal structures driving performance across marketing, content, and economics and the creator of Papilon, an open-source framework for complex system optimization, Daedalus, a platform for causal and economic analysis, and the scientist behind Memory-Node Encapsulation (MNE): A Revolutionary Data Structure for Artificial Episodic Memory and Cognitive AI Reasoning.

© 2025 Brian Curry. All rights reserved. Original models, frameworks, and conceptual structures are proprietary intellectual property. Do not reproduce without permission.

Citation

If you use Papilon in your research or practice, please cite:

Curry, B. (2025). End-to-End Media Mix Optimization with Papilon: 
A Rigorous Simulation-Driven Approach to Causal Marketing Measurement. 
Vector1 Research Technical Report.

This article is intended for technical practitioners, data scientists, and researchers working in marketing analytics, econometrics, and causal inference. For implementation support, visit the Papilon GitHub repository or contact the author.