Using equity premium survey data to estimate future wealth

Abstract

We present the first systematic methods for combining different experts’ responses to equity premium surveys. These techniques are based on the observation that the survey data are approximately gamma distributed. This distribution has convenient analytical properties that enable us to address three important problems that investment managers must face. First, we construct probability density functions for the future values of equity index tracker funds. Second, we calculate unbiased and minimum least square error estimators of the future value of these funds. Third, we derive optimal asset allocation weights between equities and the risk-free asset for risk-averse investors. Our analysis allows for both herding and biasedness in expert responses. We show that, unless investors are highly uncertain about expert biases or forecasts are very highly correlated, many investment decisions can be based solely on the mean of the survey data minus any expected bias. We also make recommendations for the design of future equity premium surveys.

Appendix

In this appendix, we present further details of the Kulkarni and Powar (2010) OPNAM method for estimating the density function of the mean of a gamma distribution conditional on a finite sample, the Winkler (1981) method for combining correlated expert forecasts, and the application of these two methods to the problem at hand. We also, by simulation, demonstrate the accuracy of these two techniques.

1.1 Transforming the data

The OPNAM method relies on the fact that if the underlying data $\lambda _{i}\sim \varGamma \left( \alpha ,\beta \right)$ then $Y_{i}=\lambda _{i}^{p}$ will be approximately normally distributed for $p=0.246$ provided that $\alpha >1.5$. The first four moments of $Y_{i}$ are known in closed form—we derive the first two here and note that the equivalent equations for skewness and kurtosis are given on page 437 of Kulkarni and Powar (2010). The pdf of a gamma distribution is given by:

$$\begin{aligned} f\left( \lambda _{i}\right) \sim \varGamma \left( \alpha,\beta \right) =\frac{\beta ^{\alpha }}{\varGamma \left( \alpha\right) }\lambda _{i}^{\alpha -1}e^{-\beta \lambda _{i}} \end{aligned}$$

(24)

So, the expectation of $Y_{p}=E\left[ \lambda _{i}^{p}\right]$ is given by

$$\begin{aligned} Y_{p}&= \int _{0}^{\infty }\lambda _{i}^{p}\frac{\beta ^{\alpha }}{\varGamma \left( \alpha\right) }\lambda _{i}^{\alpha -1}e^{-\beta \lambda _{i}}d\lambda _{i}=\int _{0}^{\infty }\frac{\beta ^{\alpha }}{\varGamma \left(\alpha\right)} \lambda _{i}^{p+\alpha -1}e^{-\beta \lambda _{i}}d\lambda _{i} \\&= \frac{\varGamma \left( \alpha+p\right) }{\beta ^{p}\varGamma \left( \alpha\right) } \int _{0}^{\infty }\frac{\beta ^{p+\alpha }}{\varGamma \left( \alpha+p\right) } \lambda _{i}^{p+\alpha -1}e^{-\beta \lambda _{i}}d\lambda _{i} \nonumber \end{aligned}$$

(25)

The transformation in the last line is ‘trivial’ in the sense that it is just multiplying both the numerator and denominator by $\beta ^{p}\varGamma \left( \alpha+p \right)$ and then rearranging. The purpose of this trivial transformation is to ensure that the term within the integral is identical to the probability density function of a gamma distribution with parameters $\alpha +p$ and $\beta$. Clearly as this is the integral over the whole support of well-defined pdf, it must equal one. Therefore:

$$\begin{aligned} Y_{p}=\frac{\varGamma \left( \alpha+p\right) }{\beta ^{p}\varGamma \left( \alpha\right) } =\left( \frac{\alpha }{\beta }\right) ^{p}\frac{\varGamma \left( \alpha+p\right) }{ \alpha^{p}\varGamma \left( \alpha\right) }=E\left[ \lambda _{i}\right] ^{p}\frac{\varGamma \left( \alpha+p\right) }{\alpha^{p}\varGamma \left( \alpha\right) } \end{aligned}$$

(26)

and this is Eq. 2 in the body of this paper.

The cross-sectional variance of the forecasts is given by $\sigma _{Y}^{2}=E \left[ \lambda _{i}^{2p}\right] -E^{2}\left[ \lambda _{i}^{p}\right]$, which follows directly from the preceding argument:

$$\begin{aligned} \sigma _{Y}^{2}=\left( \frac{1}{\beta ^{p}\varGamma \left( \alpha \right) } \right) ^{2}\left[ \varGamma \left( \alpha +2p\right) \varGamma \left( \alpha\right) -\varGamma ^{2}\left( \alpha +p\right) \right] \end{aligned}$$

(27)

Similar expressions follow for the skewness and excess kurtosis and Kulkarni and Powar note that these depend on $\alpha$ and $p$ but not $\beta$. Therefore the choice of $p$ depends on $\alpha$ only. For $\alpha >1.5$, they select $p=0.246$ so as to get the skewness and excess kurtosis jointly close to zero. Their method is an approximation, though, and therefore the $Y_{i}s$ are only approximately normally distributed even if the $\lambda _{i}s$ are perfectly gamma distributed.

A key feature to notice is that this transformation from gamma to normal does not depend on the generating process for the underlying $\lambda _{i}s$. Therefore this applies equally well to correlated as uncorrelated forecasts. To explicitly demonstrate this, we simulated correlated gamma distributed random variables as follows. First, using a Cholesky decomposition approach, we generate 131 normally distributed $N(0,1)$ and exponentially correlated random variables. We then convert these to correlated uniform random variables using the cumulative density function of the standard normal distribution. Finally, we take an Inverse Gamma function to convert again into correlated gamma distributed random variables, $\lambda _{i}$.¹⁴ We then construct the transformations $Y_{i}=\lambda _{i}^{0.246}$ and calculate the first four moments of the transformed forecasts.

In Table 4, we report our results. In Panel A we present information on the initial correlation structure that is put into the Cholesky decomposition and the estimated correlations between $\lambda _{i}$ and $\lambda _{i+1}$ and $Y_{i}$ and $Y_{i+1}$ across the simulations (where $i=50$ is chosen at random for illustrative purposes). This clearly demonstrates that the correlation structure remains almost entirely unaffected by the transformations that we undertake—we discuss this point further below. In Panel B we present the median values of the first four moments of the $Y_{i}$s across these simulations, along with lower 2.5 % and upper 97.5 % simulated values. These are compared against the theoretical values given above and in Kulkarni and Powar (2010).

Table 4 The first four moments of the transformed forecasts, $Y_i$

For all values of $\rho$ considered, the theoretical predictions lie within the 95 % confidence interval and close to the observed median values. This shows the robustness of the power transformation to correlation in the experts’ forecasts.

1.2 Constructing confidence intervals for $\alpha /\beta$

The purpose of using the OPNAM method is to construct confidence intervals, $L_{x}^{*}$, for the population mean (rather than survey mean), $\lambda _{p}$, of expert responses. To do this, we first construct confidence intervals, $L_{x}$ for the population mean, $Y_{p}$ (as given in Eq. 26) of the transformed responses and then use Eq. 4 to construct $L_{x}^{*}$. We do this using the “Unknown $\Sigma$” method described in Winkler (1981). While we assume that the correlation between experts, $\rho ^{|i-j|}$, is known perfectly, the sample cross-sectional variance of transformed forecasts, $\widehat{\sigma }_{Y}^{2},$ will differ from the population cross-sectional variance of transformed forecasts, $\sigma _{Y}^{2}$. As $\Sigma _{u|ij}=\sigma _{Y}^{2}\rho ^{|i-j|}$, this introduces sampling error into our estimate $\widehat{\varvec{\Sigma }}_{u}$.

To overcome this problem, following Winkler (1981), we form a Bayesian prior that the “true” variance-covariance matrix is Inverse Wishart distributed with parameter value $\delta$:

$$\begin{aligned} f\left( \varvec{\Sigma }|\widehat{\varvec{\Sigma }},\delta \right) \propto \left| \varvec{\Sigma }^{-1}\right| ^{\frac{\delta +2n}{2}}\exp \left( -\frac{\delta }{2}tr\left( \varvec{\Sigma }^{-1}\widehat{\mathbf { \Sigma }}\right) \right) \end{aligned}$$

(28)

