How to generalize shrinkage to fully-bayesian models?

Question

I've got a time series that I'm modelling as an exponential; growth rate, with the rate following a logistic distribution:

$$ y_t = e^{x_t r_t} $$ where $$ r_t = \frac{L}{1-e^{-k(x_t-x_0)}} $$

I've got priors on $L$, $k$, and $x_0$, and they're distributed $\beta$, $\Gamma_1$, and $\Gamma_2$, respectively.

I fit the model using Metropolis-Hastings/MCMC. I get a good fit to observed data. Here's the problem: when I try to validate it by holding out a few most recent days of data to see how well it does predictively, it becomes apparent that I'm overfitting horribly.

If I were doing this in a frequentist context, I'd just do a gridsearch over the parameters using look-forward cross validation. But that doesn't work for two reasons:

1. I need credible intervals
2. By definition, my problem is non-stationary.

So, how do I control the bias/variance tradeoff in a bayesian context? I had an idea that I could add a criterion in the MH algorithm that requires a proposed position to fit both the past (training) set and the future (holdout set). But that feels like polluting the training set with the testing set.

I'd appreciate ideas and references if this is a standard/solved problem.

Answer 1

To answer my own question: I don't know if there is a canonical answer, but I came up with a working solution by sort of generalizing a $L_n$ penalty by penalizing the quantiles of the marginal priors.

Here's the main function:

def do_shrinkage(pos, shrinkage):
    densities = sps.beta.pdf(pos, a = shrinkage[0], b = shrinkage[1])
    regularization_penalty = -np.sum(np.log(densities))
    return regularization_penalty

pos is a vector of quantiles, in the space of the prior, of the given MCMC iterate, for each variable. Shrinkage is a couple of beta parameters that I get from this function:

from scipy.optimize import fmin
from scipy.stats import gamma, beta
import numpy as np
def beta_from_q(l, u, quantiles_percent=0.95):
    def loss(params):
        a, b = params
        lq = (1 - quantiles_percent) / 2
        uq = 1 - lq
        return ( (beta.cdf(l, a, b) - lq)**2 + (beta.cdf(u, a, b) - uq)**2 )

    start_params = (1, 1)
    fit = fmin(loss, start_params, disp = 0)
    return fit

I loop through a bunch of $\beta$ regularizers with 95% quantiles going from $[.05,.95]$ to $[.45, .55]$, subtracting off the regularization_penalty from the posterior, and thereby making it less likely that MH will pick heavily-regularized cases. I pick the one that works best on the week-ahead holdout set.

It works pretty well! Maybe this'll be useful to someone someday. I'd love to know of other approaches people have taken to this problem.