CDF of maximum of $n$ correlated normal random variables The maximum of $n$ normal i.i.d. random variables
$$Y=\max\{x_1,...,x_n\},$$
$$x_i \sim N[0,1]$$
has the CDF
$$P(Y\le y)=\Phi(y)^n  $$
but how does the CDF look like, if the variables are identically correlated?
$$\text{corr}(x_i, x_j)=\rho,$$
$$0\le\rho\le1.$$
I am looking for the analytical formula
$$P(Y\le y)=f(y,n,\rho).$$

This answer is essentially the same method as the answer by Thomas Lumley, but I'll give you a bit more detail on the derivations.  Suppose we take independent values:
$$\tilde{X}_i \sim \text{N}(0, 1-\rho^2)
\quad \quad \quad \quad \quad 
Z \sim \text{N}(0, \rho^2).$$
Setting $X_i = \tilde{X}_i+Z$ gives a random vector containing equicorrelated values with unit variance:
$$\mathbf{X} \sim \text{N}(\mathbf{0}, \mathbf{\Sigma})
\quad \quad \quad \quad \quad 
\mathbf{\Sigma} = \begin{bmatrix}
1 & \rho & \cdots & \rho \\
\rho & 1 & \cdots & \rho \\
\vdots & \vdots & \ddots & \vdots \\
\rho & \rho & \cdots & 1 \\
\end{bmatrix}.$$
Let $\phi$ and $\Phi$ denote the PDF and CDF of the standard normal distribution.  We can then write the CDF of the random variable of interest as:
$$\begin{align}
F_Y(y)
&\equiv \mathbb{P}(Y \leqslant y) \\[16pt]
&= \mathbb{P}(\max (X_1,...,X_n) \leqslant y) \\[16pt]
&= \mathbb{P}(Z + \max (\tilde{X}_1,...,\tilde{X}_n) \leqslant y) \\[12pt]
&= \int \limits_{-\infty}^\infty \mathbb{P}(Z + \max (\tilde{X}_1,...,\tilde{X}_n) \leqslant y |Z=z) \cdot \text{N}(z|0,\rho^2) \ dz \\[6pt]
&= \int \limits_{-\infty}^\infty \mathbb{P}(\max (\tilde{X}_1,...,\tilde{X}_n) \leqslant y-z) \cdot \text{N}(z|0,\rho^2) \ dz \\[6pt]
&= \int \limits_{-\infty}^\infty \mathbb{P}\Bigg( \frac{\max (\tilde{X}_1,...,\tilde{X}_n)}{\sqrt{1-\rho^2}} \leqslant \frac{y-z}{\sqrt{1-\rho^2}} \Bigg) \cdot \text{N}(z|0,\rho^2) \ dz \\[6pt]
&= \int \limits_{-\infty}^\infty \mathbb{P} \Bigg( \frac{\tilde{X}_1}{\sqrt{1-\rho^2}} \leqslant \frac{y-z}{\sqrt{1-\rho^2}}, ..., \frac{\tilde{X}_n}{\sqrt{1-\rho^2}} \leqslant \frac{y-z}{\sqrt{1-\rho^2}} \Bigg) \cdot \text{N}(z|0,\rho^2) \ dz \\[6pt]
&= \int \limits_{-\infty}^\infty \Phi \bigg( \frac{y-z}{\sqrt{1-\rho^2}} \bigg)^n \cdot \text{N}(z|0,\rho^2) \ dz \\[6pt]
&= \frac{1}{|\rho|} \int \limits_{-\infty}^\infty \Phi \bigg( \frac{y-z}{\sqrt{1-\rho^2}} \bigg)^n \cdot \phi \bigg( \frac{z}{\rho} \bigg) \ dz. \\[6pt]
\end{align}$$
The corresponding density is:
$$\begin{align}
f_Y(y)
&= \frac{dF_Y}{dy}(y) \\[16pt]
&= \frac{1}{|\rho|} \frac{d}{dy} \int \limits_{-\infty}^\infty \Phi \bigg( \frac{y-z}{\sqrt{1-\rho^2}} \bigg)^n \cdot \phi \bigg( \frac{z}{\rho} \bigg) \ dz \\[6pt]
&= \frac{1}{|\rho|} \int \limits_{-\infty}^\infty \frac{\partial}{\partial y} \Phi \bigg( \frac{y-z}{\sqrt{1-\rho^2}} \bigg)^n \cdot \phi \bigg( \frac{z}{\rho} \bigg) \ dz \\[6pt]
&= \frac{n}{|\rho| \sqrt{1-\rho^2}} \int \limits_{-\infty}^\infty \Phi \bigg( \frac{y-z}{\sqrt{1-\rho^2}} \bigg)^{n-1} \cdot \phi \bigg( \frac{y-z}{\sqrt{1-\rho^2}} \bigg) \cdot \phi \bigg( \frac{z}{\rho} \bigg) \ dz. \\[6pt]
\end{align}$$
This function has no closed form expression (except in the trivial case where $n=1$) but it can be evaluated using numerical methods.
