Transformation from the Black-Scholes differential equation to the diffusion equation - and back

I know the derivation of the Black-Scholes differential equation and I understand (most of) the solution of the diffusion equation. What I am missing is the transformation from the Black-Scholes differential equation to the diffusion equation (with all the conditions) and back to the original problem.

All the transformations I have seen so far are not very clear or technically demanding (at least by my standards).

My question:
Could you provide me references for a very easily understood, step-by-step solution?

Answer

One starts with the Black-Scholes equation $$\frac{\partial C}{\partial t}+\frac{1}{2}\sigma^2S^2\frac{\partial^2 C}{\partial S^2}+ rS\frac{\partial C}{\partial S}-rC=0,\qquad\qquad\qquad\qquad\qquad(1)$$ supplemented with the terminal and boundary conditions (in the case of a European call) $$C(S,T)=\max(S-K,0),\qquad\qquad\qquad\qquad\qquad\qquad\qquad\qquad(2)$$ $$C(0,t)=0,\qquad C(S,t)\sim S\ \mbox{ as } S\to\infty.\qquad\qquad\qquad\qquad\qquad\qquad$$ The option value $C(S,t)$ is defined over the domain $0< S < \infty$, $0\leq t\leq T$.

Step 1. The equation can be rewritten in the equivalent form $$\frac{\partial C}{\partial t}+\frac{1}{2}\sigma^2\left(S\frac{\partial }{\partial S}\right)^2C+\left(r-\frac{1}{2}\sigma^2\right)S\frac{\partial C}{\partial S}-rC=0.$$ The change of independent variables $$S=e^y,\qquad t=T-\tau$$ results in $$S\frac{\partial }{\partial S}\to\frac{\partial}{\partial y},\qquad \frac{\partial}{\partial t}\to - \frac{\partial}{\partial \tau},$$ so one gets the constant coefficient equation $$\frac{\partial C}{\partial \tau}-\frac{1}{2}\sigma^2\frac{\partial^2 C}{\partial y^2}-\left(r-\frac{1}{2}\sigma^2\right)\frac{\partial C}{\partial y}+rC=0.\qquad\qquad\qquad(3)$$

Step 2. If we replace $C(y,\tau)$ in equation (3) with $u=e^{r\tau}C$, we will obtain that $$\frac{\partial u}{\partial \tau}-\frac{1}{2}\sigma^2\frac{\partial^2 u}{\partial y^2}-\left(r-\frac{1}{2}\sigma^2\right)\frac{\partial u}{\partial y}=0.$$

Step 3. Finally, the substitution $x=y+(r-\sigma^2/2)\tau$ allows us to eliminate the first order term and to reduce the preceding equation to the form $$\frac{\partial u}{\partial \tau}=\frac{1}{2}\sigma^2\frac{\partial^2 u}{\partial x^2}$$ which is the standard heat equation. The function $u(x,\tau)$ is defined for $-\infty < x < \infty$, $0\leq\tau\leq T$. The terminal condition (2) turns into the initial condition $$u(x,0)=u_0(x)=\max(e^{\frac{1}{2}(a+1)x}-e^{\frac{1}{2}(a-1)x},0),$$ where $a=2r/\sigma^2$. The solution of the heat equation is given by the well-known formula $$u(x,\tau)=\frac{1}{\sigma\sqrt{2\pi \tau}}\int_{-\infty}^{\infty} u_0(s)\exp\left(-\frac{(x-s)^2}{2\sigma^2 \tau}\right)ds.$$

Now, if we evaluate the integral with our specific function $u_0$ and return to the old variables $(x,\tau,u)\to(S,t,C)$, we will arrive at the usual Black–Merton-Scholes formula for the value of a European call. The details of the calculation can be found e.g. in The Mathematics of Financial Derivatives by Wilmott, Howison, and Dewynne (see Section 5.4).