# Bayesian variable selection — does it really work?

I thought I might toy with some Bayesian variable selection, following a nice blog post and the linked papers therein. I wrote a program in rjags (where I am quite a rookie) and fetched price data for Exxon Mobil, along with some things that are unlikely to explain its returns (e.g. palladium prices) and other things that should be highly correlated (like the SP500).

Running lm(), we see that there strong evidence of an overparameterized model, but that palladium should definitely be excluded:

Call:
lm(formula = Exxon ~ 0 + SP + Palladium + Russell + OilETF +
EnergyStks, data = chkr)

Residuals:
Min         1Q     Median         3Q        Max
-1.663e-03 -4.419e-04  3.099e-05  3.991e-04  1.677e-03

Coefficients:
Estimate Std. Error t value Pr(>|t|)
SP          0.51913    0.19772   2.626 0.010588 *
Palladium   0.01620    0.03744   0.433 0.666469
Russell    -0.34577    0.09946  -3.476 0.000871 ***
OilETF     -0.17327    0.08285  -2.091 0.040082 *
EnergyStks  0.79219    0.11418   6.938 1.53e-09 ***


After converting to returns, I tried running a simple model like this

  model {
for (i in 1:n) {
mean[i]<-inprod(X[i,],beta)
y[i]~dnorm(mean[i],tau)
}
for (j in 1:p) {
indicator[j]~dbern(probindicator)
betaifincluded[j]~dnorm(0,taubeta)
beta[j] <- indicator[j]*betaifincluded[j]
}
tau~dgamma(1,0.01)
taubeta~dgamma(1,0.01)
probindicator~dbeta(2,8)
}


but I found that, pretty much regardless of the parameters to the chosen gamma distributions, I got pretty nonsensical answers, such as an unvarying 20% inclusion probability for each variable.

I also got tiny, tiny regression coefficients, which I am willing to tolerate since this supposed to be a selection model, but that still seemed weird.

                              Mean        SD  Naive SE Time-series SE
SP         beta       -4.484e-03   0.10999  0.003478       0.007273
Palladium  beta        1.422e-02   0.16646  0.005264       0.011106
Russell    beta       -2.406e-03   0.08440  0.002669       0.003236
OilETF     beta       -4.539e-03   0.14706  0.004651       0.005430
EnergyStks beta       -1.106e-03   0.07907  0.002500       0.002647
SP         indicator   1.980e-01   0.39869  0.012608       0.014786
Palladium  indicator   1.960e-01   0.39717  0.012560       0.014550
Russell    indicator   1.830e-01   0.38686  0.012234       0.013398
OilETF     indicator   1.930e-01   0.39485  0.012486       0.013229
EnergyStks indicator   2.070e-01   0.40536  0.012819       0.014505
probindicator  1.952e-01   0.11981  0.003789       0.005625
tau            3.845e+03 632.18562 19.991465      19.991465
taubeta        1.119e+02 107.34143  3.394434       7.926577


Is Bayesian variable selection really that bad/sensitive? Or am I making some glaring error?

• Pardon my ignorance; but what was the evidence for overfitting that you refer to? – curious_cat Mar 6 '13 at 5:06
• You should explain which variables are which in the second output. I've used Bayesian variable selection on a variety of problems and in a number of situations (including regression) it usually works reasonably well. But your results - especially the estimates - look weird to me. – Glen_b Mar 6 '13 at 8:58
• @curious_cat The evidence for overfitting is, e.g., in the negative coefficient between Exxon (an oil company) and oil price. It arises because I have deliberately made this model victim to multicollinearity. (Perhaps "overfitting" is the wrong word to describe it -- I suppose overparameterized is more accurate). – Brian B Mar 6 '13 at 13:54
• @BrianB Does that coefficient become positive if you drop all explanatory variables except oil? Just curious. – curious_cat Mar 6 '13 at 14:16
• @curious_cat Yes, certainly (roughly 0.7). This is a classic case of multicollinearity (another ugly word). – Brian B Mar 6 '13 at 14:33

## 3 Answers

In the BUGS code, mean[i]<-inprod(X[i,],beta) should be mean[i]<-inprod(X[i,],beta[]).

