An Equiprobable Approximation to the Bivariate Lognormal

Economic agents face risks of many kinds, which may mutually covary. A stock broker, for example, is likely to earn a salary bonus that is positively related to the performance of the stock market; if that broker also has personal stock investments, his financial wealth and labor income will be positively correlated.

The first part of this section presents a convenient (and empirically realistic) formulation in which a consumer faces two shocks (which can be interpreted as a shock to noncapital income and a shock to the rate of return) that are distributed according to a multivariate lognormal that allows for correlation between them. The second part describes a computationally simple and convenient method for approximating that joint distribution.

1Theory¶

Consider a consumer who faces both a risk to transitory noncapital income^[1]

and a risky log rate-of-return that is affected by following factors: the riskless rate $\rfree$;

a risk premium $\EpremLog$; an additional constant $\zeta$ (whose purpose will become clear below); a component that is linearly related to $\ShkMeanOneLog_{1,t+1}$; and an independent shock $\ShkMeanOneLog_{2} \sim \mathcal{N}(-0.5 \sigma^{2}_{2},\sigma^{2}_{2})$:

for some constant $\omega$. Since $(\sigma_{2}/\sigma_{1}) \omega \ShkMeanOneLog_{1,t+1}$ is the only component of $\risky_{t+1}$ that covaries with $\ShkMeanOneLog_{1,t+1}$,

Equation (2) yields a description of the return process in which the parameter $\omega$ controls the correlation between the risky log return shock and the risky log labor income shock. If $\omega = 0$ the processes are independent.

Now we want to find the value of $\zeta$ such that the mean risky return is unaffected by $\sigma^{2}_{1}$ (so that we will be able to understand clearly the distinct effects of labor income risk, the independent component of rate-of-return risk $\sigma^{2}_{2}$, and the correlation between labor income risk and rate-of-return risk, $\omega$). Thus, we want to find the $\zeta$ such that

regardless of the values of $\sigma^{2}_{1}$ and $\sigma^{2}_{2}$. We therefore need:

Using standard facts about lognormals (cf. Math Facts), and for convenience defining $\hat{\omega}= (\sigma_{2}/\sigma_{1}) \omega$, we have

Because the two shocks are independent by construction, the joint distribution takes a convenient form. Stacking $\ShkMeanOneLog_{1,t+1}$ and $\ShkMeanOneLog_{2,t+1}$ into a vector $\vec{\ShkMeanOneLog}$, we can write

where $\vec{\mu}= \{-0.5 \sigma^{2}_{1},-0.5 \sigma^{2}_{2}\}^{\prime}$ and the covariance matrix $\Sigma$ is diagonal, carrying $\sigma^{2}_{1}$ and $\sigma^{2}_{2}$.^[2]

2Computation¶

A key step in the computational solution of any model with uncertainty is the calculation of expectations. Writing $\tilde{\ShkMeanOne}_{1} \equiv \tilde{\ShkMeanOne}_{1,t+1}$ and $\tilde{\Risky} \equiv \Risky_{t+1}$ and $\Ex[\bullet] = \Ex_{t}[\bullet_{t+1}]$, the expectation of some function $\hFunc$ that depends on the realization of the risky return $\tilde{\Risky}$ and the labor income shock is:

where $\FFunc(\tilde{\ShkMeanOne}_{1},\tilde{\Risky})$ is the joint cumulative distribution function. Standard numerical computation software can compute this double integral, but at such a slow speed as to be almost unusable.

Computation of the expectation can be massively speeded up by advance construction of a numerical approximation to $\FFunc(\tilde{\ShkMeanOne}_{1},\tilde{\Risky})$.

Such approximations generally take the approach of replacing the distribution function with a discretized approximation to it; appropriate weights $w_{i,j}$ are attached to each of a finite set of points indexed by $i$ and $j$, and the approximation to the integral is given by:

where the $\hat{\ShkMeanOne}_{1}$ and $\hat{\Risky}$ matrices contain the conditional means of the two variables in each of the $\{i,j\}$ regions. Various methods are used for constructing the weights $w[i,j]$ and the nodes (the $i$ and $j$ points for $\ShkMeanOne_{1}$ and $\Risky$).

Perhaps the most popular such method is Gauss-Hermite interpolation (see Judd (1998) for an exposition, or Kopecky & Suen (2010) for some alternatives). Here, we will pursue a particularly intuitive alternative: Equiprobable discretization. In this method, $m=n$ and boundaries on the joint CDF are determined in such a way as to divide up the total probability mass into submasses of equal size (each of which therefore has a mass of $n^{-2}$). This is conceptually easier if we represent the underlying shocks as statistically independent, as with $\ShkMeanOneLog_{1,t+1}$ and $\ShkMeanOneLog_{2,t+1}$ above; in that case, each submass is a square region in the $\ShkMeanOne_{1}$ and $\ShkMeanOne_{2}$ grid. We then compute the average value of $\ShkMeanOne_{1}$ and $\Risky$ conditional on their being located in each of the subdivisions of the range of the CDF. Since, in this specification, $\Risky$ is a function of $\ShkMeanOne_{1}$, the $\Risky$ values are indexed by both $i$ and $j$, but since we have written $\ShkMeanOne_{1}$ as IID, the representation of the approximating summation is even simpler than (9):

where the function $\Risky(\ShkMeanOne_{1},\ShkMeanOne_{2})$ is implicitly defined by (2).

Details can be found in the Mathematica notebook associated with this section. A particular example, in which $\sigma^{2}_{2} = \sigma^{2}_{1}$ and $\omega = 0.5$, is illustrated in Figure 1; the red dots reflect the height of the approximation to the CDF above the conditional mean values for $\ShkMeanOne_{1}$ and $\Risky$ within each of the equiprobable regions.