Why is Brownian Motion Almost Surely Continuous?

A user from the Quant StackExchange recently asked why the regularity of condition of Brownian motion, namely almost sure continuity, is what it is: almost sure?  Why can't this be upgraded to Brownian motion being surely continuous?

The answer to the latter question is that, actually, it can and very often is.  The answer to the former question is that stipulating almost sure continuity is required in order to make the defining conditions of Brownian motion axiomatic, while still encompassing all of the methods of construction.

Indeed, the most common construction of Brownian motion (or at least the most direct) is through an application of Kolmogorov's extension theorem (the details of this approach can be found in Durrett).  But due to technical issues arising from measure theory (which are actually quite natural), the resulting construction leads to realizations of Brownian motions that are discontinuous.

On the other hand, the approach to constructing Brownian motion from the limit of scaled random walks actually leads to surely continuous realizations.  There are two available routes one can go when having this approach in mind: (a) construct Brownian motion paths directly (i.e. pointwise) from scaled random walks (one common way to do this is by appropriately specifying Brownian motion on the dyadic intervals, interpolating between, and taking limits) or (b) construct Brownian motion by obtaining the Wiener measure on the space of continuous functions beginning at the origin from the induced measures on this space obtained from the scaled random walks on $\mathbb{Z}^{\infty}_{2},$ the space of sequences with values of $-1$ or $1$.

The user also asked whether an explicit example of a discontinuous Brownian motion path could be exhibited.  The following is my complete answer to this and the above questions.

____________________________

Exhibiting a counter-example is straight-forward enough.  For example, let $B_{t}(\omega)$ be a Brownian motion and $\mathcal{T}(\omega)$ a stopping time on $(\Omega,\mathbb{P})$ with a continuous distribution.
Then with
$$B'_{t}(\omega)=\left\{\begin{array}{ll}B_{t}(\omega),&t\neq\mathcal{T}(\omega)\\B_{t}(\omega)+1,&t=\mathcal{T}(\omega),\end{array}\right.$$
$B'_{t}(\omega)$ satisfies (1) and (2) below, but is discontinuous precisely when $t=T(\omega)$.  Therefore, $B_{t}(\omega)$ is a particular realization of Brownian motion that is not everywhere continuous.

There are lots of other ways to obtain a "bad" Brownian motion.  Another example is
$$B'_{t}(\omega)=B_{t}(\omega)\mathbb{1}_{\{B_{t}(\omega)\;\text{irrational}\}},$$
but this is less straight-forward to prove.

____________________________

The reason for stipulating almost sure continuity has to do with the way one constructs Brownian motion, and the issue can be completely dispensed with dependent on one's approach.
The usual presentation in finance texts is the abstract one, namely given a probably space $(\Omega,\mathbb{P})$, one has a Brownian motion $B_{t}(\omega)$ on this space if

1. For every set of times $0\leq t_{1}<t_{2}<\ldots<t_{n}$ the increments $B_{t_{1}},B_{t_{2}}-B_{t_{1}},\ldots,B_{t_{k}}-B_{t_{k-1}}$ form a mutually independent set of random variables on $(\Omega,\mathbb{P}).$
2.  The increments above are normally distributed with mean $0$ and variance $\Delta t$.
3. For almost every $\omega\in\Omega$ the path $t\mapsto B_{t}(\omega)$ is continuous.

Most texts also include a section that sketches a concrete realization of Brownian motion as the limit of scaled random walks.  If one does this rigorously, one sees that (3) upgrades to for every $\omega\in\Omega$

Indeed, if we start with $(\Omega,\mathbb{P})$ satisfying the above and let $\mathcal{P}$ denote the collection of continuous functions $[0,\infty)\to\mathbb{R}$ with $p(0)=0$, then we get from (3) the inclusion map
$$\mathcal{i}:\Omega\to\mathcal{P},$$
defined on a set $\Omega'\subset\Omega$ of full measure, and the push-forward measure of $\mathbb{P}$ onto $\mathcal{P}$ under this inclusion map turns out to be equal to the Wiener measure $\mathbb{W}$ on $\mathcal{P}$, which is unique.

Conversely, one can construct $(\mathcal{P},\mathbb{W})$ directly by starting with the set $\mathcal{P}$ (where every element of this set is continuous a priori) and demonstrating that the measures $\mu_{N}$ on $\mathbb{Z}^{\infty}_{2}$ arising from the appropriately scaled random walks $S_{t}^{N}(\omega)$ ($\omega\in\mathbb{Z}^{\infty}_{2})$ induce a collection of tight measures on $\mathcal{P}$ which converge weakly to $\mathbb{W}$:
$$\mu_{N}\Longrightarrow\mathbb{W}\;\text{(weakly)}$$
One then defines
$$\tilde{B}_{t}(\omega):=p(t)\in\mathcal{P}$$
and readily shows that under $\mathbb{W}$, $\tilde{B}_{t}$ satisfies (1)-(3) and that therefore
$$\tilde{B}_{t}(\omega)=B_{t}(\omega),$$
but that now *every* Brownian motion is continuous.

The equivalence of the implications above show the existence of Brownian motion is essentially tantamount to the existence of a Wiener measure on $\mathbb{W}$ arising from the sequence of measures arising naturally from the scaled random walks.  If one starts from the goal of obtaining this measure, one gets continuity for *every* Brownian motion $p(t)=B_{t}(\omega)$.

____________________________

Other constructions of Brownian motion require us stipulate almost sure continuity due to technicalities arising from measure theory on product spaces.  The quickest construction of Brownian motion in this direction is by applying Kolmogorov's extension theorem on a suitable class of processes; details can be found in Durrett.

A Simple Monte Carlo Simulator for European Call Options

In this post I will present a procedural C++ implementation of a simple Monte Carlo simulator for the pricing of a European call option.  Subsequent articles will make significant improvements such as the pricing of puts and different types of options, improved sampling, incorporation of jump processes, etc.

