Smoothly Mixing Regressions

Definition

A smoothly mixing regression (SMR) model is a Bayesian approach to inference about the full conditional distribution $p(y|x)$ — not just the conditional mean — by modeling $y$ as a finite mixture of normal regressions whose state probabilities are smooth (continuous) functions of observed covariates. The key extension over a standard finite mixture is that mixing weights depend on $x$ through a multinomial probit model, so the entire shape of the conditional distribution (location, scale, skewness, kurtosis) can vary flexibly with covariates.

Key Ideas

Why full conditional distribution? Quantile regression models only specific quantiles and can violate monotonicity in finite samples. SMR recovers the entire $p(y|x)$ , from which any quantile, tail probability, or functional $\int f(y)\,p(y|x)\,dy$ can be computed as a by-product.
State probabilities from multinomial probit: A latent vector $\tilde w_t \sim \mathcal{N}(\Gamma z_t, I_m)$ determines state $\tilde s_t = j$ iff $\tilde w_{tj} \geq \tilde w_{ti}$ for all $i$ . This gives $P(\tilde s_t = j | z_t)$ as a continuous function of covariates $z_t$ . The $I_m$ variance identification avoids the multimodal likelihood of a full multinomial probit (MNP) (Keane 1992); multinomial logit produces nonconjugate conditionals with inferior MCMC mixing.
Five model configurations (Table 1): A–E vary which covariate vectors $(u_t, v_t, z_t)$ $(u_{t}, v_{t}, z_{t})$ are active:
- A: $q=1$ , $k>1$ , $p=1$ — fixed mixture of disturbances in regression
- B: $q=1$ , $k=1$ , $p>1$ — fixed mixture of linear regressions
- C: $q>1$ , $k=1$ , $p=1$ — smooth mixture of normal distributions (for asset returns)
- D: $q>1$ , $k>1$ , $p=1$ , $u_t=z_t$ — smooth mixture of disturbances
- E: $q>1$ , $k=1$ , $p>1$ , $v_t=z_t$ — smooth mixture of linear regressions (best for earnings)
Variance decomposition: $\text{Var}(y_t|\tilde s_t=j) = \sigma^2\sigma_j^2$ ; common scale $\sigma^2$ governs overall magnitude, state-specific $\sigma_j^2$ govern shape and relative scale.
Grid-based priors: Polynomial basis functions map two covariates $(a,b)$ into $u_t, v_t, z_t$ . Priors $C\beta \sim \mathcal{N}(\iota_r\mu, \tau_\beta^2 r I_r)$ where $C$ is an $r$ -point evaluation grid — each grid point has the weight of roughly one observation.
Label switching: States are exchangeable because the conditional distribution $p(y|x)$ is invariant to permutation of state labels. Functions of interest (quantiles, moments) depend only on $p(y|x)$ , not on which label is attached to which component.

How It Works

Markov Chain Monte Carlo (MCMC) Algorithm (Five-Block Gibbs)

Given parameters $(\beta, A, \Gamma^*, \sigma^2, \sigma_j^2)$ and latent states $(\tilde W, \tilde s)$ :

$(\beta, A)$ block: Conditional on $\tilde s$ , this is a generalized Gaussian regression — conjugate Normal posterior.
$\sigma^2$ block: Inverse-gamma conjugate, updated from residuals $y_t - \alpha_{\tilde s_t}'v_t - \beta'u_t$ .
$\sigma_j^2$ block ( $j=1,\ldots,m$ ): Conditionally independent inverse-gamma draws, each from observations assigned to state $j$ .
$\Gamma^*$ block: Conditional on $\tilde W^*$ , the rows $\gamma_j^*$ are independent Normal — fully conjugate.
$\tilde W$ block (Metropolis-within-Gibbs): For each $t$ , draw $\tilde w_t^{\text{prop}} \sim \mathcal{N}(\Gamma z_t, I_m)$ ; the proposal selects state $j^*$ (the max-coordinate index). Accept with probability $\min\left\{\frac{\sigma_{j^*}^{-1}\exp[-(y_t - \alpha_{j^*}'v_t - \beta'u_t)^2/(2\sigma^2\sigma_{j^*}^2)]}{\sigma_j^{-1}\exp[-(y_t - \alpha_j'v_t - \beta'u_t)^2/(2\sigma^2\sigma_j^2)]},\; 1\right\}$ where $j$ is the current state. If $j^* = j$ , acceptance is automatic.

The algorithm runs 12,000 iterations; first 2,000 discarded; every 100th used for inference. One–two minutes on desktop hardware for $T \approx 2{,}500$ .

State Probability Computation

$P(\tilde s_t = j | \Gamma, z_t) = \int_{-\infty}^{\infty} \phi(y - \gamma_j'z_t)\prod_{i \neq j}\Phi(y - \gamma_i'z_t)\,dy$

This integral has only one dimension regardless of $m$ — the outer integral over the maximum. Conventional quadrature suffices; no Geweke-Hajivassiliou-Keane (GHK) simulation needed.

Posterior Quantile Computation

For quantile $q$ at covariate values $(a,b)$ : solve $P(y \leq c | a, b, \text{data}) = q$ iteratively, where

$P(y \leq c | a, b, \text{data}) \approx \frac{1}{M}\sum_{\text{iter}} \sum_{j=1}^m \Phi\!\left(\frac{c - \beta'u - \alpha_j'v}{\sigma\sigma_j}\right) P(\tilde s_t = j | \Gamma, z)$

Root-finding via inverse cumulative distribution function (CDF); computation for 450 grid points × 7 quantiles takes ~3 minutes.

Why It Matters

Provides a Bayesian alternative to quantile regression that (a) guarantees monotonicity of quantiles by construction, (b) gives the full conditional distribution rather than isolated quantile functions, and (c) enables computation of any functional $\int f(y)\,p(y|x)\,dy$ (Gini coefficients, option prices, exceedance probabilities).
Outperforms t-generalized autoregressive conditional heteroskedasticity (t-GARCH)(1,1) on S&P 500 out-of-sample predictive likelihood by 57 log-points (1990–1999 data) — a very large margin for time-series forecasting.
The single-dimension state probability integral makes the model computationally tractable without simulation-based integration of state probabilities.
Generalizes threshold autoregression (Geweke-Terui 1993) and related regime models as limiting cases: as $\|\Gamma\| \to \infty$ , state assignment becomes deterministic and piecewise, recovering Voronoi tessellations and threshold models.

Open Questions

Choosing $m$ requires cross-validation or predictive likelihoods; no fully Bayesian approach (e.g., reversible-jump or Dirichlet process) is implemented.
Multivariate $y$ raises the curse of dimensionality in specifying $(u_t, v_t, z_t)$ from lagged $y$ .
Theoretical consistency of the posterior distribution of $\int f(y)\,p(y|x|\theta)\,dy$ around its truth — asymptotic theory is open.
Comparison with Dirichlet process mixtures: Dirichlet process (DP) approaches grow $m$ endogenously but face harder computation; SMR fixes $m$ but achieves competitive performance at lower computational cost.