Bayesian Inference for Trading: Probabilistic Modeling
Traditional statistical approaches provide point estimates—a single number for expected return or volatility. Bayesian inference offers a richer framework: full probability distributions that quantify uncertainty, update systematically as new data arrives, and incorporate prior knowledge. For traders navigating uncertain markets, Bayesian methods provide probabilistic forecasts, adaptive models, and principled uncertainty quantification.Understanding Bayesian Inference
Bayesian inference is based on Bayes' theorem, which describes how to update beliefs given new evidence.Bayes' Theorem
``
P(θ|D) = P(D|θ)·P(θ) / P(D)
`
Where:
- P(θ|D): Posterior (updated belief about parameter θ given data D)
- P(D|θ): Likelihood (probability of observing data given θ)
- P(θ): Prior (initial belief about θ before seeing data)
- P(D): Evidence (normalizing constant)
In practical terms:
`
Posterior ∝ Likelihood × Prior
`
Sequential Learning
Bayesian inference naturally handles sequential data:
`
Prior_today = Posterior_yesterday
`
Each day's posterior becomes next day's prior—continuous learning without catastrophic forgetting.
Key Takeaways
- Bayesian methods provide full probability distributions, not just point estimates
- Naturally incorporate prior knowledge and expert judgment
- Update beliefs systematically as new data arrives
- Quantify uncertainty in predictions and parameters
- Enable probabilistic forecasts for risk management
- Regularization through informative priors prevents overfitting
Why Bayesian Inference Matters for Trading
1. Probabilistic Forecasting
Get full distributions of future returns:
`python
import numpy as np
import pymc3 as pm
def bayesian_return_forecast(returns, forecast_horizon=21):
"""
Bayesian forecasting with uncertainty quantification
Returns full posterior distribution of future returns
"""
with pm.Model() as model:
# Prior for mean return (weakly informative)
mu = pm.Normal('mu', mu=0, sigma=0.01)
# Prior for volatility (half-normal)
sigma = pm.HalfNormal('sigma', sigma=0.02)
# Likelihood
obs = pm.Normal('obs', mu=mu, sigma=sigma, observed=returns)
# Sample posterior
trace = pm.sample(2000, tune=1000, return_inferencedata=False, progressbar=False)
# Forecast future returns
mu_posterior = trace['mu']
sigma_posterior = trace['sigma']
# Simulate future returns from posterior predictive
future_returns = np.random.normal(
mu_posterior[:, np.newaxis],
sigma_posterior[:, np.newaxis],
size=(len(mu_posterior), forecast_horizon)
)
# Cumulative returns
cumulative = np.cumprod(1 + future_returns, axis=1) - 1
return {
'forecast_dist': cumulative[:, -1], # Distribution of return at horizon
'forecast_median': np.median(cumulative[:, -1]),
'forecast_95_ci': np.percentile(cumulative[:, -1], [2.5, 97.5]),
'probability_positive': (cumulative[:, -1] > 0).mean()
}
Usage
spy_returns = np.log(spy_prices / spy_prices.shift(1)).dropna()
forecast = bayesian_return_forecast(spy_returns, forecast_horizon=21)
print(f"21-day forecast: {forecast['forecast_median']:.2%}")
print(f"95% CI: [{forecast['forecast_95_ci'][0]:.2%}, {forecast['forecast_95_ci'][1]:.2%}]")
print(f"P(positive return): {forecast['probability_positive']:.1%}")
`
2. Adaptive Parameters
Update strategy parameters as market conditions change:
`python
def bayesian_momentum_strategy(prices, lookback=20, update_freq=5):
"""
Momentum strategy with Bayesian parameter updating
Adapts to regime changes automatically
"""
returns = prices.pct_change().dropna()
signals = []
posterior_means = []
# Initialize prior
prior_mean_momentum = 0.0
prior_std_momentum = 0.5
for i in range(lookback, len(returns), update_freq):
# Historical window
window = returns.iloc[i-lookback:i]
# Calculate momentum
momentum = window.mean()
# Bayesian update
likelihood_mean = momentum
likelihood_std = window.std() / np.sqrt(len(window)) # SEM
# Posterior (conjugate normal-normal)
posterior_var = 1 / (1/prior_std_momentum2 + 1/likelihood_std2)
posterior_mean = posterior_var * (
prior_mean_momentum / prior_std_momentum**2 +
likelihood_mean / likelihood_std**2
)
posterior_std = np.sqrt(posterior_var)
# Trading signal based on posterior
if posterior_mean > 2 * posterior_std:
signal = 1 # Strong positive momentum
elif posterior_mean < -2 * posterior_std:
signal = -1 # Strong negative momentum
else:
signal = 0 # Uncertain
signals.extend([signal] * update_freq)
posterior_means.append(posterior_mean)
# Update prior for next iteration
prior_mean_momentum = posterior_mean
prior_std_momentum = posterior_std * 1.1 # Slight forgetting
return pd.Series(signals[:len(returns)], index=returns.index)
`
3. Hierarchical Models
Model relationships across assets:
`python
def bayesian_factor_model(returns_df):
"""
Bayesian factor model with hierarchical priors
Models: r_i = α_i + β_i·f + ε_i
"""
n_assets = len(returns_df.columns)
with pm.Model() as model:
# Hyperpriors for alpha distribution
mu_alpha = pm.Normal('mu_alpha', mu=0, sigma=0.01)
sigma_alpha = pm.HalfNormal('sigma_alpha', sigma=0.01)
# Hyperpriors for beta distribution
mu_beta = pm.Normal('mu_beta', mu=1.0, sigma=0.5)
sigma_beta = pm.HalfNormal('sigma_beta', sigma=0.5)
# Asset-specific parameters (hierarchical)
alpha = pm.Normal('alpha', mu=mu_alpha, sigma=sigma_alpha, shape=n_assets)
beta = pm.Normal('beta', mu=mu_beta, sigma=sigma_beta, shape=n_assets)
# Factor (market)
market_returns = returns_df.mean(axis=1) # Equal-weighted market
factor = pm.Data('factor', market_returns.values)
# Idiosyncratic volatility
sigma_idio = pm.HalfNormal('sigma_idio', sigma=0.02, shape=n_assets)
# Likelihood for each asset
for i, col in enumerate(returns_df.columns):
pm.Normal(f'returns_{i}',
mu=alpha[i] + beta[i] * factor,
sigma=sigma_idio[i],
observed=returns_df[col].values)
# Sample posterior
trace = pm.sample(1000, tune=500, return_inferencedata=False, progressbar=False)
# Extract posterior distributions
results = {
'alpha_mean': trace['alpha'].mean(axis=0),
'beta_mean': trace['beta'].mean(axis=0),
'alpha_std': trace['alpha'].std(axis=0),
'beta_std': trace['beta'].std(axis=0)
}
return results
`
Practical Trading Applications
Bayesian A/B Testing for Strategies
Compare strategies with proper uncertainty quantification:
`python
def bayesian_strategy_comparison(returns_A, returns_B):
"""
Bayesian comparison of two trading strategies
Returns probability that A outperforms B
"""
with pm.Model() as model:
# Priors for each strategy's Sharpe ratio
sharpe_A = pm.Normal('sharpe_A', mu=0, sigma=1)
sharpe_B = pm.Normal('sharpe_B', mu=0, sigma=1)
# Observed Sharpe ratios (simplified)
obs_sharpe_A = returns_A.mean() / returns_A.std() * np.sqrt(252)
obs_sharpe_B = returns_B.mean() / returns_B.std() * np.sqrt(252)
# Likelihood (treating observed Sharpe as noisy measurement)
n_A = len(returns_A)
n_B = len(returns_B)
pm.Normal('obs_A', mu=sharpe_A, sigma=1/np.sqrt(n_A), observed=obs_sharpe_A)
pm.Normal('obs_B', mu=sharpe_B, sigma=1/np.sqrt(n_B), observed=obs_sharpe_B)
# Difference
diff = pm.Deterministic('diff', sharpe_A - sharpe_B)
# Sample
trace = pm.sample(2000, tune=1000, return_inferencedata=False, progressbar=False)
# Probability A > B
prob_A_better = (trace['diff'] > 0).mean()
# Expected difference
expected_diff = trace['diff'].mean()
# 95% credible interval for difference
ci_diff = np.percentile(trace['diff'], [2.5, 97.5])
return {
'prob_A_better': prob_A_better,
'expected_sharpe_diff': expected_diff,
'ci_95': ci_diff,
'decision': 'Choose A' if prob_A_better > 0.95 else \
'Choose B' if prob_A_better < 0.05 else \
'Inconclusive'
}
`
Bayesian Online Learning
Update predictions continuously:
`python
class BayesianOnlineLearner:
def __init__(self, prior_mean=0, prior_var=1):
self.mu = prior_mean
self.var = prior_var
self.predictions = []
def update(self, observation, obs_variance):
"""
Bayesian update with new observation (Kalman filter-like)
"""
# Posterior variance
posterior_var = 1 / (1/self.var + 1/obs_variance)
# Posterior mean
posterior_mean = posterior_var * (self.mu/self.var + observation/obs_variance)
# Update
self.mu = posterior_mean
self.var = posterior_var
return posterior_mean, np.sqrt(posterior_var)
def predict(self):
"""Return current belief (mean and std)"""
return self.mu, np.sqrt(self.var)
Usage
learner = BayesianOnlineLearner(prior_mean=0, prior_var=0.01)
for return_t in spy_returns:
# Update with observed return
obs_var = 0.01 # Observation noise
posterior_mean, posterior_std = learner.update(return_t, obs_var)
# Prediction for next period
pred_mean, pred_std = learner.predict()
# Trading signal
if pred_mean > 2 * pred_std:
signal = 1 # Long
elif pred_mean < -2 * pred_std:
signal = -1 # Short
else:
signal = 0 # Flat
`
Bayesian Portfolio Optimization
Robust portfolio weights with uncertainty:
`python
def bayesian_portfolio_optimization(returns_df, target_return):
"""
Portfolio optimization with Bayesian posterior distributions of returns
Accounts for parameter uncertainty in optimization
"""
with pm.Model() as model:
# Prior for mean returns (shrunk toward zero)
mu = pm.Normal('mu', mu=0, sigma=0.01, shape=len(returns_df.columns))
# Prior for covariance (Wishart or empirical)
cov_empirical = returns_df.cov().values
cov = pm.Data('cov', cov_empirical)
# Likelihood
for i, col in enumerate(returns_df.columns):
pm.Normal(f'obs_{i}', mu=mu[i],
sigma=np.sqrt(cov_empirical[i, i]),
observed=returns_df[col].values)
# Sample posterior
trace = pm.sample(1000, tune=500, return_inferencedata=False, progressbar=False)
# Portfolio optimization using posterior samples
from scipy.optimize import minimize
weights_samples = []
for i in range(len(trace['mu'])):
mu_sample = trace['mu'][i]
# Mean-variance optimization with this sample
def portfolio_vol(w):
return np.sqrt(w @ cov_empirical @ w)
def portfolio_return(w):
return w @ mu_sample
constraints = [
{'type': 'eq', 'fun': lambda w: np.sum(w) - 1},
{'type': 'eq', 'fun': lambda w: portfolio_return(w) - target_return}
]
bounds = [(0, 1)] * len(returns_df.columns)
result = minimize(
portfolio_vol,
x0=np.ones(len(returns_df.columns)) / len(returns_df.columns),
method='SLSQP',
bounds=bounds,
constraints=constraints
)
if result.success:
weights_samples.append(result.x)
# Average weights across posterior
optimal_weights = np.mean(weights_samples, axis=0)
return pd.Series(optimal_weights, index=returns_df.columns)
`
Advanced Techniques
MCMC for Complex Models
Use Markov Chain Monte Carlo when posteriors aren't analytical:
`python
def regime_switching_bayesian(returns):
"""
Bayesian regime-switching model using MCMC
"""
with pm.Model() as model:
# Regime probabilities
p = pm.Dirichlet('p', a=np.ones(2)) # 2 regimes
# Regime-specific parameters
mu_regime = pm.Normal('mu_regime', mu=0, sigma=0.02, shape=2)
sigma_regime = pm.HalfNormal('sigma_regime', sigma=0.02, shape=2)
# Latent regime variable
regime = pm.Categorical('regime', p=p, shape=len(returns))
# Likelihood (regime-dependent)
pm.Normal('obs',
mu=mu_regime[regime],
sigma=sigma_regime[regime],
observed=returns)
# Sample
trace = pm.sample(1000, tune=500, return_inferencedata=False, progressbar=False)
return trace
`
Variational Inference for Speed
Approximate Bayesian inference for large datasets:
`python
def variational_bayes_fast(returns_df):
"""
Variational inference: faster than MCMC
Good for production systems needing real-time updates
"""
with pm.Model() as model:
# Model specification (same as before)
mu = pm.Normal('mu', mu=0, sigma=0.01, shape=len(returns_df.columns))
sigma = pm.HalfNormal('sigma', sigma=0.02, shape=len(returns_df.columns))
for i, col in enumerate(returns_df.columns):
pm.Normal(f'obs_{i}', mu=mu[i], sigma=sigma[i],
observed=returns_df[col].values)
# Variational inference (ADVI)
approx = pm.fit(n=10000, method='advi')
# Sample from variational posterior
trace = approx.sample(1000)
return trace
`
Bayesian Model Selection
Compare models using Bayes factors:
`python
def bayesian_model_comparison(returns, models):
"""
Compare multiple models using marginal likelihood
models: list of PyMC3 models
"""
import arviz as az
model_traces = {}
model_scores = {}
for name, model in models.items():
with model:
trace = pm.sample(1000, tune=500, return_inferencedata=True)
# WAIC (Widely Applicable Information Criterion)
waic = az.waic(trace)
model_traces[name] = trace
model_scores[name] = waic.waic
# Best model (lowest WAIC)
best_model = min(model_scores.items(), key=lambda x: x[1])[0]
return {
'best_model': best_model,
'scores': model_scores,
'traces': model_traces
}
`
Implementation Best Practices
1. Choosing Priors
Select informative but not overpowering priors:
`python
def construct_prior(historical_mean, historical_std, confidence=0.5):
"""
Construct weakly informative prior from historical data
confidence: 0 = uninformative, 1 = very informative
"""
# Prior mean = historical mean
prior_mean = historical_mean
# Prior std reflects confidence
# Low confidence → wider prior
prior_std = historical_std / (confidence + 0.1)
return {
'prior_distribution': 'Normal',
'prior_mean': prior_mean,
'prior_std': prior_std,
'rationale': f'{confidence*100:.0f}% confidence in historical estimate'
}
`
2. Convergence Diagnostics
Ensure MCMC has converged:
`python
def check_mcmc_convergence(trace):
"""
Check MCMC convergence using R-hat and effective sample size
"""
import arviz as az
# Convert to InferenceData if needed
if not isinstance(trace, az.InferenceData):
trace = az.from_dict(trace)
# R-hat (should be < 1.01)
rhat = az.rhat(trace)
# Effective sample size (should be > 100 per chain)
ess = az.ess(trace)
# Summary
summary = az.summary(trace)
diagnostics = {
'rhat_ok': (summary['r_hat'] < 1.01).all(),
'ess_ok': (summary['ess_bulk'] > 100).all(),
'summary': summary
}
return diagnostics
`
3. Posterior Predictive Checks
Validate model fit:
`python
def posterior_predictive_check(model, trace, observed_data):
"""
Check if model can reproduce observed data characteristics
"""
with model:
# Generate posterior predictive samples
post_pred = pm.sample_posterior_predictive(trace, samples=500)
# Compare distributions
from scipy import stats
# KS test: observed vs posterior predictive
ks_statistic, p_value = stats.ks_2samp(
observed_data.flatten(),
post_pred['obs'].flatten()
)
return {
'ks_statistic': ks_statistic,
'p_value': p_value,
'model_adequate': p_value > 0.05
}
``