A question concerning the joint probability distribution Here is the original question:

Given a stochastic process $X(t)=Y_1+tY_2$, where $Y_1,Y_2$ are i.i.d satisfying $Y_1 \sim N(0,1)$. Derive the joint probability distribution for $(X(t),X(s))$ where $t\not = s$.

For single distribution, I know that $X(t)\sim N(0,1)+N(0,t^2)=N(0,1+t^2)$. I also know that the covariance is $\mathrm{cov}(X(t),X(s))=1+ts$. So I guess that the joint distribution is the normal distribution with two variables. But how can I prove it? The integration is too sophisticated. Is there any good way to deal with this problem?
 A: Why not use the good ol'change of variables method? To compute the joint distribution of $(X(t),X(s))$ is to be able to compute $E[A(X(t),X(s))]$ for every measurable and (say) bounded function $A$. But $(Y_1,Y_2)$ is standard normal hence
$$
E[A(X(t),X(s))]=E[A(Y_1+tY_2,Y_1+sY_2)]=\iint A(x+ty,x+sy)\mathrm e^{-(x^2+y^2)/2}\frac{\mathrm dx\mathrm dy}{2\pi}.
$$
A suitable change of variable (strongly suggested by the setting!) is $u=x+ty$, $v=x+sy$. Then $y=(u-v)/(t-s)$, $x=(tv-su)/(t-s)$, and, assuming without loss of generality that $t\gt s$, $\mathrm dx\mathrm dy=\mathrm du\mathrm dv/(t-s)$, hence
$$
E[A(X(t),X(s))]=\iint A(u,v)\mathrm e^{-B(u,v)/2}\frac{\mathrm du\mathrm dv}{2\pi(t-s)},
$$
where 
$$
B(u,v)=\frac{(tv-su)^2+(u-v)^2}{(t-s)^2}.
$$
This proves that for every $t\gt s$, the distribution of $(X(t),X(s))$ has density $f_{t,s}$ where
$$
f_{t,s}(u,v)=\frac{\mathrm e^{-B(u,v)/2}}{2\pi(t-s)}.
$$
Developing the $2$-homogenous polynomial $B(u,v)$ yields what you seek, namely,
$$
B(u,v)=\frac{(s^2+1)u^2-2(st+1)uv+(t^2+1)v^2}{(t-s)^2}.
$$
A: With some calculations, I have actually reached a solution. However, I'm not quite sure it is a true one. 
To calculate the density function, we consider 
\begin{align}
f(x_t,x_s)&=P(X(t)=x_t,X(s)=x_s)\\&=P(Y_1+tY_2=x_t,Y_1+sY_2=x_s)\\&=P\left(Y_1=\frac{tx_s-sx_t}{t-s},Y_2=\frac{x_t-x_s}{t-s}\right)\\
&=\frac{\partial(Y_1,Y_2)}{\partial(X_1,X_2)}\frac{1}{\sqrt{2\pi}}e^{-\frac{\left(\frac{tx_s-sx_t}{t-s}\right)^2}{2}}\frac{1}{\sqrt{2\pi}}e^{-\frac{\left(\frac{x_t-x_s}{t-s}\right)^2}{2}}\\
&=\frac{1}{2\pi|t-s|}e^{-\frac{(1+s^2)x_t-2(1+ts)x_tx_s+(1+t^2)x_s}{2(t-s)^2}}\\
&=\frac{1}{2\pi \sigma_1\sigma_2 \sqrt{1-\rho^2}} e^{-\frac{\frac{x_t^2}{\sigma_1^2}-\frac{2\rho x_t x_s}{\sigma_1 \sigma_2}+\frac{x_s^2}{\sigma_2^2}}{2(1-\rho^2)}}
\end{align}
where $\sigma_1=1+t^2,\sigma_2=1+s^2,\rho=\frac{1+ts}{\sqrt{1+t^2}\sqrt{1+s^2}}.$
Thus $(X_t,X_s)\sim N(0,0,\sigma_1,\sigma_2,\rho).$
