Consider a market composed by two stocks whose prices $X$ and $Y$ are given by B&S diffusion:
$$dX_t= \mu X_t dt+ \sigma X_tdW_t$$
$$dY_t= \mu Y_t dt+ \sigma Y_tdB_t$$
Supposing the market is complete, how to evaluate the fair price of an option whose payoff is $$ \phi(X_T,Y_T)=(X_T-Y_T)_+$$
My idea was to apply a change of numeraire technique and so obtain price as a function of the B&S formula. However, I was not able to find it.
Any advice would be appreciated.
Measure change is still the most natural approach for such problems. We assume that, under the measure $P$, \begin{align*} dX_t &= \mu X_t dt + \sigma X_t dW_t^1,\\ dY_t &= \mu Y_t dt + \sigma Y_t \left(\rho dW_t^1 + \sqrt{1-\rho^2} dW_t^2 \right), \end{align*} based on the Cholesky decomposition, where $\{W_t^1, t \ge 0\}$ and $\{W_t^2, t \ge 0\}$ are two independent standard Brownian motions, and $\rho$ ($|\rho|<1$) is the correlation. We define the measure $Q$ such that \begin{align*} \frac{dQ}{dP}\big|_t = \exp\left(-\frac{1}{2}\sigma^2 t + \sigma\left(\rho W_t^1 + \sqrt{1-\rho^2} W_t^2 \right) \right). \end{align*} Then, $\{\widehat{W}_t^1, t \ge 0\}$ and $\{\widehat{W}_t^2, t \ge 0\}$ are two independent standard Brownian motions under $Q$, where \begin{align*} \widehat{W}_t^1 &= W_t^1 - \rho\sigma t,\\ \widehat{W}_t^2 &= W_t^2 - \sqrt{1-\rho^2}\sigma t. \end{align*} Moreover, \begin{align*} dX_t &= (\mu +\rho\sigma^2) X_t dt + \sigma X_t d\widehat{W}_t^1,\\ dY_t &= (\mu + \sigma^2) Y_t dt + \sigma Y_t \left(\rho d\widehat{W}_t^1 + \sqrt{1-\rho^2} d\widehat{W}_t^2 \right), \end{align*} and then \begin{align*} Y_T &= Y_0 \exp\left(\Big(\mu + \frac{1}{2}\sigma^2 \Big)T + \sigma\left(\rho \widehat{W}_T^1 + \sqrt{1-\rho^2} \widehat{W}_T^2 \right) \right),\\ \frac{X_T}{Y_T} &= \frac{X_0}{Y_0}\exp\left((\rho-1)\sigma^2T +\sigma(1-\rho)\widehat{W}_T^1-\sigma\sqrt{1-\rho^2} \widehat{W}_T^2 \right)\\ &=\frac{X_0}{Y_0}\exp\left((\rho-1)\sigma^2T +\sqrt{2(1-\rho)}\sigma \frac{\sigma(1-\rho)\widehat{W}_T^1-\sigma\sqrt{1-\rho^2} \widehat{W}_T^2}{\sqrt{2(1-\rho)}\sigma} \right)\\ &=\frac{X_0}{Y_0}\exp\left(-\frac{1}{2}\hat{\sigma}^2T +\hat{\sigma} W_T \right), \end{align*} where $\hat{\sigma} = \sqrt{2(1-\rho)}\sigma$, and \begin{align*} W_t =\frac{\sigma(1-\rho)\widehat{W}_t^1-\sigma\sqrt{1-\rho^2} \widehat{W}_t^2}{\sqrt{2(1-\rho)}\sigma} \end{align*} is a standard Brownina motion, by Levy's characterization. Therefore, \begin{align*} E_P\left( (X_T-Y_T)^+\right) &= E_P\left(Y_T \left(\frac{X_T}{Y_T}-1\right)^+\right)\\ &=E_Q\left( \left( \frac{dQ}{dP}\big|_T\right)^{-1}Y_T \left(\frac{X_T}{Y_T}-1\right)^+\right)\\ &=Y_0e^{(u+\sigma^2)T}E_Q\left(\left(\frac{X_T}{Y_T}-1\right)^+ \right)\\ &= Y_0e^{(u+\sigma^2)T} \left[\frac{X_0}{Y_0}\Phi(d_1) - \Phi(d_2) \right]\\ &= e^{(u+\sigma^2)T} \Big[X_0\Phi(d_1) - Y_0 \Phi(d_2) \Big]. \end{align*} where \begin{align*} d_1 &= \frac{\ln \frac{X_0}{Y_0} + \frac{1}{2}\hat{\sigma}^2T}{\hat{\sigma}\sqrt{T}}\\ &=\frac{\ln \frac{X_0}{Y_0} + (1-\rho)\sigma^2T}{\sqrt{2(1-\rho)}\sigma\sqrt{T}},\\ d_2 &= d_1 - \hat{\sigma}\sqrt{T}\\ &= d_1 -\sqrt{2(1-\rho)}\sigma\sqrt{T}. \end{align*}
I think there are 2 ways to get the answer. First way is what Gordon said. But when I first saw his answer, I didn't know why he defined Radon–Nikodym like that, so I thought about it for a long time, trying to give my understanding here which was inspired by this answer.
We want to define a new measure $\mathbb{Q}^{1} \sim \mathbb{Q}$ which uses the $Y_{T}$ to be the numeraire. Then we define,
$$ \frac{\mathrm{d} \mathbb{Q}^1}{\mathrm{~d} \mathbb{Q}}=\frac{Y_T}{Y_0} \frac{B_0}{B_T}=\frac{Y_T}{Y_0} e^{-r T} . $$
The reason why we define this Radon–Nikodym are
2.$\frac{\mathrm{d} \mathbb{Q}^1}{\mathrm{~d} \mathbb{Q}} > 0$ under BS model can make sure $\mathbb{Q}^1 \sim \mathbb{Q}$ and use ratio of 2 numeraire can make sure $\mathbb{Q}^1$ an equivalent martingale measure.($Y_{T}$/numeraire is a martingale under corresponding measure)
$$ \frac{Y_{0}}{B_{0}} = E^{\mathbb{Q}}[\frac{Y_{T}}{B_{T}}] = E^{\mathbb{Q}^1}[\frac{\mathrm{~d} \mathbb{Q}}{\mathrm{d} \mathbb{Q}^1}\frac{Y_{T}}{B_{T}}] = E^{\mathbb{Q}^1}[\frac{Y_{T}}{Y_{T}}\frac{Y_{0}}{B_{0}}] $$ Thus we have $$ E^{\mathbb{Q}^1}[\frac{Y_{T}}{Y_{T}}] = \frac{Y_{0}}{Y_{0}} $$
Back to the question, under measure $\mathbb{Q}$, $$ \begin{gathered} d X_t=r X_t d t+\sigma X_t d W^{1}_t \\ d Y_t= r Y_t d t+\sigma Y_t d W^{2}_t \end{gathered} $$ where $W^{1}_{t}$ and $W^{2}_t$ are 2 Brownian Motion under $\mathbb{Q}$ and $dW^{1}_{t}dW^{2}_t = \rho t$.
We know that $Y_{T} = Y_{0}e^{(r - \frac{1}{2}\sigma^2)T - \sigma W^{2}_{T}}$. If we just substitude this $Y_{T}$ formula into $\frac{\mathrm{d} \mathbb{Q}^1}{\mathrm{~d} \mathbb{Q}}$, this Radon–Nikodym only contatins one Brownian Motion $W^{2}_{T}$, we can not use multi-dimension Girsanov. Thus here is the trick, we use cholesky decompostion to obtain that $$ \begin{aligned} d X_t&=r X_t d t+\sigma X_t d W^{'}_t \\ d Y_t&= r Y_t d t+\sigma Y_t d (\rho W^{'}_t + \sqrt{1 - \rho^2} W^{''}_t) \end{aligned} $$ where $W^{'}_t$ and $W^{''}_t$ are 2 independent Brownian Motion.Thus we have $Y_{T} = Y_{0}e^{(r - \frac{1}{2}\sigma^2)T - \sigma(\rho W^{'}_T + \sqrt{1 - \rho^2} W^{''}_T)}$. Now we substitude $Y_{T}$ into $\frac{\mathrm{d} \mathbb{Q}^1}{\mathrm{~d} \mathbb{Q}}$, then $$ \frac{\mathrm{d} \mathbb{Q}^1}{\mathrm{~d} \mathbb{Q}}=\frac{Y_T}{Y_0} e^{-r T} = \frac{Y_{0}e^{(r - \frac{1}{2}\sigma^2)T - \sigma(\rho W^{'}_T + \sqrt{1 - \rho^2} W^{''}_T)}}{Y_0} e^{-r T} = e^{-\frac{1}{2}\sigma^2T - \sigma(\rho W^{'}_T + \sqrt{1 - \rho^2} W^{''}_T)} $$ which is coninside with how Gordon defined.
The second way is we investigate the dynamic of $\frac{X_{T}}{Y_{T}}$. Use ito formula to $f:= \frac{x}{y}$ we can derive that $$ \begin{aligned} df &= \left(\sigma^2-\rho \sigma^2\right) f d t-\sqrt{2\sigma^2-2 \rho \sigma^2} f d W_{t} \end{aligned} $$ Thus to make $\frac{X_{T}}{Y_{T}}$ a martingale, we need to have $\bar{W}_{t}:= W_{t} - \frac{\sigma^2-\rho \sigma^2}{\sqrt{2\sigma^2-2 \rho \sigma^2} } t$. Thus by girsanov, we define the similar Radon–Nikodym.(Actually I didn't use cholesky in this method, If we use cholesky first and do the ito formula, I think it will result in the same Radon-Nikodym by the reason maybe the uniqueness of Martingale Measure.)
By the way, I think in Gordon answer of $$\begin{aligned} E_P\left(\left(X_T-Y_T\right)^{+}\right) & =E_P\left(Y_T\left(\frac{X_T}{Y_T}-1\right)^{+}\right) \\ & =E_Q\left(\left(\left.\frac{d Q}{d P}\right|_T\right)^{-1} Y_T\left(\frac{X_T}{Y_T}-1\right)^{+}\right) \\ & =Y_0 e^{\left(u+\sigma^2\right) T} E_Q\left(\left(\frac{X_T}{Y_T}-1\right)^{+}\right) \\ & =Y_0 e^{\left(u+\sigma^2\right) T}\left[\frac{X_0}{Y_0} \Phi\left(d_1\right)-\Phi\left(d_2\right)\right] \\ & =e^{\left(u+\sigma^2\right) T}\left[X_0 \Phi\left(d_1\right)-Y_0 \Phi\left(d_2\right)\right] . \end{aligned}$$ should be $$ \begin{aligned} E_P\left(\left(X_T-Y_T\right)^{+}\right) & =E_P\left(Y_T\left(\frac{X_T}{Y_T}-1\right)^{+}\right) \\ & =E_Q\left(\left(\left.\frac{d Q}{d P}\right|_T\right)^{-1} Y_T\left(\frac{X_T}{Y_T}-1\right)^{+}\right) \\ & =Y_0 e^{\mu T} E_Q\left(\left(\frac{X_T}{Y_T}-1\right)^{+}\right) \\ & =Y_0 e^{\mu T}\left[\frac{X_0}{Y_0} \Phi\left(d_1\right)-\Phi\left(d_2\right)\right] \\ & =e^{\mu T}\left[X_0 \Phi\left(d_1\right)-Y_0 \Phi\left(d_2\right)\right] . \end{aligned} $$
Relatively quick Solution
If $U$ and $V$ be normally distributed with means $\mu_u\,,\,\mu_v$, variances $\sigma^2_u\,,\,\sigma^2_v$ and correlation $\rho$ then we can show ( by definition of expectation and apply joint density function ) $$\mathbb{E}\left[\left(e^U-e^V\right)^+\right]={\large{e^{\mu_u+\frac{1}{2}\sigma_u^2}}}\Phi\left(d_1\right)-{\large{e^{\mu_v+\frac{1}{2}\sigma_v^2}}}\Phi\left(d_2\right)\qquad (1)$$ where $$d_1=\frac{\mu_u-\mu_v+\sigma_u^2-\rho\sigma_u\sigma_v}{\sqrt{\sigma_u^2-2\rho\sigma_u\sigma_v+\sigma_v^2}}$$ and $$d_2=\frac{\mu_v-\mu_u+\sigma_v^2-\rho\sigma_v\sigma_u}{\sqrt{\sigma_v^2-2\rho\sigma_v\sigma_u+\sigma_u^2}}$$ By application of Ito's lemma and Girsanov theorem we have $$\ln X_T=\ln X_t+\left(r-\frac{1}{2}{\sigma_x}^2\right)(T-t)+{\sigma_x}({W_T}^{\mathbb{Q}}-{W_t}^{\mathbb{Q}})$$ $$\ln Y_T=\ln Y_t+\left(r-\frac{1}{2}{\sigma_y}^2\right)(T-t)+{\sigma_y}({B_T}^{\mathbb{Q}}-{B_t}^{\mathbb{Q}})$$ let $U=\ln X_T$ and $V=\ln Y_T$ apply $(1)$
Edit for SRKX
We know $$ {{f}_{U,V}}(u,v)=\frac{1}{2\pi {{\sigma }_{u}}{{\sigma }_{v}}\sqrt{1-{{\rho }^{2}}}}{{e}^{\large{-\frac{1}{2(1-{{\rho }^{2}})}\left[ {{\left( \frac{u-{{\mu }_{u}}}{{{\sigma }_{u}}} \right)}^{2}}-2\rho \left( \frac{u-{{\mu }_{u}}}{{{\sigma }_{u}}} \right)\left( \frac{v-{{\mu }_{v}}}{{{\sigma }_{v}}} \right)+{{\left( \frac{v-{{\mu }_{v}}}{{{\sigma }_{v}}} \right)}^{2}} \right]}}} $$ let $$I=\int_{-\infty }^{+\infty }{\int_{v }^{\infty }{{{e}^{u}}}}{{f}_{U,\,V}}(u,v)dudv$$ and define $$g(u,v)=\frac{1}{{{\sigma }_{u}}\sqrt{2\pi (1-{{\rho }^{2}})}}{{e}^{\large{-\frac{1}{2(1-{{\rho }^{2}})}{{\left[ u-\left( {{\mu }_{u}}+\rho \sigma_u \left( \frac{v-{{\mu }_{v}}}{{{\sigma }_{v}}} \right) \right) \right]}^{2}}}}}$$ Indeed $g$ is probability density function normal distribution. we have $$I=\int_{-\infty }^{\infty }{\frac{1}{{{\sigma }_{v}}\sqrt{2\pi }}{{e}^{\large{-\frac{1}{2}{{\left( \frac{v-{{\mu }_{v}}}{{{\sigma }_{v}}} \right)}^{2}}}}}}\left[ \int_{v }^{\infty }{{{e}^{u}}}g(u,v)du \right]dv$$ thus we can write $$I={{e}^{{{\mu }_{u}}+\frac{1}{2}(1-{{\rho }^{2}}){{\sigma }_{u}}^{2}}} \int_{-\infty }^{\infty }{\frac{1}{{{\sigma }_{v}}\sqrt{2\pi }}{{e}^{\large{-\frac{1}{2}{{\left[ \left(\frac{v-{{\mu }_{v}}}{{{\sigma }_{v}}}\right)^2-2\rho\sigma_v\left(\frac{v-\mu_v}{\sigma_v}\right) \right]}}}}}}\Phi \left[ \Lambda \right]dv$$ where $$\Lambda =\frac{{{\mu }_{v}}+\rho {{\sigma }_{u}}\left( \frac{v-{{\mu }_{v}}}{{{\sigma }_{v}}} \right)+{{(1-\rho )}^{2}}{{\sigma }_{u}}^{2}-v}{{{\sigma }_{u}}\sqrt{1-{{\rho }^{2}}}}$$ let $y=\large \frac{v-{{\mu }_{v}}-\rho {{\sigma }_{u}}{{\sigma }_{v}}}{{{\sigma }_{v}}}$ by Change of variable and setting good variable we have $$I={{e}^{{{\mu }_{u}}+\frac{1}{2}\sigma _{u}^{2}}}\Phi ({{d}_{1}})$$ using similar steps for the case $$J=\int_{-\infty }^{+\infty }{\int_{v }^{\infty }{{{e}^{v}}}}{{f}_{U,\,V}}(u,v)dudv$$ as a result
$$J={{e}^{{{\mu }_{V}}+\frac{1}{2}\sigma _{V}^{2}}}\Phi ({{d}_{2}})$$ Now we should use $I$ and $J$ and calculate $$E\left[(e^U-e^V)^+\right]=\int_{-\infty}^{\infty}\int_{v}^{\infty}(e^u-e^v)f_{U,V}(u,v)dudv$$
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.
You should study the dynamic of $X_t-Y_t$, don't forget about correlation and that the brownian motion is not the same.
I am pretty sure that there is some big flaws in your model (as taking same interest rate). You should really take a look at this: John C. Hull Options, Futures and Other Derivatives