where $\left| \cdot \right|$ is a matrix determinant and $tr\left( \cdot \right)$ is a matrix trace. In this case, as discussed in the body of the paper, $Y_{p}$ is Student’s $t$-distributed conditional on the transformed sample data, with $\delta +n-1$ degrees of freedom and mean $m^{*}$ and variance $s^{*^{2}},$ as given in Eq. 3. In the case of exponential correlation, $\widehat{\varvec{\Sigma }} _{u}^{-1}$ is given in Eq. 5 and, from here, the expressions for $m^{*}$ and $s^{*^{2}}$can be simplified by noting that:

$$\begin{aligned} {\mathbf {1}}^{\prime }\widehat{\varvec{\Sigma }}_{u}^{-1}{\mathbf {1}}&= \kappa \widehat{\sigma }_{Y}^{-2}\frac{n(1-\rho )+2\rho }{1+\rho },\quad {\mathbf {Y}} ^{\prime }\widehat{\varvec{\Sigma }}_{u}^{-1}{\mathbf {1}} = \kappa \widehat{ \sigma }_{Y}^{-2}\frac{n\overline{Y}-\rho \left( n-2\right) Y_{^{R}}}{1+\rho } \nonumber \\ {\mathbf {Y}}^{\prime }\widehat{\varvec{\Sigma }}{\mathbf {Y}}&= \kappa \widehat{ \sigma }_{Y}^{-2}\frac{nV+\rho ^{2}\left( n-2\right) V_{^{R}}-2\rho \left( n-1\right) Z}{1-\rho ^{2}} \\ V&= \frac{1}{n}\sum \limits _{i=1}^{n}Y_{i}^{2},\quad V_{R}=\frac{1}{ n-2}\sum \limits _{i=2}^{n-1}Y_{i}^{2},\quad Z=\frac{1}{n-1} \sum \limits _{i=1}^{n-1}Y_{i}Y_{i+1} \nonumber \end{aligned}$$

(29)

Therefore the only complexity that arises when operationalizing Eq. 3 is identifying the parameter $\delta$.

1.3 Estimating $\delta$

To estimate $\delta ,$ we draw upon the commonly observed parallel between the Inverse Wishart distribution and the Inverse Gamma distribution. We start with the standard assumption that the relationship between the cross-sectional sample variance and cross-sectional population variance is given by a Chi-squared distribution with $v-1$ degrees of freedom; $\left( v-1\right) \left( \widehat{\sigma }_{Y}^{2}/\sigma _{Y}^{2}\right) \sim \chi _{v-1}^{2}$. We explain below how $v$ is estimated. This means that:

$$\begin{aligned} \frac{\widehat{\sigma }_{Y}^{2}}{\sigma _{Y}^{2}}\sim \frac{\chi _{v-1}^{2}}{ v-1}=\varGamma \left( \frac{v-1}{2},\frac{v-1}{2}\right) \end{aligned}$$

(30)

which is a gamma distribution with both shape and rate parameters equal to $\left( v-1\right) /2$. It then follows that $\sigma _{Y}^{2}\sim \varGamma ^{-1}\left( \frac{v-1}{2},\frac{v-1}{2}\widehat{\sigma }_{Y}^{2}\right)$, where $\varGamma ^{-1}\left( \cdot ,\cdot \right)$ is an inverse gamma function and the second parameter is now the scale parameter. Then, by the properties of the Inverse Gamma distribution:

$$\begin{aligned} Var\left[ \sigma _{Y}^{2}\right] =\frac{2(v-1)\widehat{\sigma }_{Y}^{2}}{ \left( v-3\right) ^{2}\left( v-5\right) } \end{aligned}$$

(31)

To link this to the Inverse Wishart distribution that forms the basis for the “Unknown $\varSigma$” method of Winker, we note $\sigma _{Y}^{2}=\varSigma _{u|ii}$, and therefore the variance of this variable is given by $Var\left[ \varSigma _{u|ii}\right]$. For an Inverse Wishart distribution, the relationship between the precision of our estimate of $\Sigma _{u|ii}$ and $\widehat{\sigma }_{Y}^{2}$ is given by:

$$\begin{aligned} Var\left[ \varSigma _{u|ii}\right] =\frac{2\delta \widehat{\sigma }_{Y}^{2}}{ \left( \delta -2\right) ^{2}\left( \delta -4\right) } \end{aligned}$$

(32)

and by comparison of Eqs. 31 and 32 it is clear that $\delta =v-1$.