Recall that we model a stock price $S$ as a stochastic process $\{S_{t}:0\leq t\leq T\}$ governed by the stochastic differential equation (geometric Brownian motion)
$$(1)\;\;\;\;dS_{t}=\mu S_{t}dt+\sigma S_{t}dW_{t}$$
where $\mu$ is the drift or expected return on the stock, $W_{t}$ is a Weiner process $(dW_{t}=N(0,1)\sqrt{dt}$ where $N(0,1)$ is the standard normal distribution), and $\sigma$ is the volatility of the stock (a measure of the variance since $\sigma N(0,1)=N(\sigma^{2},1)$).  Using arbitrage arguments (Merton's method) we arrive at the (linear) Black-Scholes-Merton partial differential equation
$$(2)\;\;\;\;V_{t}+\frac{1}{2}\sigma^{2}S^{2}V_{ss}+rSV_{s}-rV=0.$$ The price $V$ of any derivative must satisfy this equation, and conversely, any solution to (2) gives the price of some derivative based on the boundary conditions used (for European call options, the boundary conditions at the expiry time $T$ are $V(\cdot,T)=\max(S_{T}-K,0)$ where $K$ is the strike price of the option and $S_{T}$ is the value of the (random) stock price at time $T$).  Pricing methods for derivatives based on (2) are known as PDE methods in mathematical finance, and usually consist of numerically solving (2) using finite difference methods (since the domain of definition is rectangular in the $(S,t)$ coordinate system).  This can be difficult, however, since the varied boundary conditions must be taken into account and issues of convergence, stability, etc. enter in.

A different (and more recently popular) approach is probabilistic and involves sampling the randomness in the geometric Brownian motion model (1) of the stock price multiple times and taking an average.  Indeed, an application of the Feynman-Kac formula shows that when discounted appropriately, solutions to (1) are martingales and this implies that $V$ is simply the expected value of the discounted payoff of the derivative.  For a European call option, the payoff is $f(S)=\max(S-K,0)$ and so we get
$$(3)\;\;\;\;V_{\text{EuroCall}}=e^{-rT}\mathcal{E}[f(S_{T})]$$
where $r$ is the riskless return, $\mathcal{E}$ is taken under the risk-neutral probability measure, and $S_{T}$ is the final value of the stochastic process $\{S_{t}\}$ at expiry of the option $T$.  (From now on, $V$ is for the price of a European call option.)  Thus, beginning with (1), passing to the logarithm, applying Ito's Lemma and then applying (3) we get
$$\begin{align*}
&dS_{t}=rS_{t}\;dt+\sigma S_{t}dW_{t}\\
\Longrightarrow&d\log S_{t}=(r-\frac{1}{2}\sigma^{2})dt+\sigma dW_{t}\\
\Longrightarrow&\log S_{t}=\log S_{0}+(r-\frac{1}{2}\sigma^{2})t+\sigma W_{t}\\
\Longrightarrow&S_{t}=S_{0}\exp\left\{\left(r-\frac{1}{2}\sigma^{2}\right)t+\sigma\sqrt{T}N(0,1)\right\}\\
\Longrightarrow&V=e^{-rT}\mathcal{E}\left[f\left(S_{0}\exp\left\{\left(r-\frac{1}{2}\sigma^{2}\right)t+\sigma\sqrt{T}N(0,1)\right\}\right)\right]\\
(4) \Longrightarrow&V=e^{-rT}\mathcal{E}\left[\max\left(S_{0}\exp\left\{\left(r-\frac{1}{2}\sigma^{2}\right)t+\sigma\sqrt{T}N(0,1)\right\}-K,0\right)\right]\\
\end{align*}$$

(Recall that $dW_{t}=N(0,1)\sqrt{dt}$ and so $W_{T}=N(0,1)\sqrt{T}$).  The above computations are somewhat complicated when carried out in detail, and while the explanation for how (2) is derived from (1) is relatively straight-forward, the derivation of (3) from (1) (and thus the previous computations) is significantly more subtle but far more important in the theory of mathematical finance.  A full explanation of the concepts involved (properties of Weiner processes, stochastic integration and Ito's lemma, risk-neutral valuation and the risk-neutral measure, change of numeraire, etc.) is beyond the scope of this post but will be elaborated on in subsequent posts.  In any event, we only need to take faith that (4) is equivalent to (2) for pricing European call options.

The idea of Monte Carlo simulation is now evident. We simulate $N$ paths of $S_{T}$ by sampling from the standard normal distribution $N(0,1)$ and then compute $f(S^{n}_{T})$.  Since each sampling is independent, the random variables $\{f(S_{T}^{n})\}_{n}$ are independent and identically distributed and so the law of large numbers implies
$$\frac{1}{N}\sum_{n=1}^{N}f(S^{n}_{T})\to\mathcal{E}[f(S_{T})]\;\text{as}\;N\to\infty$$
pointwise (and of course, in probability).  We then then discount at $e^{-rT}$ and this is our estimate on the price $V$.

The following implementation is coded procedurally in C++ and prices a European call option using the method explained above.  The only part which may require explanation is the method in which $N(0,1)$ is sampled.  While there are several ways to do this, a simple yet accurate/efficient method is known as the Box Muller algorithm, implemented here.

MonteCarloMethod_Source.cpp

#include <iostream>
#include <cmath>

using namespace std;

double BoxMullerGaussian();
double MonteCarloSimulator(double Expiry, double Strike, double Spot, double Vol, double Riskless, unsigned long N);

int main(void) {
 double Expiry, Strike, Spot, Vol, Riskless, N = 0;
 cout << "Enter: Expiry, Strike, Spot, Vol, Riskless, # Trials" << endl;
 cin >> Expiry >> Strike >> Spot >> Vol >> Riskless >> N;

 cout <<  "The price of the option is "
  << MonteCarloSimulator(Expiry, Strike, Spot, Vol, Riskless, N)
  << endl;
 
 double pause;
 cin >> pause;

 return 0;
}

double BoxMullerGaussian() {
 double x,y;
 
 double sizeSquared;
 do {
  x = 2.0*rand()/(double)(RAND_MAX) - 1;
  y = 2.0*rand()/(double)(RAND_MAX) - 1;
  sizeSquared = x*x + y*y;
 } while(sizeSquared >= 1.0);

 return x*sqrt(-2*log(sizeSquared)/sizeSquared);
}

double MonteCarloSimulator(double Expiry, double Strike, double Spot, double Vol, double Riskless, unsigned long N) {
 double var = Vol*Vol*Expiry;
 double std = sqrt(var);
 double itoCorrection = -0.5*var;
 double movedSpot = Spot*exp(Riskless*Expiry + itoCorrection);
 double _Spot = 0;
 double runningSum = 0;

 for (unsigned long i=0; i<N; i++) {
  double _Gaussian = GetOneGaussianByBoxMuller();
  _Spot = movedSpot*exp(std*_Gaussian);
  double _Payoff = _Spot - Strike;
  _Payoff = _Payoff>0 ? _Payoff : 0;
  runningSum += _Payoff;
 }

 return exp(-Riskless*Expiry)*(runningSum/N); // expectation
}

As an example, we consider Microsoft's stock. Today it closed at 41.35, so $\text{Spot}=41.35$ (the spot price in option pricing is $S_{0}$, which is not random, but $S_{t}$ is for $t>0$ upto $t=T=\text{Expiry}$). We write the option with a strike price $\text{Strike}=40.00$ so that the option is in the money initially (for the buyer). If the time to expiry is 6 mo., then $\text{Expiry}=0.5$. The 6mo. risk-free rate is $\text{Riskless}=0.33$ (taken to be the LIBOR rate quoted this week). The implied 6mo. volatility for Microsoft options is around $\text{Vol}=0.22$. With these parameters and $N=100$ trials, we get $$V_{\text{MSFTEuroCall}}=8.98$$

Put-Call Parity

Put-Call Parity for European Options.  Fix $t>0$ and let $T>t$ be a fixed future time. Denote the continuously compounded risk-free interest rate of tenor $T-t$ at time $t$ by $r_{t}(t,T)$, and let $K$ be the strike price on some asset $S$ negotiated at time $t=0$, whose price at time $t\geq0$ is denoted by $S_{t}$.  Then if $c_{t}(K,T)$ and $p_{t}(K,T)$ are the respective prices of European call and put options on $S$ with strike $K$ and expiry $T$, then we have the following result: $$c_{t}(K,T)-p_{t}(K,T)=e^{-r_{t}(t,T)(T-t)}\left(F_{t}(T)-K\right),$$ where $F_{t}(T)$ is the price of a forward contract on $S$ expiring at time $T$, calculated at time $t$. Theoretically, $$F_{t}(T)=S_{t}e^{-r_{t}(t,T)(T-t)}.$$ We are assuming throughout that the asset $S$ pays no income and has no further funding/carrying cost over the period $[t,T]$ beyond the risk-free rate $r_{t}(t,T)$.  The usual adjustments apply if there are dividends, storage costs, etc., with income benefiting the short party and funding/carrying costs benefiting the long party.

Proof.  The proof consists of a simple replication argument.  Let $V_{t}=c_{t}-p_{t}$ be the value of the portfolio at time $t$ consisting of a long position on $c$ and a short position on $p$.  At time $T$ the payoff from our portfolio is
$$V_{T}=c_{T}-p_{T}=\max(S_{T}-K,0)-\max(K-S_{T},0)=S_{T}-K.$$
Therefore, the payoff from our contract is equal to a long forward position on $S$.  If we took opposite positions, then we would have a short forward position.

Thus our portfolio replicates the payoff of a forward contract on $S$.  It follows that

If the value of a derivative is known at time $T$ with certainty, then the value at any previous time $t$ is equal to the value at time $T$ discounted back to time $t$.  In particular, since $V_{T}$ is known with certainty,
$$V_{t}=e^{-r_{t}(T-t)}(S_{T}-K)\;\;\;\;0\leq t\leq T.$$
This assertion follows from the fact that forward contracts have a certain known payoff at their time of expiration, the fact that our portfolio is equal to this payoff, and the principle of Rational Pricing.  It follows that
$$V_{t}=e^{-r_{t}(T-t)}(F_{t}(T)-K)\;\;\;\;0\leq t\leq T.$$

Overview of the Black-Scholes Model and PDE

In this post we take the PDE approach to pricing derivatives in the Black-Scholes universe.  In a subsequent approach we will cover the risk-neutral valuation approach; the two are essentially equivalent by the Feynman-Kac formula.

I. ASSUMPTIONS

We will assume that our market consists soley of an equity (stock ) $S$, a risk-free money market account (bond) $B$ with return $r\geq0$, and any number of derivatives with $S$ as the underlying.  We will make several non-technical assumptions about $S$ and our market, collectively termed the "Black Scholes market."

1. Infinite liquidity.  Market participants can buy or sell $S$ at any time.

2. Infinite depth.  The buying and selling of $S$ does not affect the price of $S$, no matter the transaction size.

3.  No friction.  It costs nothing to buy or sell $S$ (i.e. trading $S$ incurs no transaction costs).  In particular, the price paid by the buyer and seller is the same (i.e. bid-offer spread it $0$).

4. Constant risk-free rate. $r\equiv\text{const}.$

5. No arbitrage.  There do not exist portfolios of assets where the portfolio is riskless and earns more than the risk-free money account $B$.

6. Infinite divisibility. Market participants can buy $S$ in any amount $\Delta\in\mathbb{R}$.

7. Short selling is possible.  Market participants can short (borrow) $S$ at no cost.

8. No storage costs.  Market participants can hold $S$ at no cost.

It is possible to relax several of these assumptions in various ways.  Short selling is generally permitted in US markets (cf. SEC up-tick rule) and there are of course no storage costs for equities except for the payment of dividends if one has shorted a dividend paying stock.

One can factor variable risk-free interest rates directly into the model by allowing $r$ to be a deterministic function of time, or even to follow a stochastic process.  There is a tremendous amount of on-going research in dealing with  transaction costs, but let us just say that most institutional investors in derivative contracts (e.g. market makers/banks) take huge positions and therefore the effect of transaction costs is minimized (as we will see below, it is the frequency of trades in $S$ in a process termed dynamic hedging that is the source of transaction costs, not the quantity/volume of a single trade).  For similar reasons, the divisibility assumption is also relatively minor when large numbers are involved, since the fractional component of the quantity in a transaction is small relative to the whole number quantity.

The assumption of infinite liquidity is generally not an issue for plain vanilla contracts that actively trade on exchanges, as there are always counter-parties available to take opposing positions.  The assumption becomes more dubious when one moves to the over-the-counter market where exotic contracts are traded; however, the existence of market makers willing to take positions in essentially any contract that can be hedged essentially validates the assumption.  The non-zero bid-offer spread of course then violates the no friction assumption.

The assumption which is most arguable, however, is that of infinite depth, which violates the law of supply and demand.  At the end of the day however, our goal is to develop a model which is both simple and provides a good approximation to reality, and these assumptions will lead us to such a model which is generally good enough, and can be perfected in various ways as needed.

II. ASSET PRICE MODEL - QUANTITATIVE ANALYSIS APPROACH

In order to price derivatives dependent on a stock $S$ (or any underlying for that matter), it is necessary to first come up with a mathematical model for the price movements of $S$.  This is where the approaches of fundamental and quantitative analysis diverge markedly.  Whereas fundamental analysis attempts to predict stock (asset) price movements through a careful analysis of a company's financial statements and other qualitative sources of information like press releases, general market sentiment, new product launches, etc. (this process is often termed equity research and is the method used by traditional long-term value investors), quantitative analysis models stock (asset) price movements with a stochastic model which we now discuss, and attempts to make predictions on future stock price movements through statistical significance analysis/econometrics, which we will not discuss in any detail here.  Incidentally, a third approach, known as technical analysis, attempts to use past data and trends in order to predict the price in the future.  As we will see, the validity of such an approach represents a direct contradiction to the quantitative model we are about to develop, and we therefore assume it is invalid.  There is in fact, little evidence to suggest that technical analysis leads to reliably predictable results of any degree.

We assume the existence of a probability space $(\Omega,\mathbb{P},\mathcal{F})$ and an associated filtration $\{\mathcal{F}\}_{t\geq0}$ that supports a Brownian motion $\{W(\omega,t)\}_{\omega\in\Omega,t\geq0}.$  We will assume that the reader is familiar with the construction and basic properties of $W(\omega,t)$.

Associated to any process $\Delta(\omega,t)\in L^{2}(\Omega\times\mathbb{R}_{t})$ is an integral called the Ito integral

$$I(\omega,t)=\int_{0}^{t}\Delta(\omega,t)\;dW(\omega,t),$$

which we also assume the reader is familiar with (it is defined for each fixed $\omega\in\Omega$ just like the usual Riemann-Stieltjes integral, except that the sample point in the approximating sum is always taken to be the left-end of each partition interval).  In what follows, we will generally suppress reference to $\omega\in\Omega$ and sometimes use the more customary subscript notation $W_{t}, I_{t}, \Delta_{t},$ etc.

An Ito process is a stochastic process $\{X_{t}\}_{t\geq0}$ defined by the stochastic integral equation (SDE)

$$X(t)=X(0)+\int_{0}^{t}a(X(s),s)\;ds+\int_{0}^{t}b(X(s),s)\;dW(s).$$

Note that the terminology stochastic refers to the fact that the equation involves an Ito integral; however, the solution to an SDE is a random process, and so the method of solution of an SDE in principle does not involve any probability theory.  It is the fact that $W(s)$ is almost surely non-differentiable that complicates any proposed solution method (this is expressed by Ito's lemma, which is a resultantly more complicated version of the chain rule).

For $\mu,\sigma\geq0$, if $a=\mu X$ and $b=\sigma X$, we get (with $S=X$)

$$S(t)=S(0)+\mu\int_{0}^{t}S(s)\;ds+\sigma\int_{0}^{t}S(s)\;dW(s).$$

This is our model for stock-price movements, and it is often called geometric Brownian motion.  It is often expressed locally in "differential" form as

$$dS=\mu Sdt+\sigma SdW$$

or

$$\frac{dS}{S}=\mu dt+\sigma dW.$$

There is no harm in doing this as long as one understands vividly that the differential form has no rigorous meaning attached to it, whereas the integral form is well-defined mathematically.  In fact, the entire topics of SDE's should be replaced by SIE's, since $W$ is no-where differentiable and so differential equations involving it really make no sense.  Nonetheless, the convention pervades the subject, and in particular mathematical finance, and so one must get accustomed to it.

There are two reasons why the differential representation of geometric Brownian motion is used.  First, it provides better intuition for why the model is a good model for stock price movements, as we will explain shortly.  Second, as we will see through the remainder of this article, the it facilitates the computations which come up in quantitative finance, in particular those involving Ito's lemma.

Why is geometric Brownian motion a good model for stock price movements?  The two parameters $\mu$ and $\sigma$ represent drift (expected return) and standard deviation per unit time (volatility), respectively.

Geometric Brownian motion captures the intuitive idea that asset prices should drift according to the expected return $\mu$; the riskier the asset, the higher the expected return demanded by investors, and therefore the greater drift in the asset price, regardless of any random fluctuations.  The expected return on an asset is also independent of the stock price, i.e. investors will demand $\mu$ whether the asset trades at $50$ or $5$.  This is modeled by the $\mu\int_{0}^{t}S(s)ds$ term, or $\mu Sdt.$  If there was no stochastic component to the model, then we would have

$$S(t)=S(0)+\mu\int_{0}^{t}S(s)\;ds,$$

or by the fundamental theorem of calculus

$$S'(t)=\mu S(t),$$

which has the solution

$$S(t)=S(0)e^{\mu t}.$$

This is how we would expect the price of a riskless asset with return $\mu$ to grow.

Of course, asset prices are anything but deterministic, and it is natural to assume the presence of an unbiased (i.e. centered at $0$) white noise weighted according to the perceived volatility of the asset (volatility is not a risk-meaure; it is, among other things, directly tied to the liquidity of the asset and how trading affects its price).  The asset price swings should also be directly proportional the the price of the asset itself.  For instance, a stock that trades around $1$ will have price swings markedly lower than an asset that trades around $100.$  This combined with the volatility of the asset is modeled by the $\sigma\int_{0}^{t}S(s)\;dW(s)$ term, or $\sigma SdW.$

If $\mu=0$, then we would have

$$S(t)=S(0)+\sigma\int_{0}^{t}S(s)\;dW(s).$$

There is no corresponding fundamental theorem of calculus for the Ito integral; however, we will see how to solve this below by using Ito's lemma.  In any event, if we combine these two terms, we recover the geometric Brownian motion model

$$S(t)=S(0)+\mu\int_{0}^{t}S(s)\;ds+\sigma\int_{0}^{t}S(s)\;dW(s).$$

III.  EXTENDING THE MODEL TO DERIVATIVES

The extension of the model to derivatives, that is functions of the asset price $S$ and time $t$, involves a technical theorem known as Ito's lemma.  An expository article on its motivation and rigorous proof can be found on one of my previous posts entitled a rigorous proof of Ito's lemma.

If $V=V(S(t),t)$ is the value of a contingent claim on $S$, then we have from Ito's lemma (appropriately extended to handle geometric Brownian motion) that $V$ follows the process
$$\begin{align*}
V(t)&=V(0)+\int_{0}^{t}V_{t}(S(s),t)\;ds+\int_{0}^{t}V_{S}(S(s),s)\;dW(s)+\frac{1}{2}\int_{0}^{t}V_{SS}(S(s),s)\;ds\\
&=V(0)+\mu\int_{0}^{t}S(s)V_{S}(S(s),s)\;ds+\sigma\int_{0}^{t}S(s)V_{S}(S(s),s)\;dW(s)+\int_{0}^{t}V_{t}(S(s),s)\;ds+\frac{1}{2}\sigma^{2}\int_{0}^{t}S^{2}(s)V_{SS}(S(s),s)\;ds\\
&=V(0)+\int_{0}^{t}\left(\mu S(s)V_{S}+\frac{1}{2}\sigma^{2}S^{2}(s)V_{SS}(S(s),s)+V_{t}(S(s),s)\right)\;ds+\sigma\int_{0}^{t}S(s)V_{S}(S(s),s)\;dW(s).\end{align*}$$

In the more usual differential form, we have then
$$dV=\left(\mu SV_{S}+\frac{1}{2}\sigma^{2}S^{2}V_{SS}+V_{t}\right)dt+\sigma SV_{S}dW.$$

We again point out that the differential form is just a short-hand for the integral form above, which is the only mathematically meaningful expression.

VI. DERIVING A NO-ARBITRAGE CONDITION ON $V(S(t),t): THE BLACK-SCHOLES PDE

In this section $\Pi$ denotes the value of a portfolio consisting of $\Delta$ units of an asset $S$, whose value $\{S_{t}\}_{t\geq0}$ follows the geometric Brownian motion process discussed previously, and a short position in an derivative whose underlying is $S$.

As we will see, the exact nature of the derivative and its payoff (the value of $V$ at expiry, i.e. $V(S(T),T)$) function is unimportant, just that its value at time $0\leq t\leq T$ depends only $(S(t),t)$.  In other words, $V$ is path independent. The PDE $V$ must satisfy is the same whether $V$ is the value function for an option, a forward agreement, future, whatever.  In the subsequent posts we will discuss briefly on how one might extend the PDE to cover more exotic types of derivatives with path dependent payoffs or other exotic features written into their contracts like barriers and knock-outs.  We will also assume that $S$ provides no income over the life of the derivative.  In subsequent posts we will discuss how to incorporate income (i.e. dividends if $S$ is a stock), and in any event the modification is trivial (just replace $r$ with $r-q$ below, $q$ being the

The idea of the derivation is to determine the initial capital (cost) that must be put up to construct the portfolio $\Pi$ and hedge the position (risk) so as to ensure a riskless payoff, which must be equal to the payoff of the risk-free money market account $P_{0}e^{rt}$; otherwise there would exist simple arbitrage strategies.  This puts a condition on $\Delta$, the quantity of the underlying we must hold at time $t$ in order to hedge against a short position in $V$.  The strategy is known as continuous time delta hedging.  The hedging eliminates the risk of the portfolio in real time, and as mentioned, means the portfolio must earn the risk free rate.  It is really quite phenomenal that this is possible, and the basic reason that it is has to do with the fact that both $S$ and $V$ are affected by the same underlying uncertainty, namely the Brownian motion process $\{W_{t}\}_{t\geq0}$.  Using Ito's lemma, we can eliminate the terms involving $W$ and therefore eliminate uncertainty.  Indeed, it is impossible to do this if we assume more complicated asset price models different from geometric Brownian motion (this is an active area of research since there is much evidence to suggest that the true process followed by an asset $S$ is a Levy flight, or a fat tailed Brownian motion).

From the way we have constructed our portfolio and the previous sections, we have at time $0\leq t\leq T$ we have
$$\begin{align*}
\Pi(t)&=\Delta(t)S(t)-V(S(t),t)\\
&=\Delta(t)\left\{S(0)+\int_{0}^{t}\mu S(s)\;ds+\int_{0}^{t}\sigma S(s)\;dW(s)\right\}\\
&\;\;\;\;-\left\{V(0)+\int_{0}^{t}\left(\mu S(s)V_{S}(S(s),s)+\frac{1}{2}\sigma^{2}S^{2}(s)V_{SS}(S(s),s)+V_{t}(S(s),s)\right)\;ds+\int_{0}^{t}\sigma S(s)V_{S}(S(s),s)\;dW(s)\right\}\\
&=\Delta(t)S(0)-V(0)\\
&\;\;\;\;+\int_{0}^{t}\left[\Delta(t)\mu S(s)-\left(\mu S(s)V_{S}(S(s),s)+\frac{1}{2}\sigma^{2}S^{2}(s)V_{SS}(S(s),s)+V_{t}(S(s),s)\right)\right]ds\\\
&\;\;\;\;+\int_{0}^{t}\left[\Delta(t)\sigma S(s)-\sigma S(s)V_{S}(S(s),s)\right]dW(s).
\end{align*}$$

We seek to eliminate the uncertainty in $\Pi$, so we match for each time
$$\Delta(t)=V_{S}(S(t),t).$$
(Observe very carefully below how this substitution works!)

Thus, $\Pi$ being riskless, we must have in order to preclude arbitrage opportunities with this portfolio,
$$\begin{align*}
\Pi(t)&=V_{S}(S(0),0)S(0)-V(0)\\
&\;\;\;\;+\int_{0}^{t}\left[\left(\mu V_{S}(S(s),s)S(s)-\mu S(s)V_{S}(S(s),s)\right)-\left(\frac{1}{2}\sigma^{2}S^{2}(s)V_{SS}(S(s),s)+V_{t}(S(s),s)\right)\right]ds\\\
&\;\;\;\;+\int_{0}^{t}\left[\sigma V_{S}(S(s),s)S(s)-\sigma S(s)V_{S}(S(s),s)\right]dW(s).\\
&=V_{S}(S(0),0)S(0)-V(0)-\int_{0}^{t}\left[\frac{1}{2}\sigma^{2}S^{2}(s)V_{SS}(S(s),s)+V_{t}(S(s),s)\right]ds\\
&=\Pi(0)e^{rt}\\

&=\text{value of risk free money market account at time}\;t\;\text{with initial investment}\;\Pi(0).
\end{align*}$$

Note that with the non-differentiable Brownian motion now gone, we can legitimately pass to the differential form of this expression by using the fundamental theorem of calculus.  Differentiating with respect to $t$ we have the localized version of $\Pi$ given by
$$\Pi_{t}=-\left(V_{t}+\frac{1}{2}\sigma^{2}S^{2}V_{SS}\right)=r\Pi(0)e^{rt}=r\Pi=r(SV_{S}-V).$$
Consequently,
$$V_{t}+\frac{1}{2}\sigma^{2}S^{2}V_{SS}+rSV=rV.$$
Note that we also have
$$V_{t}+\frac{1}{2}\sigma^{2}S^{2}V_{SS}=-r\Pi_{0}e^{rt}.$$
Of coruse we don't know $\Pi_{0}$ a priori, so we appeal to the previous PDE, which is known as the Black-Scholes PDE.  Observe that despite $S$ being present, there is nothing random about the PDE and $S$ can be regarded as as dummy variable.  It is the no-arbitrage requirement that leads to a deterministic outcome, and at each time/asset price $(S(t),t)$ we have the price given by the solution to the PDE, with remaining time to expiry beign $T-t$.

It is backwards parabolic, and so requires a terminal condition $f$ and time $t=T$, the payoff for the specific derivative under consideration.  Additional boundary conditions will lead to prices for other more exotic options discussed in subsequent posts.

Analysis of the Black Scholes PDE

In this post we conduct a cursory analysis of the Black-Scholes (B-S) partial differential equation (PDE), including existence and uniqueness of solutions, well-posedness, and in certain special circumstances, analytical solutions.

The B-S PDE is defined the space positive half space $\mathbb{R}^{+}_{S}\times\mathbb{R}^{+}_{t}=\{(S,t): S\in(0,\infty),t\in(0,\infty)\}$ by the equation
$$(1)\;\;\;\;V_{t}+\frac{1}{2}\sigma^{2}S^{2}V_{SS}+rSV_{S}-rV=0.$$

For various pay-off functions $f=f(S,t)$ (where the $S$ in the argument is interpreted as $\{S_{t}\}_{t\in\mathbb{R}}$; in many cases it will just be $t=T$, the value of $S$ at expiration), boundary conditions may be instituted and the domain of definition becomes $(0,S_{\text{max}})\times(0,T)$.  We shall discuss boundary conditions in more detail later; for now we will work on the positive plane.

For $\xi\in\mathbb{R}$ we have
$$-\frac{1}{2}\sigma^{2}\xi^{2}<\theta\xi^{2}$$
for some appropriate positive parameter $\theta$ depending only on $\sigma$.  Thus, (1) is uniformly backward parabolic.

I.  Converting the B-S PDE to the heat equation (systematic procedure)

We begin our analysis by performing a sequence of change of independent and dependent variables in order to reduce the equation to a more simple one.  We first deal with the variable co-efficients, whose special structure suggests the change of variables $S\mapsto e^{-S}$, or using a new variable name, $x=\log S$.  The differential operators then become
$$\frac{\partial}{\partial S}=\frac{\partial}{\partial x}\frac{\partial{x}}{\partial S}=\frac{1}{S}\frac{\partial}{\partial x}$$
and
$$\frac{\partial^{2}}{\partial S^{2}}=\frac{\partial}{\partial S}\left[\frac{1}{S}\frac{\partial}{\partial x}\right]=\frac{1}{S^{2}}\left(\frac{\partial}{\partial x^{2}}-\frac{\partial}{\partial x}\right).$$
Substituting these into (1) yields
$$V_{t}+\frac{1}{2}\sigma^{2}(V_{xx}-V_{x})+rV_{x}-rV=0,$$
or
$$V_{t}+\frac{1}{2}\sigma^{2}V_{xx}+(r-\frac{1}{2}\sigma^{2})V_{x}-rV=0.$$

The reaction term $rV$ will lead to an exponential increase in the solution, and so we assume $V$ has the form $e^{rt}u$ for some to be determined function $u$.  This leads to the change of variables $V\mapsto e^{rt}u$ and the equation becomes
$$e^{rt}\left(ru+u_{t}+\frac{1}{2}\sigma^{2}u_{xx}+(r-\frac{1}{2}\sigma^{2})u_{x}-ru\right)=0,$$
or
$$u_{t}+\frac{1}{2}\sigma^{2}u_{xx}+(r-\frac{1}{2}\sigma^{2})u_{x}=0.$$

(If $r=r(t)$ is not constant, but a deterministic function of $t$ alone, then we can solve the ODE $w_{t}-r(t)w=0$ and make the change of variables $V=w(t)u$ in order to eliminate the source term; we will assume however that $r\equiv\text{const}$ for the remainder of this article.)

We now deal with the drift term $(r-\frac{1}{2}\sigma^{2})u_{x}$ by switching to the moving frame $x'=x-(r-\frac{1}{2}\sigma^{2})t$ (i.e. switching to the characteristic coordinate system).  We also at this stage change the structure of the equation to forward parabolic by changing the direction of evolution to forward time through the change of variables $t\mapsto -t$, or $t'=-t.$  We compute as before the new differential operators
$$\frac{\partial}{\partial x}=\frac{\partial}{\partial x'}\frac{\partial x'}{\partial x}+\frac{\partial}{\partial t'}\frac{\partial t'}{\partial x}=\frac{\partial}{\partial x'},$$
$$\frac{\partial^{2}}{\partial x^{2}}=\frac{\partial^{2}}{\partial x'^{2}},$$
$$\frac{\partial}{\partial t}=\frac{\partial}{\partial x'}\frac{\partial x'}{\partial t}+\frac{\partial}{\partial t'}\frac{\partial t'}{\partial t}=-(r-\frac{1}{2}\sigma^{2})\frac{\partial}{\partial x'}-\frac{\partial}{\partial t'}.$$

Substitution then yields
$$-(r-\frac{1}{2}\sigma^{2})u_{x'}-u_{t'}+\frac{1}{2}\sigma^{2}u_{x'x'}+(r-\frac{1}{2}\sigma^{2})u_{x'}=0$$
or
$$(2)\;\;\;\;u_{t'}-ku_{x'x'}=0.$$

This is the heat equation, with $k=\frac{1}{2}\sigma^{2}.$  It therefore suffices to study the properties of (1) by studying (2).

II. Uniqueness and Stability

We shall restrict our attention to the domain $D=\{(x,t):0\leq t\leq T,0\leq x\leq \ell\}.$ This is consistent with domain used in practice with solving (1).  We have the following theorem

Maximum Principle.  If $u$ solves the heat equation $u_{t}=ku_{xx}$, then $u$ achieves its maximum on the boundary lines $t=0$, $x=0$, or $x=\ell$.  Moreover, if $u$ achieves its maximum in the interior of $D$ or on the line $t=T$, then $u$ is constant.

The first assertion is what is known as the weak maximum principle, and this is what we shall prove.  The second assertion is called the strong maximum principle since it implies the weak version; however, its proof requires tools from the theory of harmonic functions (e.g. the mean-value property) and some non-trivial analysis (cf. Evans).  Analogous statements hold for the minimum by considering $-u$.

Proof.  Recall from calculus that if a function $u$ obtains its maximum at an interior point of a set $E$, then $D_{i}u=0$ and $D^{2}_{i}u\leq0$.  Therefore we have $u_{t}=0$ and $u_{xx}\leq0$.  The idea is to show that $u_{xx}<0$ in order to obtain a contradiction, and we will do this with a perturbation.

Let $\epsilon>0$, put $M=\sup_{(x,t)\in\partial D-\{(x,t):t=T\}}u$, and define $v:=u+\epsilon x^{2}.$ The definition of $v$ shows
$$(3)\;\;\;\;v_{t}-kv_{xx}=u_{t}-ku_{xx}-2\epsilon k=-2\epsilon k<0.$$
Applying the argument in the first paragraph shows that $v$ does not attain its maximum in $D^{\circ}$.  Suppose now $u$ obtains its maximum at a point $(x_{0},T)$ on the line $t=T$.  Then $v_{x}=0,v_{xx}\leq0$ as before and
$$v_{t}(x_{0},T)=\lim_{\delta\to0^{+}}\frac{v(x_{0},T)-v(x_{0},T-\delta)}{\delta}\geq0,$$
a contradiction.  Since $D$ is compact, the continuous function $v$ must obtain its maximum on $\partial D-\{(x,t):t=T\}.$  Therefore $v\leq M+\epsilon\ell^{2}$ on $D$.  By sending $\epsilon\to0$, we obtain the assertion of the weak maximum principle for $u$.

This theorem immediately settles the uniqueness of the boundary value problem for (2), for suppose $u$ and $v$ are both solutions to (2) with identical boundary data. Then $w:=u-v$ is also a solution by linearity has vanishes on the boundary $x=0,x=\ell,t=0$.  Thus $w\equiv0$ in all of $D$ by combining the assertions of the maximum and minimum principles.  This implies $u\equiv v$ and the claim follows.

Another important route to establishing uniqueness is through the concept of conservation of energy.

Consider the "energy" integral
$$E(t)=\int_{0}^{\ell}|u(x,t)|^{2}\;dx$$
where $u$ vanishes on the spatial boundary $0\leq x\leq\ell$ and at the intial time $t=0$.  Differentiating with respect to $t$ yields
$$E'(t)=\frac{d}{dt}\int_{0}^{\ell}u^{2}\;dx=\int_{0}^{\ell}2uu_{t}\;dx.$$
where differentiation under the integral sign is valid, since as will be clear in the next section when we solve the heat equation, the mapping $x\mapsto 2u(x,t)u'(x,t)$ is absolutely integrable on $[0,\ell]$ for every $t\in[0,T)$.  Substituting $u_{t}=ku_{xx}$ and integrating by parts yields
$$E'(t)=2k\int_{0}^{\ell}uu_{xx}\;dx=2kuu_{x}\Big|_{x=0}^{x=\ell}-2k\int_{0}^{\ell}u^{2}_{x}\;dx=-2k\int_{0}^{\ell}u_{x}^{2}\;dx\leq0.$$
However, $E(0)=0$, and therefore $E\equiv0$ for all $t$.  Consequently, $u\equiv0$ and by considering the difference of two solutions with identical boundary data, we recover uniqueness (initially, we have uniqueness only almost everywhere, but since solutions to the heat equation are smooth, we recover pointwise uniqueness.)

We continue our efforts to establish well-posedness by showing that the PDE is stable with respect to the initial and boundary data.

III. Existence and the fundamental solution

There are by now multiple approaches to solving the heat equation recorded in the literature.  One of the simplest makes use of the Fourier transform.  We begin by extending our spatial domain to all of $\mathbb{R}$ and enforcing a decay condition so that $u\in L^{2}(dx).$  Our convention for the (one-dimensional) Fourier transform is
$$\mathcal{F}(u)(\xi):=\hat{u}(\xi)=(2\pi)^{-1/2}\int_{-\infty}^{\infty}u(x)e^{-ix\xi}\;dx$$
so that
$$||\hat{u}||_{2}=||u||_{2}\;\;\;\;(\text{Plancheral})$$
along with the associated inversion formula
$$u(x)=\int_{-\infty}^{\infty}\hat{u}(\xi)e^{ix\xi}\;d\xi.$$

The Fourier transform drops polynomial weights in exchange for anti-differentiation, so we apply it to the spatial derivative $x$ in order to get

$$\hat{u}_{t}=-k|\xi|^{2}\hat{u}.$$

Note that we have used the relation $\mathcal{F}(\partial_{t}u)=\partial_{t}\mathcal{F}u,$ which holds because of our condition on $u$.

Solving this ODE (in $t$) yields
$$\hat{u}=\hat{u}(0)e^{-k\xi^{2}t}.$$

A simple computation involving substituting the above result into the inversion formula and applying the convolution theorem then yields

$$u(x,t)=(4\pi kt)^{-\frac{1}{2}}\int_{-\infty}^{\infty}\phi(x)e^{-(x-y)^{2}/4kt}\;dy,$$
where $\phi(x)=u(x,0).$  This implies that our fundamental solution is
$$\Phi(x,t)=(4\pi kt)^{-\frac{1}{2}}e^{-x^{2}/4kt}.$$

IV. Boundary conditions and special solutions to the B-S PDE

Let us examine some typical boundary conditions in the Black-Scholes PDE.

Explicit Finite Difference Method for Black-Scholes-Merton PDE (European Calls)

Here we will implement a basic finite difference method to the solution of the Black-Scholes-Merton PDE
$$(1)\;\;\;\;\frac{\partial V}{\partial t}+\frac{1}{2}\sigma^{2}S^{2}\frac{\partial^{2}V}{\partial S^{2}}+rS\frac{\partial V}{\partial S}-rV=0,$$
which of course models the value $V$ of any derivative contract in the absence of arbitrage (see the Wikipedia article for a more comprehensive list of assumptions under which the Black-Scholes-Merton model is valid).

This PDE is a backwards diffusion (parabolic) equation: it is ill-posed in forward time.  This means that "initial conditions" are prescribed at $t=T$, the expiry time of the derivative being priced (in the case of this article, an option).  This is analogous to the prototypical heat equation $u_{t}=\Delta u,$ which is ill-posed in backward time, and hence requires initial conditions to be prescribed at time $t=0$.  It makes sense that a reverse-time evolution equation would model the value of a derivative since the payoff of the conract is known at the time of expiry $t=T$, and it is at the initial time $t=0$ that we are interested in knowing the "correct" value of the contract.  Indeed, we assume that the underlying of the derivative is a stock which follows a geometric Browning motion $\{S_{t}\}_{t\geq0}$ as in the Black-Scholes-Merton model.  This is a stochastic process for which we know $S_{0}$ (the time that the contract is being sold) and so the Black-Scholes equation allows us to compute $V(S(t),t)\big|_{t=0}=V(S_{0},0)$, which is the value of the contract today with observed stock price $S_{0}$.

The above discussion gives us an indication as to what boundary conditions must be prescribed in order that the solution $V(S,t)$ corresponds to the particular derivative contract being valued.  Indeed, (1) gives the value of any derivative; it is the boundary conditions which specify which derivative.  We shall consider the PDE to be defined on the $(S,t)$ space $[0,\bar{S}]\times[0,T]$ where $\bar{S}$ is the maximal stock price (typically $\bar{S}$ is taken to be $\bar{S}\approx 3K$, where $K$ is the strike price of the option).  (Ideally, we would have $\bar{S}=\infty$ so that we could apply the natural bounary condition $V(S,t)\to S$ as $S\to\infty$, but we must use a finite domain when applying finite differences.)  We then have the following boundary-value problem corresponding to (1)

$$\left\{\begin{array}{ll}
\frac{\partial V}{\partial t}+\frac{1}{2}\sigma^{2}S^{2}\frac{\partial^{2}V}{\partial S^{2}}+rS\frac{\partial V}{\partial S}-rV=0,&0\leq t< T, 0<S<\bar{S}\\
V(S,T)=\max(S-K,0),&t=T\\
V(0,t)=0,&S=0\\
V(\bar{S},t)=\bar{S}-Ke^{-r(T-t)},&S=\bar{S}.\end{array}\right.
$$

We now discretize the domain into increments $\Delta S$ and $\Delta t$.  The following notation is used: $k=\Delta t$, $h=\Delta S,$ $V_{m}^{n}=V(mh, nk),$ $N=T/k$, $M=\bar{S}/h.$  For simplicity we start with forward differences in time and centered differences in space to get (after some algebra)
$$V_{m}^{n}=A_{m}V_{m-1}^{n+1}+B_{m}V_{m}^{n+1}+C_{m}V_{m+1}^{n+1}$$
where
$$\begin{array}{l}
A_{m}=\frac{k}{2}\left(\sigma^{2}m^{2}-rm\right)\\
B_{m}=(1-k)(r\sigma^{2}m^{2})\\
C_{m}=\frac{k}{2}\left(\sigma^{2}m^{2}+rm\right)\end{array}
$$

Recall that the evolution is backward in time, so that this method is indeed explicit.

The following method was coded in Excel VBA using the method discussed above.  In this implementation we have used parameters $T=0.25$ (expiry in 3 months), $K=10$ (strike price 10 USD), $r=0.1$ and $\sigma=0.4.$

BS_FDM_Explicit.xlsm (Coded in Excel VBA)

Sub BSFiniteDifferenceScheme()

    ' User defined parameters
    Dim Riskless As Double
    Dim Vol As Double
    Dim Expiry As Double
    Dim Strike As Double
    
    Riskless = 0.1
    Vol = 0.4
    Expiry = 0.25
    Strike = 10
    
    ' Finite difference parameters
    Dim T As Double ' final time
    Dim S As Double ' stock price upper limit
    Dim N As Integer ' temporal steps
    Dim M As Integer ' stock steps
    Dim k As Double ' time step
    Dim h As Double ' stock step
    Dim V() As Double ' solution V(s,t)
    Dim T_Array() As Double
    Dim S_Array() As Double
    
    ' Set preliminary parameters
    T = Expiry
    S = 30
    k = 0.001
    h = 0.5
    N = T / k
    M = S / h
    ReDim V(N, M)
    
    For i = N To 0 Step -1
        Cells(i + 2, 1) = i * k
    Next i
    For j = 0 To M
        Cells(1, j + 2) = j * h
    Next j
    
    ' Set left boundary condition (S=0)
    For i = N To 0 Step -1
        V(i, 0) = 0
        Cells(i + 2, 2) = V(i, 0)
    Next i
    
    ' Set upper boundary condition (t=T)
    For j = 0 To M
        If j * h > Strike Then
            V(N, j) = j * h - Strike
        Else
            V(N, j) = 0
        End If
        Cells(N + 2, j + 2) = V(N, j)
    Next j
    
    ' Set right boundary condition (S=S_max)
    For i = N To 0 Step -1
        V(i, M) = S - Strike * Exp(-Riskless * (T - i * k))
        'V(i, M) = S
        Cells(i + 2, M + 2) = V(i, M)
    Next i
    
    ' Computation Loop
    Dim A As Double
    Dim B As Double
    Dim C As Double
    For i = N - 1 To 0 Step -1
        For j = 1 To M - 1
            A = (k / 2) * (Vol ^ 2 * j ^ 2 - Riskless * j)
            B = 1 - k * (Riskless + Vol ^ 2 * j ^ 2)
            C = (k / 2) * (Vol ^ 2 * j ^ 2 + Riskless * j)
            V(i, j) = A * V(i + 1, j - 1) + B * V(i + 1, j) + C * V(i + 1, j + 1)
            Cells(i + 2, j + 2) = V(i, j)
        Next j
    Next i
   
End Sub

The output has been setup to give the value of $V(S,t)$ for any pair $(S,t)$.  The table of data is of course to large to be given here, and of course we are only interested in the data $V(S_{0},0)$ where $S_{0}$ is the spot price at $t=0$.  For example, when $S_{0}=20$ the computation gives $V(20,0)=10.25$ USD.  This is the price at time $t=0$ for the option to buy a stock at time $T=0.25$ years when the current spot price is $20$ USD and the strike price is $10$ USD.  This rather high price (relative to the stock price) makes sense considering that the option is deep in the money and likely to be exercised at expiry.  To check that our value is indeed correct, we run the Monte Carlo simulation discussed in this post with identical parameters; the computed value with $N=10000$ trials is $V=10.21$ USD.  If we increase $N$ to $N=1000000$ (one million), we get with much greater consistency $V=10.25$, which is in agreement with our finite difference computation.

A graph produced from the output in Excel is given below.

The orientation has been adjusted so that the solution curve at $t=0$ is clearly visible (the one of interest).

Exact pricing formula for a binary put or call

Consider a contract that pays $B=1$ if $S(T)\geq K$ at the terminal time $T$ (and $B=0$ if $S(T)<K$), where $S$ follows the geometric brownian motion process

$$S(t)-S_{0}=\mu\int_{0}^{t}S(t)\;dt+\sigma\int_{0}^{t}S(t)\;dW(t).$$

In "differential" form,

$$dS_{t}=\mu S_{t}dt+\sigma S_{t}dW_{t}.$$

This contract is known as a binary call option.  We seek to determine its price $V$ at the initial time $t=0$.  

The Black-Scholes theory implies (through its connection to parabolic partial differential equations and the Feyman-Kac formula) that this price is equal to the risk-neutral expected payoff $\mathbb{E}_{\mathbb{Q}}$ discounted at the risk-free rate $r$.  That is,

$$V=e^{-rT}\mathbb{E}_{\mathbb{Q}}[B].$$

A simple computation reveals that this is equal to $e^{-rT}\mathbb{Q}(S_{T}>K)$.  In detail,

$$e^{rT}V=\mathbb{E}_{\mathbb{Q}}[B]=1\cdot\mathbb{Q}(B=1)+0\cdot\mathbb{Q}(B=0)=\mathbb{Q}(B=1)=\mathbb{Q}(S_{T}>K).$$

Now, as covered in the Black-Scholes theory, an application of Girsonov's theorem leads us to the equivalent process for $S$
$$dS=rSdt+\sigma Sd\tilde{W_{t}}$$
where $\tilde{W}$ is a Brownian motion under the equivalent martingale measure (risk-neutral measure) $\mathbb{Q}$.  We now proceed as usual taking the natural logarithm to get
$$d\log S=(r-\frac{1}{2}\sigma^{2})dt+\sigma d\tilde{W_{t}},$$
which, in light of the risk-neutral probability computation above, immediately leads us to
$$\mathbb{E}_{\mathbb{Q}}[B]=N(d_{2}),$$
where $N=\Phi(0,1)$, the standard normal density, and
$$d_{2}=\frac{\ln(S_{0}/K)+(r-\frac{1}{2}\sigma^{2})T}{\sigma\sqrt{T}}.$$
Thus,
$$V=e^{-rT}N(d_{2}).$$
The price of the put can be determined by following the steps above, but where the contract has a payoff of $1$ when $S(T)\leq K$.  The result is
$$V=e^{-rT}N(d_{1}),$$
or equivalently
$$V=e^{-rt}N(-d_{2}).$$

As an aside from this, we can compute the price of this option using only the known prices of European call options.  Observe that the payoff function $P$ of this option is given by
$$P=\frac{\partial C_{K}}{\partial S}$$
where $C$ is the payoff function for European call option with strike price $K$.
This observation leads us to the conclusion
$$V=\lim_{\delta\to0}\frac{{U(K-\delta})-U(K)}{\delta}=-\frac{\partial U}{\partial S},$$
where $U(K)$ is the value of the call option with strike $K$.