Bayesian Analysis: Introduction

0. Introduction. After some pause, we start a new series of posts about Bayesian inference with regard to our financial modeling. This modeling involves many linear regressions fit using ordinary least squares (OLS). We tried very hard to make these residuals to be multivariate Gaussian independent identically distributed for stock and bond markets. They have zero mean vector but nontrivial covariance matrix. In many cases, as we can see from our previous posts, we succeeded. As many readers know, Bayesians introduce prior distribution of coefficients and standard errors for regressions. Then they use the observed financial data to infer the posterior distributions for said coefficients and standard errors.

1. Methodology. After making Bayesian inference, we can make Monte Carlo simulations of time series of market factors and asset prices. These time series have multivariate Gaussian residuals (innovations). However, we shall do this differently from the classic frequentist approach.

Frequentist approach: Fit the values of coefficients and standard errors, which are the same for all simulations. Then run as many simulations as you wish.

Bayesian approach: Fix the number of simulations, for example 10000. From the posterior distributions, simulate 10000 values of coefficients and standard errors. For each such simulation out of these 10000, fix these values of coefficients and standard errors, and simulate the time series with the same fixed standard errors, using regression equations with these same coefficients. We have two versions of this Bayesian approach:

Bayesian-lite: We treat the covariance matrix of residuals as fixed, as in the classic approach. We use these to simulate residuals (innovations) as multivariate Gaussian. But we use Bayesian inference only for regression coefficients.

Bayesian-full: We treat the correlation matrix of residuals as fixed, as in the classic approach. But we use Bayesian inference upon standard errors of residuals. Conditional on these standard errors, we run Bayesian inference upon regression coefficients. Together with the correlation matrix, this gives us the covariance matrix to simulate residuals (innovations) as multivariate Gaussian.

Actually, the true full Bayesian inference would involve making inference and imposing priors upon the entire covariance matrix, not just its diagonal elements which are squares of standard errors. But this would involve complicated matrix-valued Wishart distributions (multivariate versions of  \chi^2 distributions). Furthermore, in the frequentist approach in all of the above research we fit each regression using OLS as univariate. We did not use combined OLS for all regressions. If we were concerned about these correlated residuals, we would have a combined OLS/GLS model. But we did not do this, instead computing the empirical covariance matrix after fitting these regressions in an individual and univariate way.

Effectively, we did not much care about interaction between regressions in the frequentist approach already. Thus it does not make sense for us to care about this in the Bayesian approach. We simply consider the correlation matrix to be fixed and estimated in the classic frequentist sense. This is not ideal, but we have only limited time and enthusiasm.

2. A Simple Model. Consider only S&P returns with only one factor: volatility. The data is the same as in the previous posts, nominal total returns and empirical annual volatility 1928–2025. This is taken from the standard data file. We have:

  •  Q(t) = \theta + \gamma V(t) + V(t)Z(t) for returns (98 data points)
  •  \ln V(t) = \alpha + \beta \ln V(t-1) + W(t) for volatility (97 data points)

Assume for now that  Z(t) \sim \mathcal N(0, \sigma^2) and  W(t) \sim \mathcal N(0, \rho^2) are independent Gaussian random variables for each  t. Also,  (Z(1), W(1)), \ldots, (Z(T), W(T)) Point estimates are: and

  •  \alpha = 0.84785
  •  \beta = 0.620146
  •  \theta = 0.227336
  •  \gamma = -0.012476
  •  \sigma = 0.0162389
  •  \rho = 0.364353

By the way, in this case the two series of residuals have (close to) zero correlation, and the Jarque-Bera test gives fail to reject results for each of them. Thus we assume they are bivariate Gaussian with independent components. Thus the discussion about interactions between regressions from Section 1 is not applicable here.

3. Noninformative Priors Controversy. These are priors which assume zero known information about coefficients. These are very different from the usual informative priors, which are based on some existing information or beliefs and are later updated to get posteriors using new data. Usually, we take Jeffrey’s noninformative priors. To be fair, these priors are improper: They do not represent real probability measures, but infinite measures, and their densities integrate up to infinity. But, as we see later, their corresponding posteriors are proper: They are real probability measures. This paradox is discussed in textbooks.

We must also mention that some Bayesian statisticians consider it unacceptable to use Jeffrey’s priors. They insist that any prior must necessarily include some existing beliefs, otherwise this contradicts the very spirit and meaning of Bayesian inference. However, sometimes we do have little information. In this case, weakly informative priors are needed. For example, they could be uniform on some fixed but large interval. Or they could be exponential but with large means. Different coefficients can have independent weakly informative priors. Again, this is discussed in textbooks.

