Conditional expectation of two identical marginal normal random variables 
Let $Y_0$ and $Y_1$ be both identically (not necessarily independent)
   normally distributed with mean $\mu$ and $\sigma^2$, i.e., $Y_i \sim
 N(\mu, \sigma^2)$ for $i = 1, 2$. Let $\rho$ denote the correlation
   between $Y_0$ and $Y_1$. Show that the conditional expectation
   $E[Y_1|Y_0 = y_0] = (1-\rho)\mu + \rho y_0$.

I know that in general if $(X, Y)$ are bivariate normally distributed, then $E[X|Y] = E[X] + \rho \frac{\sigma_X}{\sigma_Y}(Y-E[Y])$. But in this question, there is no indication that $Y_0$ and $Y_1$ are jointly normally distributed; they are simply marginally normally distributed. If they were joint normal, the result falls out simply from the result I just mentioned. Any ideas on how I can prove the result?
 A: I don't think the formula given in the question can be correct in all cases, it is developed using joint normality.  Without joint normality we can use copulas.  For $X,Y$ random variables with joint distribution with cumulative distribution function $F(x,y)$ and joint density $f(x,y)$ define the transformed random variables (rv) $U=F_X(X), V=F_Y(Y)$ where $F_X, F_Y$ denotes the marginal cumulative distribution functions (cdf). Then the joint distribution of $U;V$
$$ \DeclareMathOperator{\P}{\mathbb{P}}
  \P(U \le u, V \le v)=C(u,v)
$$
is the copula of $X$ and $Y$, with copula density $c(u,v)$ (when it exists).  So, in your setup let us assume that the copula density exists, and for simplicity I will take both $X$ and $Y$ as standard normals.  So what possibility exists for the conditional expectation of $X$ given $Y=y$ ? Using Sklar's theorem we can write the joint density as
$$
   f(x,y) = c(\Phi(x),\Phi(y)) \phi(x) \phi(y)
$$
where $\Phi, \phi$ are the standard normal cdf, pdf, respectively.  Then the conditional density is given by
$$
   f(x \mid y) = \frac{f(x,y)}{f(y)}=\frac{c(\Phi(x),\Phi(y))\phi(x)\phi(y)}{\phi(y)}= c(\Phi(x),\Phi(y))\phi(x)
$$
Then we can look at this with various copula functions, see https://en.wikipedia.org/wiki/Copula_(probability_theory) .
There is a general inequality for copulas
$$
  W(u,v)=\max(u+v-1,0) \le C(u,v) \le M(u,w)=\min(u,v)
$$
where both upper and lower limits are copulas (This is the Frechet-Hoeffding bounds). The upper limit isn't very interesting, since it gives $\P(U=V)=1$ so gives correlation equal 1. The lower limit similarly corresponds to $U=1-V$ with probability one, so correlation is -1.  But for these two extremal copula the conditional expectation function certainly is linear!
Lets look at some intermediate cases.  I will use the R package copula and some numerical integration to find the conditional expectation function, for the case of the gumbel copula. The code can be simply adapted for other copulas.
    library(copula)
    C  <-  gumbelCopula(2)
    make_cond <-  Vectorize(function(y) {
        function(x) dCopula( cbind(pnorm(x), pnorm(y)), C)*dnorm(x) 
    }  )

the last command makes a function representing the conditional density given $Y=y$. Let us look at how this looks like for three different values of $y$:
    cond_dens <- make_cond(c(-1, 0, 1))
    plot(cond_dens[[1]],  from=-3,  to=3, col="blue", ylim=c(0, 0.7))
    plot(cond_dens[[2]],  from=-3,  to=3,  col="orange",  add=TRUE)
    plot(cond_dens[[3]],  from=-3,  to=3,  col="red",  add=TRUE)
    title("conditional densities for y=-1, 0, 1")


showing clearly that the conditional distributions now are non-normal. We can also see clearly that the conditional variance is non-constant.
For more examples using the copula package see https://www.r-bloggers.com/modelling-dependence-with-copulas-in-r/ and Generating values from copula using copula package in R . Then we can find the conditional expectation function using numerical integration:
    plot(function(y) sapply(make_cond(y), FUN=function(fun)
                   integrate(function(x) x*fun(x) ,
                             lower=-Inf,  upper=Inf)$value), 
                              from=-3,  to=3,
         ylab="conditional expectation given y", xlab="y")
    title("conditional expectation of X given Y=y")


and it is quite clear that the conditional expectation function is not linear!
A: Here is an example using the Ali-Mikhail-Haq copula: with cdf
$$F_\theta(x,y)=\dfrac{xy}{1-\theta(1-x)(1-y)}\qquad\theta\in(-1,1)$$
and associated joint density
$$f_\theta(x,y)=\dfrac{1+\theta[(1+x)(1+y)-3]+\theta^2(1-x)(1-y))
}{[1-\theta(1-x)(1-y)]^3}$$
Simulating a sample from the conditional distribution associated with this joint distribution (by Metropolis-Hastings) and applying the Normal quantile function transform to the sample returns the following graph of $\mathbb{E}[X|Y]$ for $\theta=0.5, -0.9$, clearly not a linear regression:


A: Consider $X$ and $Y$ with joint pdf
$$f_{X,Y}(x,y) = \begin{cases}2\phi(x)\phi(y), & x, y 
\geq 0 ~~ \text{or}~~ x, y < 0,\\ 0, & \text{otherwise.}\end{cases}$$
It is readily verified that $X$ and $Y$ are marginally standard normal random variables, but they are not jointly normal random variables, and they are not independent random variables; their correlation coefficient $\rho$ is positive and the linear minimum mean-square error (MMSE) estimate of $Y$ given that $X=x$ is just $\rho x$, a straight line with slope $\rho$ passing through the origin. On the other hand, for $x \geq 0$, the 
conditional density of $Y$ given that $X=x$ is 
$$f_{Y\mid X}(y\mid X=x) = \frac{f_{X,Y}(x,y)}{f_X(x)} = 2\phi(y)\mathbf 1_{y\geq 0}$$ and so
\begin{align}
E[Y\mid X = x] &= \int_0^\infty y \cdot 2\phi(y) \,\mathrm dy = \sqrt{\frac{2}{\pi}}.
\end{align}
Similarly, for $x<0$, $E[Y \mid X=x] = -\sqrt{\frac{2}{\pi}}$, that is,
regarded as a function of $x$, $E[Y \mid X=x]$ is a piecewise constant function (or step function) that jumps abruptly at $x=0$ from value $-\sqrt{\frac{2}{\pi}}$ on the negative axis to value $+\sqrt{\frac{2}{\pi}}$ on the positive axis.  This is quite different from the linear MMSE estimate of $\rho x$. Indeed, the two estimates are equal only at two points in the first and third quadrants respectively.
In short, I don't believe that the result that the OP wishes to prove is correct. Of course, as the OP says, if $X$ and $Y$ are jointly normal, then the result is equivalent to the fact that for jointly normal random variables, the linear MMSE estimator is the same as the MMSE estimator $E[Y\mid X=x]$.