Your priors on tau and taubeta are too informative.

You need a non-informative prior on betaifincluded, use e.g. a gamma(0.1,0.1) on taubeta. This may explain why you get tiny regression coefficients.

• Thanks for noting that. Unfortunately it did not improve the situation. – Brian B Mar 11 '14 at 3:18

It does work, but you gave all the variable inclusion indicators the same underlying distribution.

  model {
for (i in 1:n) {
mean[i]<-inprod(X[i,],beta)
y[i]~dnorm(mean[i],tau)
}
for (j in 1:p) {
indicator[j]~dbern(probindicator[j])
probindicator[j]~dbeta(2,8)
betaifincluded[j]~dnorm(0,taubeta)
beta[j] <- indicator[j]*betaifincluded[j]
}
tau~dgamma(1,0.01)
taubeta~dgamma(1,0.01)

}


might work better with a limited number of variables.

• Trying this recipe worked no better, at least at 10000 samples. – Brian B Mar 11 '14 at 3:19

If you used log returns, then you made a slightly biasing error but if you used future value divided by present value then your likelihood is wrong. Actually, your likelihood is wrong in either case. It is wrong enough to matter.

Consider that a statistic is any function of the data. Returns are not data, they are transformations of data. They are a future value divided by a present value. Prices are data. Prices must have a distribution function, but the distribution function for returns must depend solely on the nature of the prices.

For securities in a double auction, there is no "winner's curse." The rational behavior is to bid your expectation. With many buyers and many sellers, the limit book should converge to the normal distribution as it is a distribution of expectations. So $p_t$ should be normally distributed. Also $p_{t+1}$ should be normally distributed. Therefore returns should be the ratio of $$\frac{p_{t+1}}{p_t}-1.$$

The likelihood function for your regression should have been $$\frac{1}{\pi}\frac{\sigma}{\sigma^2+(y-\beta_1x_1-\beta_2x_2\dots-\beta_nx_n-\alpha)^2}.$$

OLS forces a best fit to the observed data even if it is the wrong solution. Bayesian methods attempt to find the data generating function through the likelihood. You had the likelihood wrong, so it could not find it.

I have a paper out on this if you need additional information.

EDIT I think you have misunderstood. If you would convert the likelihood to a density function and take the expectation, you would find that it has none. By proof by Augustin Cauchy in 1852 or maybe 1851, any form of least squares solution is perfectly imprecise. It will always fail. It isn't that you should use standard regression because the Bayesian is sensitive to the likelihood, it is that Bayes is the only available solution that is admissible, with some special exceptions for some unusual special cases.

In doing the empirical testing on this, and before I had read enough of the math, I naively thought that the Bayesian and the Frequentist solution should match. There is, approximately, a theorem that says that as the sample becomes large enough, the two will converge. I used all end of day trades in the CRSP universe from 1925-2013 to test it. That isn't what the theorem says though. I was misunderstanding the rules.

I also tried the problem in logs, and it still did not match. So I realized something, all distributions are shapes, and so I constructed a geometric solution to determine which solution was correct. I treated it as a pure geometry problem to determine which algebraic answer matched the data.

The Bayesian one matched. This lead me down a very mathematical path because I couldn't figure out why the unbiased estimator was so wrong. Just for the record, using disaggregated returns over the period 1925-2013 and removing shell companies, closed-end funds and so forth, the discrepancy between the center of location is 2% and the measure of risk is understated by 4% for annual returns. This discrepancy holds under log transformation, but for a different reason. It may be different for individual indices or subsets of the data.

The reason for the discrepancy is two-fold. The first is that the distributions involved lack a sufficient statistic. For certain types of problems, this doesn't matter. For projective purposes, such as prediction or allocation, however, they matter quite a bit. The second reason is that the unbiased estimator is always a version of the mean, but the distribution has no mean.

The density above is not a member of the exponential family as the normal or the gamma distribution is. By the Pitman–Koopman–Darmois theorem, no sufficient point statistic exists for the parameters. This implies that any attempt to create a point estimator must throw away information. This is not a problem for Bayesian solutions because the posterior is an entire density and if you did need a point estimate, you could find the predictive density and minimize a cost function over it to reduce it to a single point. The Bayesian likelihood is always minimally sufficient.

