Suppose $x\sim N(\mu_x,\sigma_x^2)$ and we are given that $y\mid x \sim N(a+bx,\sigma^2)$, where $a$ and $b$ are some constants.
It is a fact that the joint distribution of $x$ and $y$ is a bivariate normal. How can one explicitly calculate the bivariate normal parameters from first principles?
In particular I wonder if there is a "clever way", since tedious calculation can produce a messy answer to this question after multiplying the density of $x$ with the conditional density of $y\mid x$.
The normal bivariate distribution is characterized by the parameters $$E(X),\qquad E(Y),\qquad\mathrm{var}(X),\qquad\mathrm{cov}(X,Y),\qquad\mathrm{var}(Y).$$ To compute these, note that $X=\mu_X+\sigma_XU$ and $Y=a+bX+\sigma V$ where $(U,V)$ is bivariate standard normal, that is, such that. $$E(U)=E(V)=\mathrm{cov}(U,V)=0,\qquad\mathrm{var}(U)=\mathrm{var}(V)=1.$$