Personally, I do not feel qualified to judge this. However, Jeffrey’s priors are attractive to me because they have very simple posteriors, as we see below. Furthermore, the linear regression works in the following way: Credible intervals for given percentage level (for example, [2.5%, 97.5%]) for each coefficient coincide with confidence intervals (from the classic frequentist approach) for the same confidence percentage level and for the same coefficient. This remarkable matching of frequentist and Bayesian statistics appeals to me since this is a way to take care of the uncertainty in coefficient estimates (especially with small data) but stay objective and not use any prior beliefs. This is why this approach is called objective Bayes inference.

4. Jeffrey’s Priors and Posteriors. Take the univariate but maybe multiple linear regression  \mathbf{y} = \mathbf{X}\mathbf{\beta}| + \delta for  \delta = (\delta_1, \ldots, \delta_n) with  \delta_i \sim \mathcal N(0, \sigma^2) IID. Here  \mathbf{y} = (y_1, \ldots, y_n) is the output and  \mathbf{X} is an  n\times m design matrix. We assume there are  m regression factors, and  n > m data points.

Frequentist approach: The OLS estimate for the coefficient is

 \hat{\beta}_{OLS} = (\mathbf{X}^T\mathbf{X})^{-1}\mathbf{X}^T\mathbf{y}

and the unbiased estimate for the variance (square of the standard error) is

 s^2 = \frac1{n-m}|\!|\mathbf{y} - \mathbf{X}\hat{\beta}_{OLS}|\!|^2

Consider both Bayesian cases mentioned in Section 1. Of course, we do not perform any derivations, referring to textbooks.

Bayesian-lite: Consider the linear regression with fixed standard error  \sigma estimated using classic OLS as  \hat{\sigma} = s (the square root of the empirical variance) which, in turn, gives us the estimate  \hat{\beta}_{OLS}. Then the prior is flat:  \pi(\beta) \propto 1 and the posterior is  \mathcal N(\hat{\beta}_{OLS}, (\mathbf{X}^T\mathbf{X})^{-1}\sigma^2) multivariate normal.

Bayesian-full: Consider the linear regression with unknown standard error  \sigma and related unknown precision  p = \sigma^{-2} which is more convenient for Bayesian analysis. The prior  \pi(p) \propto 1/p is also improper, since the integral of this function from zero to infinity is infinite. But the corresponding posterior is scaled Gamma, distributed as  s^{-2}Y where  Y \sim \Gamma((n-m)/2, (n-m)s^2/2).

Next, the conditional prior for  \beta given precision  p (or, equivalently, the standard error  sigma = p^{-0.5} ) is flat, as in Bayesian-lite. The conditional posterior for  \beta given precision  p is multivariate Gaussian, similarly to the Bayesian-lite version:  \beta \mid p \sim \mathcal N(\hat{\beta}_{OLS},  (\mathbf{X}^T\mathbf{X})^{-1}/p). In fact, the inverse covariance matrix (the precision matrix) is  p\mathbf{X}^T\mathbf{X}.

Here we define the Gamma distribution  \Gamma(\alpha, \lambda) with shape  \alpha and rate  \lambda given by density  (1/\Gamma(\alpha))\lambda^{\alpha}y^{\alpha - 1}e^{-\lambda y} and which has mean  \alpha/\lambda.

5. Bayesian Simple Linear Regression. Trying a particular case  y_i = a + bx_i + \delta_i for  i = 1, \ldots, n. This is important in our simple model in this post, since both regressions in Section 2 are, in fact, simple. We get in the notation of the previous section:  m = 2, and the design matrix is

 \mathbf{X} = \begin{bmatrix} 1 & x_1\\ 1 & x_2 \\ \vdots & \vdots \\ 1 & x_n \end{bmatrix}

 \mathbf{X}^T\mathbf{X} = \begin{bmatrix} n & x_1 + \ldots + x_n \\ x_1 + \ldots + x_n & x_1^2 + \ldots + x_n^2 \end{bmatrix}

 s^2 = \frac1{n-2}\sum\limits_{i=1}^n(y_i - \hat{a}_{OLS} - \hat{b}_{OLS}x_i)^2

6. Simulation Results. We take the time horizon of 30 years. We start with volatility  V(0) = 20 which is close to the long-term average. We simulate 10000 times and compare classic, Bayesian-lite, and Bayesian-full approaches. We then consider the mean return for each simulated path out of these 10000. We take the mean, median, standard deviation, 10%, 30%, 70%, 90% quantiles, for each of these 10000 simulated paths.

7. Withdrawal Rule. Try to modify this by withdrawing 4% of the initial capital at the first year (classic withdrawal rule), and increasing this amount by 4% to adjust for inflation. Remember, returns were nominal, not real.

Published by


Leave a comment