Managing a bond fund: Simulating the maximum of correlated normal variates

Question

Two rating agencies score the safety of bonds in a particular population on separate standard normal scales. Because the two agencies take some of the same factors into account in their ratings, the correlation between the scores of the two agencies is $\rho = 0.8.$

A bond fund manager will consider a particular bond from this population for her bond fund only if at least one of the agencies gives it a score of 1.25 or greater. According to this criterion, what proportion of the bonds is eligible for consideration?

Our initial answer to this question implements the method mentioned in another Problem on this site (see using conditional distribution to simulate joint). Although some method of simulation or numerical integration is required, this is only one possible method among many; demonstrations of alternative methods are welcome.

Answer 1

To model this situation explicitly, suppose $(X, Y)$ have a bivariate normal with zero means, unit variances, and correlation $\rho = 0.8.$ With $W = \max(X, Y)$ we wish to use simulation to evaluate $P\{W > 1.25\}$.

We can use the marginal distribution $X \sim Norm(0, 1)$ together with the conditional distribution $Y | X = x \sim Norm(\rho x, \sigma),$ where $\rho = 0.8$ and $\sigma = \sqrt{1 - \rho^2},$ to simulate the bivariate distribution. From there, we can use a simple program in R to simulate the distribution of $W,$ and approximate the desired probability.

 m = 75000;  rho = .8;  sgm = sqrt(1 - rho^2)
 x = rnorm(m);  y = rnorm(m, rho*x, sgm)
 mean(x);  sd(x);  mean(y);  sd(y);  cor(x,y)  # reality check
 ## 0.0001833525  # approximates E(X) = 0
 ## 0.9974767     # approximates SD(X) = 1
 ## 0.002128830   # approximates E(Y) = 0
 ## 0.996366      # approximates SD(Y) = 1
 ## 0.7983816     # approximates Cor(X,Y) = 0.8
 w = pmax(x,y);  mean(w > 1.25)
 ## 0.1507067     # Answer: approx P(W > 1.25)

Answers using $m = 75\,000$ iterations are accurate to about two places. A larger number of iterations would provide greater of accuracy. The current value of $m$ was chosen for maximum clarity in the plot below. Simulated values of $(X,Y)$ in the area representing $P(W > 1.25) \approx 0.15$ are plotted in dark blue.

Once we have the joint normal distribution, it is easy to estimate the proportion of the bonds that satisfy other criteria: for example (i) average score exceeds 1.25 (a little more than 9%) and (ii) maximum score exceeds 1.25 (about 6%).

 mean((x+y)/2 > 1.25)
 ## 0.09342667
 mean(pmin(x,y) > 1.25)
 ## 0.05930667