To estimate $v$ we take a numerical approach. We run 100,000 simulations for all $\rho \in \{0,0.25,0.5,0.75,0.9\}$. In each case, we draw 131 exponentially correlated Normal random error terms $u_{i}\sim N\left( 0,\sigma _{Y}^{2}\right)$, where for simulation purposes, $\sigma _{Y}^{2}= \widehat{\sigma }_{Y}^{2}$. For each simulation, $s$, we then calculate the cross-sectional variance of the error term, $\widehat{\sigma }_{sY}^{2}$ . We then take the maximum likelihood estimator of $\left( v-1\right) /2$, the shape parameter of the gamma distribution, across the 100,000 simulations.¹⁵ The values that emerge for $v$ for $\rho \in \{0,$ $0.25,\,0.5,\,0.75,0.9\}$ are, to the nearest integer, $v=\left\{ 131\text {, }115\text {, }81\text {, }40,\text { }17\right\}$. These results are reported in more detail in Table 5. This table also compares the first four moments of $\widehat{\sigma }_{Y}^{2}/\sigma _{Y}^{2}$ against those of a $\varGamma \left( \left( v-1\right) /2,\left( v-1\right) /2\right)$ distribution.

Table 5 Estimates of the degrees of freedom of the Chi-squared distribution

1.4 Accuracy of the method

To demonstrate the accuracy of our joint use of the OPNAM method and Winkler’s “Unknown $\Sigma$” technique for combining correlated forecasts, we run a further set of simulations. Within each of 100,000 simulations, for each value of $\rho \in \{0,0.25,0.5,0.75,0.9\}$ we draw 131 values of $Y_{i}$ from an exponentially correlated normal distribution $N\left( E[Y_{i}],Var\left[ Y_{i}\right] \right)$. As discussed above, both $E\left[ Y_{i}\right]$ and $Var\left[ Y_{i}\right]$ are known in closed form as functions of $\alpha ,\beta$. For simulation purposes, we use the values of $\alpha$ and $\beta$ that are empirically estimated from Welch’s data. For each simulation, we use a combination of the Winkler (1981) “Unknown $\Sigma$” method described above with the Kulkarni and Powar method to derive one-sided upper and lower 0.5 %, 1 %, 1.5 %,$\ldots$, 9.5 % and 10 % confidence intervals for $\alpha /\beta$. We then count the proportion of simulations where the initial value of $\alpha /\beta$ lies outside the estimated confidence interval. These results are reported in Table 6.

Table 6 Examining the confidence intervals of $g_p(\lambda _p)$

As can be see, in all cases there is very close agreement between the estimated confidence level and the proportion of simulations where $\alpha /\beta$ lies outside the estimated confidence interval. This demonstrates the robustness of the method that we use to correlation in the forecasts.

1.5 The correlation structure

Finally, we turn to the correlation structure and demonstrate again that there is close agreement between the assumed correlation structure for the untransformed and transformed estimates (see also Panel A of Table 4 above). For the case of $\rho =0.9$, we simulate 131 correlated $Y_{i}s$ using the standard Cholesky decomposition approach and then construct values of $\lambda _{i}=Y_{i}^{1/0.246}$. We repeat these simulations 10,000 times. In Fig. 3 we present the correlation between $\lambda _{1}$ and $\lambda _{j}$ for $j\in [2,131]$ across the 10,000 simulations. This is compared against Corr$\left( Y_{1},Y_{j}\right)$ and $\rho ^{|1-j|}$, which is the theoretical value under exponential correlation.

Fig. 3

The correlation between raw forecasts, $Corr(\lambda _1,\lambda _j)$. From 131 exponentially correlated normally distributed random variables, $Y_{i}$, we estimate the correlation of $\lambda _{j}=Y_{j}^{1/0.246}$ against $\lambda _{1}=Y_{1}^{1/0.246}$ for all $j\in [2,131]$ across $10,000$ simulations. These correlations are then compared against the theoretical value $\rho ^{|1-j|}$ and the correlation of $Y_{j}$ with $Y_{1}$. This figure is constructed with $\rho =0.9$

It is clear that the correlation structure of the $\lambda _{i}s$ is also very close to being exponentially correlated with $\rho =0.9$ even though we are modelling errors on the transformed values $Y_{i}.$