The minimum variance unbiased estimator for the above function is to keep the central 24.6% of the data, find its trimmed mean, and to discard the rest of the data. That means over 75% of the data is dropped, and the information is lost. Just a note, it might be 24.8%, as I am working from memory. You can find Rothenberg's paper at:

Rothenberg, T. J. and F. M. Fisher, and C. B. Tilanus, A Note on Estimation from a Cauchy Sample, Journal of the American Statistical Association, 1964, vol 59(306), pp. 460-463

The second issue was surprising to me. Until I worked through the geometry, I didn't realize what the cause was. Returns are bound on the bottom at -100%. This shifts the median by 2% and the interquartile range is shifted by 4% although the half-mass is still at the same points. The half-mass is the proper measure of scale, but the half-width is not. If there were no truncation, then the half-width and the half-mass would be at the same points. Likewise, the median and the mode would remain at the same point. The median is the return for the mean actor or at least the mean trade. As such, it is always the location of the MVUE and the log mean.

The correct understanding of the theorem is that all Bayesian estimators are admissible estimators. Frequentist estimators are admissible estimators if one of two conditions obtain. The first is that in every sample, the Frequentist and the Bayesian solution are identical. The second is that if the limiting solution of the Bayesian method matched the Frequentist solution, then the Frequentist solution is admissible.

All admissible estimators converge to the same solution once the sample size is large enough. The Frequentist estimator presumes that its model is the true model and the data is random. The Bayesian assumes the data is true, but the model is random. If you had an infinite amount of data, then the subjective model must converge to reality. If you had an infinite amount of data, but the wrong model, then the Frequentist model will converge to reality with probability zero.

In this case, the Bayesian solution, under reasonable priors, will always stochastically dominate any Frequentist estimator because of the truncation and the loss of information to create the estimator.

In logs, the likelihood function is the hyperbolic secant distribution. It has a finite variance, but no covariance. The covariance matrix found using OLS is an artifact of the data and does not point to a parameter that exists in the underlying data. As with the raw form, nothing in the log form covaries, but nothing is independent either. Instead, a far more complex relationship exists that violates the definition of covariance, but in which they can comove.

Markowitz and Usman almost found it in their work on distributions, but the hyperbolic secant distribution isn't in a Pearson family and they misinterpreted the data by not noticing that when you change a distribution from raw data to log data you also change its statistical properties. They basically found this out but missed it because they had no reason to look for it and they didn't realize the unintended consequences of using logs.

I don't have the Markowitz and Usman cite with me where I am at, but they did one of the few very good jobs at estimating the distribution that are out there.

In any case, I don't use JAGS. I have no idea how to do it. I code all my MCMC work by hand.

I have a paper that is far more complete and accurate on this topic at:

Harris, D.E. (2017) The Distribution of Returns. Journal of Mathematical Finance, 7, 769-804.

It will provide you with a method to construct distributions for any asset or liability class, also accounting ratios.

I was wordy, but I could see you were misunderstanding the relationship between Bayes and the Pearson-Neyman methods. You had them reversed. Bayes always works, but you are trapped with a prior density that will perturb your solution. With a proper prior you are guaranteed a biased estimator and for this type of likelihood function, I believe you must use a proper prior to guarantee integrability to unity. Frequentist methods are fast and usually work. They are unbiased, but may not be valid.

• The difference between logarithms and price returns at the daily level is about 1 part in 500 (I've taken a 50 basis point return as baseline). It is also worth noting that daily returns, measured either way, are quite far from normally distributed (regardless of economic theories to the contrary). In any case if the algorithm is sensitive to these issues, then it is useless for practical purposes. I would like to know how the likelihood function you cite maps to JAGS code -- can you elaborate? – Brian B Nov 30 '17 at 4:01
• @BrianB I edited the answer. You are making a mistake of thinking that I made. It took me forever to figure out why I was wrong. – Dave Harris Nov 30 '17 at 5:53