That's a great question and it is what I always wanted to try to do.
I guess I found a solution using PDE approach. Change of numeraire would be more intuitive indeed, but I am not very good in stochastic calculus.
The idea is as follows:
1) Let's consider portfolio $\Pi = V(X,Y,t) - \Delta_X X - \Delta_Y Y$.
I will found $\Delta_X$ and $\Delta_Y$ such that portfolio $\Pi$ would be riskless and earn risk-free rate of return $r$: $d\Pi = r\Pi dt$.
Assumption: $dX = \mu_X X dt + \sigma_X X dW^X$, $dY = \mu_Y Y dt + \sigma_Y Y dW^Y$ and $dW^X dW^Y = \rho dt$.
Hence, applying Ito's lemma I obtain:
$d\Pi = \frac{\partial V}{\partial t} dt + \frac{\partial V}{\partial X} dX + \frac{\partial V}{\partial Y} dY + \frac{1}{2} \sigma_X^2 X^2 \frac{\partial^2 V}{\partial X^2} dt+ \frac{1}{2} \sigma_Y^2 Y^2 \frac{\partial^2 V}{\partial Y^2} dt+ \rho \sigma_X\sigma_Y XY \frac{\partial^2 V}{\partial X\partial Y} dt - \Delta_X dX - \Delta_Y dY =$
$\left( \frac{\partial V}{\partial t} + \frac{1}{2} \sigma_X^2 X^2 \frac{\partial^2 V}{\partial X^2}+ \frac{1}{2} \sigma_Y^2 Y^2 \frac{\partial^2 V}{\partial Y^2} + \rho \sigma_X\sigma_Y XY \frac{\partial^2 V}{\partial X\partial Y} \right)dt + \left(\frac{\partial V}{\partial X} - \Delta_X \right) dX + \left(\frac{\partial V}{\partial Y} - \Delta_Y \right) dY$.
And all this is equal to $d\Pi = r\Pi dt = r\left(V - \Delta_X X - \Delta_Y Y\right)dt$
Now, set $\frac{\partial V}{\partial Y} = \Delta_Y$ and $\frac{\partial V}{\partial X} = \Delta_X$.
Left-hand side becomes $\left(\frac{\partial V}{\partial t} + \frac{1}{2} \sigma_X^2 X^2 \frac{\partial^2 V}{\partial X^2}+ \frac{1}{2} \sigma_Y^2 Y^2 \frac{\partial^2 V}{\partial Y^2} + \rho \sigma_X\sigma_Y XY \frac{\partial^2 V}{\partial X\partial Y}\right) dt$
Right-hand side is now $r\left(V - \frac{\partial V}{\partial X} X - \frac{\partial V}{\partial Y} Y\right)dt$
The PDE is now $\frac{\partial V}{\partial t} + \frac{1}{2} \sigma_X^2 X^2 \frac{\partial^2 V}{\partial X^2}+ \frac{1}{2} \sigma_Y^2 Y^2 \frac{\partial^2 V}{\partial Y^2} + \rho \sigma_X\sigma_Y XY \frac{\partial^2 V}{\partial X\partial Y} = r\left(V - \frac{\partial V}{\partial X} X - \frac{\partial V}{\partial Y} Y\right)$, or
$\frac{\partial V}{\partial t} + \frac{1}{2} \sigma_X^2 X^2 \frac{\partial^2 V}{\partial X^2}+ \frac{1}{2} \sigma_Y^2 Y^2 \frac{\partial^2 V}{\partial Y^2} + \rho \sigma_X\sigma_Y XY \frac{\partial^2 V}{\partial X\partial Y} + r\frac{\partial V}{\partial X} X + r \frac{\partial V}{\partial Y} Y = rV$
I forgot: the boundary condition is $V(X, Y, T) = (X - Y)^+$
2) Now, in order to solve this crazy PDE i will use substitution:
$Z = \frac{X}{Y}$ and $V(X,Y,t) = G(Z, t)$.
Thanks to Wolfram Alpha, I have:
$\frac{\partial V}{\partial X} = \frac{1}{Y} \frac{\partial G}{\partial Z}$
$\frac{\partial V}{\partial Y} = -\frac{X}{Y} \frac{\partial G}{\partial Z}$
$\frac{\partial^2 V}{\partial X^2} = -\frac{1}{Y^2} \frac{\partial^2 G}{\partial Z^2}$
$\frac{\partial^2 V}{\partial Y^2} = \frac{X\left(2Y\frac{\partial G}{\partial Z}+X\frac{\partial^2 G}{\partial Z^2}\right)}{Y^4} $
$\frac{\partial^2 V}{\partial XY} = -\frac{Y\frac{\partial G}{\partial Z}+X\frac{\partial^2 G}{\partial Z^2}}{Y^3} $
Substituting into previous equation and cancelling the terms out we obtain:
$\dot{G} + [\sigma_X^2-\rho \sigma_X \sigma_Y]ZG' + \frac{1}{2}[\sigma_X^2-2\rho \sigma_X \sigma_Y + \sigma_Y^2]Z^2G'' = rG$, or
$\dot{G} + \mu_GZG' + \frac{1}{2}\sigma_G^2 Z^2G'' = rG$, where
$\dot{G} = \frac{dG}{dt}$, $G' = \frac{dG}{dZ}$
$\mu_G = \sigma_X^2-\rho \sigma_X \sigma_Y$, $\sigma_G = \sqrt{\sigma_X^2-2\rho \sigma_X \sigma_Y + \sigma_Y^2}$
And the boundary condition is $G(Z,T) = Y(Z - 1)^+$
UPDATE: PREVIOUS VERSION WAS NOT COMPLETELY CORRECT
3) Now the question is what to do with that $Y$ in the equation above? I employ next change of variables: $G(Z) = YF(Z)$.
Thanks to paper and pencil, I have:
$G' = (YF)' = Y\left(F' - \frac{F}{Z}\right)$ and $G'' = \left((YF)'\right)' = \text{after some calculations} = YF''$
Plugging this into $Z$'s PDE we obtain:
$\dot{F} + \mu_G Z F' + \frac{1}{2} \sigma_G^2 Z^2 F'' = (r+ \mu_G)F$
with boundary condition $F(Z,T) = (Z-1)^+$
Now denote $r^* = r+ \mu_G$ and equation becomes: $\dot{F} + (r^* - r) Z F' + \frac{1}{2} \sigma_G^2 Z^2 F'' = r^* F$
4) Now $r^*$ works like new risk-free rate and $r$ is like $Z$'s dividend yield and we can apply well-known formula for option on asset with continiously paid dividends:
$F(Z, T) = e^{-r^*T} N(d_1) Z_0 - e^{-rT} N(d_2) $, where
$d_{1,2} = \frac{1}{\sigma_G\sqrt{T}}\left[\ln\left(Z_0\right)+\left(r^* - r \pm\frac{\sigma_G^2}{2}\right)T\right] =\frac{1}{\sigma_G\sqrt{T}}\left[\ln\left(Z_0\right)+\left(\mu_G \pm\frac{\sigma_G^2}{2}\right)T\right] $.
5) Now $V = e^{-r^*T} N(d_1) X_0 - e^{-rT} N(d_2) Y_0$, where
$d_{1,2} =\frac{1}{\sigma_G\sqrt{T}}\left[\ln\left(\frac{X_0}{Y_0}\right)+\left(\mu_G \pm\frac{\sigma_G^2}{2}\right)T\right] $, where
$r^* = r+ \mu_G$
$\mu_G = \sigma_X^2-\rho \sigma_X \sigma_Y$
$\sigma_G = \sqrt{\sigma_X^2-2\rho \sigma_X \sigma_Y + \sigma_Y^2}$
Hope I was correct.
I also hope somebody would be able to propose any better solution, maybe using martingale approach.