Zero to Black Scholes PDE

Balaena Quant Insights: Issue 33

This is a small writeup I have made after spending a week and a half trying to derive the black scholes formula in hopes to easily understand how it works for future option strategies. This writeup covers the following topics:

• Binomial trees convergence to Black Scholes
• Brownian motion basics, (dW²) == dt
• Geometric brownian motion, log normal returns
• Itos lemma
• Randomness hedging to Black Scholes PDE

Keep in mind it does not cover out of scope advanced topics such as the solution of PDE via heat, Greeks, IV smile and so on. However the core concepts covered here should be sufficient as prerequisites on further understanding of said advanced topics.

Binomial trees converge to Black Scholes

Say I have a stock price at time t, Sₜ = 100. This stock at t + 1 can move up or down by the same factor 0.2: S_up = 120, S_down = 80. Cash grows by 1.05 per period and is denoted by B (bonds).

Suppose there is a call option out there with strike price K = 100, expiring in 1 time period, and the payoff of such option is 20 in the up state and 0 in the down state. How much should the option be priced to avoid arbitrage against a structured portfolio with the same payout?

Suppose such a portfolio exists. How would it be structured? The portfolio holds Δ shares of the stock and B in bonds, where Δ is the number of shares and B is the amount invested in bonds today. Its value today (time t) is:

Vₜ = ΔSₜ + B

And its value one period later (time t + 1), after the bond grows by 1.05, is:

Vₜ₊₁ = ΔSₜ₊₁ + 1.05B

So at t, the portfolio value is:

Vₜ = 100Δ + B

So for the following time period, matching the payout requirements based on the states, the portfolio returns would be:

120Δ + 1.05B = 20   → up state
 80Δ + 1.05B = 0    → down state

With 2 unknowns, Δ and B, we can solve them directly:

(120 − 80)Δ = 20 − 0
        40Δ = 20
          Δ = 0.5

80(0.5) + 1.05B = 0
     40 + 1.05B = 0
              B = −40 / 1.05
              B = −38.0952

We can substitute the unknowns into the original Vₜ formula to figure out the portfolio value at that time, which is the price of the call option:

Vₜ = 0.5(100) + (−38.0952) ≈ 11.90

C = 11.90

The probability of stock price movement is not used in the calculation. Since we matched the payoff of every state, this replication method will work regardless of the likelihood of the outcomes.

The option’s price is determined by:

• Current stock price
• Possible future stock prices
• Option payoff (match up/down payoff difference)
• Risk free rate (adjust payoff to correct level)

which is the foundation of the binomial option pricing model.

A related concept that uses these values is the risk neutral probability:

q = ((1 + r) − d) / (u − d)

q = risk neutral probability of upward movement
r = risk free interest rate per period (1.05 − 1 = 0.05)
u = up factor (1.2)
d = down factor (0.8)

In our case, q = (1.05 − 0.8) / (1.2 − 0.8) = 0.625.

The call option price using binomial model pricing would be:

C = (q·C_u + (1 − q)·C_d) / (1 + r)
  = (0.625(20) + 0.375(0)) / 1.05
  = 11.90

C_u = up state payoff (single step)
C_d = down state payoff (single step)

Now given a second time period with the same u, d and R, where the final stock prices are 144, 96 and 64, with matured call payoffs of 44, 0 and 0 respectively, we can construct a binomial tree like so:

Since the formulas we have used so far only work for 1 step, to get the option price C at t = 0 we have to work backward from t = 1. Using the binomial option pricing model, calculate q and C for both t = 1 nodes.

The 120 node:

q = ((1 + 0.05) − 0.8) / (1.2 − 0.8) = 0.625

C = (0.625(44) + 0.375(0)) / 1.05 = 26.19   → option value for t = 1, node #1

The 80 node:

q = 0.625 (same)

C = (0 + 0) / 1.05 = 0   → option value for t = 1, node #2

From that, we can calculate the 100 node at t = 0:

q = 0.625 (same)

C = (0.625(26.19) + 0.375(0)) / 1.05 = 15.59   → option value for t = 0

The general binomial sum formula that does all that is:

C₀ = (1 / Rᴺ) · Σⱼ₌₀ᴺ C(N, j) · qʲ · (1 − q)ᴺ⁻ʲ · max(S₀·uʲ·dᴺ⁻ʲ − K, 0)

Where N is the number of periods and C(N, j) is the binomial coefficient “N choose j”. This formula assumes u, d, R and q are constant across all N periods. This creates a recombining symmetric binomial tree.

In a more practical sense, the binomial model has a limit. Over one time period, besides the up and down factors being constant, the model does not tie u and d to the volatility of the underlying stock.

In the naive tree, u and d are inputs. So the volatility of price movement depends on the step size (how deep the tree is). In the CRR tree (Cox Ross Rubinstein), the input is the volatility, and u and d are then adjusted according to the step size. CRR replaces:

u → u = e^(σ√Δt)
d → d = 1/u

Δt = T/N

T = expiry in years
N = interval count

Since CRR scales the up and down step size based on N, assuming T is fixed, it demonstrates a Brownian motion property where the variance converges as N grows; the naive tree’s variance explodes as N grows.

The Brownian motion property is more apparent at high N, which allows the CRR option price to be close to the Black Scholes option price (which also has a Brownian motion assumption). Both models have similar inputs as well. Because of this, Black Scholes is viewed as the continuous version of CRR with a different modeling approach.

CRR and Black Scholes converging

CRR error against Black Scholes as N increases

Brownian motion basics

A standard Brownian motion series W(t) has 4 properties:

1. Starts at zero: W(0) = 0
2. Normally distributed increments: for s < t, the change W(t) − W(s) follows a normal distribution with mean 0 and variance t − s, i.e. W(t) − W(s) ∼ N(0, t − s). Example: if t − s = 4, the increment W(10) − W(6) will have mean 0 and variance 4.
3. Independent increments
4. The path is continuous.

For rule 2, it is also implicitly stated that the variance of increments should ONLY scale with the time difference between observations, not with how frequent the observations are (discrete), which is why naive binomial trees are not considered Brownian motion.

To fulfill rule 2, the step size needs to be scaled with the number of steps for the whole interval. Referring to Δt = T/N, the step size of a standard Brownian motion needs to be √Δt. If we were to use, let’s say, Δt as the step size, the variance of one step would be:

Var(ΔW) = (Δt)²

After time t, we take t/Δt steps, and since the steps are independent, the total variance adds up:

Var(W(t)) = (t/Δt)·(Δt)²
          = t·Δt

As Δt → 0, the variance approaches zero. And since t·Δt ≠ t, the rule is violated. However, if we use √Δt as the step size, the variance of one step would be:

Var(ΔW) = (√Δt)²
        = Δt

And the total variance would be:

Var(W(t)) = (t/Δt)·Δt
          = t

Which is exactly the variance required by Brownian motion.

Brownian motion is random and it is not differentiable like ordinary calculus / smooth functions, which is important for Black Scholes modeling later.

The main difference between Brownian motion calculus and ordinary calculus is the identity rule dW² = dt. dW here represents a small change in the value of the Brownian motion and dt represents a small change in time. The equation

dW² = dt

expresses Brownian motion’s quadratic variation, [W]ₜ = t. Quadratic variation measures the total accumulated squared changes over time. For the quadratic variation of a normal smooth function like a parabola, for example

f(x) = x²  from x = 0 to x = 1

and we calculate its quadratic variation by dividing the function into 4 smaller parts:

We then calculate the function values:

x       f(x) = x²    Difference
0.00    0.0000       —
0.25    0.0625       0.0625
0.50    0.2500       0.1875
0.75    0.5625       0.3125
1.00    1.0000       0.4375

The changes in f(x) between consecutive points are also displayed. With 4 intervals, the approximate quadratic variation is the sum of squared changes:

QV = Σᵢ₌₁ⁿ (Δfᵢ)²

For the example: QV = 0.328125. As we divide the interval into more parts, each individual change gets smaller, so the sum of squared changes approaches zero.

[f]ₜ = 0

For Brownian motion, we can evaluate its quadratic variation over similar steps; start by dividing into smaller intervals:

Δt = T/n

Each Brownian increment satisfies:

ΔWᵢ ∼ N(0, Δt)

Now square each increment and sum them together:

Qₙ = Σᵢ₌₁ⁿ (ΔWᵢ)²

To evaluate that term, we can use the expected value of ΔWᵢ². For any random variable X:

Var(X) = E[X²] − (E[X])²
 E[X²] = Var(X) + (E[X])²

For a Brownian increment, E[ΔWᵢ] (mean) = 0 and Var(ΔWᵢ) = Δt. Thus:

E[ΔWᵢ²] = Δt + 0 = Δt

This means the expected square of each increment equals the length of the time interval. With this result, we can calculate E[Qₙ]:

E[Qₙ] = E[Σᵢ₌₁ⁿ (ΔWᵢ)²]
      = nΔt = T

The interpretation is that no matter how finely you divide the interval, the expected sum of squared increments is always T. With the expected value done, we need the variance of Qₙ as n → ∞ to complete the quadratic variation.

Get Jun Han’s stories in your inbox

Join Medium for free to get updates from this writer.

Remember me for faster sign in

For a normal random variable X ∼ N(0, σ²), its fourth moment is:

E[X⁴] = 3σ⁴

Using the variance formula:

Var(X²) = E[X⁴] − (E[X²])²
        = 3σ⁴ − (σ²)² = 2σ⁴

For a Brownian increment, σ² = Δt. So:

Var((ΔWᵢ)²) = 2(Δt)²

Since increments are independent, the variances of the squares add:

Var(Qₙ) = Σᵢ₌₁ⁿ Var((ΔWᵢ)²) = n · 2(Δt)²

Replace Δt = T/n:

Var(Qₙ) = 2n(T/n)²
        = 2T²/n

Var(Qₙ) = 2T²/n → 0  as n → ∞

As we make the time intervals smaller, the variance of the total sum of squared increments approaches zero. Given two facts:

E[Qₙ] = T,   Var(Qₙ) = 2T²/n → 0

The mean remains T, while the variance approaches zero. Therefore, the sum converges to T in mean square. In probability terms:

Σᵢ₌₁ⁿ (ΔWᵢ)² → T

Which means:

• each individual increment is random
• the square of each individual increment is random
• when we add many small squared increments, it becomes less random
• the total approaches the deterministic length of the interval

This also gives us the shorthand dW² = dt. With that being said, the main differences between Brownian motion and a smooth function are:

Property                    Smooth function    Brownian motion
typical increment           Δt                 √Δt
sum of squared increments   → 0                → T
quadratic variation         0                  T

Itô’s lemma

Why is this distinction important? When doing calculus, we have to keep a second order term that ordinary calculus mostly discards, such as in the 2nd order Taylor expansion below for the function f(Wₜ):

df = f′(Wₜ)dWₜ + ½f″(Wₜ)(dWₜ)²

Apply the Itô rule (dW)² = dt and we get:

df = f′(Wₜ)dWₜ + ½f″(Wₜ)dt

This is the core concept of Itô’s lemma. Other Itô rules are as such:

Expression      Itô rule
(dt)²           0
dt·dWₜ          0
(dWₜ)²          dt
dWₜ·dt          0
dWₜ²            dt
d(f(Wₜ))        f′(Wₜ)dWₜ + ½f″(Wₜ)dt
d(f(t, Wₜ))     (∂f/∂t)dt + (∂f/∂W)dWₜ + ½(∂²f/∂W²)dt

The purpose of a Taylor expansion is to approximate a smooth function near a point using its derivatives. E.g. for f(x) = x², if we want to expand around x = 3, its Taylor expression is f(3 + h) = 9 + 6h + h². Putting h = 1, we are essentially computing f(4) or f(3 + 1), which both equal 16. The results from Taylor’s expansion are used in other calculus methods such as the chain rule.

However, an ordinary differential (normal Taylor, normal chain rule) describes how a function changes due to deterministic inputs, while a stochastic differential describes how a function changes when the inputs include random processes like Brownian motion. This is where and why Itô’s lemma is applied. Below are some examples.

f(Wₜ)  = Wₜ²               → original function
f′(Wₜ) = 2Wₜ,  f″(Wₜ) = 2   → derivatives

df = f′(Wₜ)dWₜ + ½f″(Wₜ)(dWₜ)²
   = 2WₜdWₜ + (dWₜ)²
   = 2WₜdWₜ + dt      → Itô's lemma version of 2nd order Taylor expansion

df = f′(Wₜ)dWₜ + ½f″(Wₜ)(dWₜ)²
   = 2WₜdWₜ + (dWₜ)²
   = 2WₜdWₜ           → ordinary version of 2nd order Taylor expansion, 2nd order dropped

As we can see, the two results are different. And only Itô’s lemma’s version holds true to the Brownian motion rules upon integration with expectations.

Ordinary version (does not match Brownian motion requirements):

E[W_T²] = 2E[∫₀ᵀ Wₜ dWₜ]
E[∫₀ᵀ Wₜ dWₜ] = 0
E[W_T²] = 0

Itô’s lemma version (meets Brownian motion requirements):

E[W_T²] = 2E[∫₀ᵀ Wₜ dWₜ] + T
E[∫₀ᵀ Wₜ dWₜ] = 0
E[W_T²] = T

Geometric Brownian motion and log normal returns

The increment of a Brownian motion can be defined like so:

dWₜ = W(t + dt) − W(t) ∼ N(0, dt)

dWₜ = increment of Brownian motion at time t
W(t + dt) − W(t) ∼ N(0, dt) = the difference between the value at t + dt
    and the value at t, drawn from a normal distribution with mean 0
    and variance dt

It can be simplified to:

dWₜ = √dt · Z,   Z ∼ N(0, 1)

Suppose we want to model a stock price Sₜ. We can use the following formula for percentage change:

dSₜ / Sₜ = μdt + σdWₜ

μdt  = average return component
σdWₜ = random Brownian path component
σ    = scale of fluctuation magnitude

This assumes the final magnitude of μdt and σdWₜ scales with the current price (geometric Brownian motion).

Say we want to get the log of this series. We don’t use Itô’s lemma correction, and we want the log for future arbitrary purposes.

Ordinary calculus formula:

d(ln S) = (1/S)dS

Apply it to our percent change equation:

d(ln S) = μdt + σdW

Which integrates to (to prove that the calculus is wrong):

ln Sₜ = ln S₀ + μt + σWₜ

And exponentiating gives:

Sₜ = S₀ · e^(μt + σWₜ)

Remember this equation, which we will use later on to prove it’s incorrect. Now suppose the drift is zero, μ = 0. The original stock price model is dSₜ = σSₜdWₜ. Since there is no drift and the Brownian random component’s increments have mean 0, the model requires:

E[dSₜ | Sₜ] = 0

The expected value of price changes is zero. However, when using Sₜ = S₀·e^(μt + σWₜ) to predict E[dSₜ | Sₜ], it is positive. It is not zero, which contradicts the requirements, therefore it is wrong.

The correct version with Itô’s lemma applied is:

d(ln S) = (μ − ½σ²)dt + σdWₜ   → notice the ½σ² correction term

ln S_T = ln S₀ + (μ − ½σ²)T + σW_T

With the logarithm, the Brownian motion at time T follows a normal distribution, W_T ∼ N(0, T). Multiplying by σ gives σW_T ∼ N(0, σ²T).

Adding a constant changes the mean but not the variance, validating that geometric Brownian motion changes are log normal. However, arithmetic Brownian motion changes aren’t log normal because they do not scale with price.

By solving the logarithm by exponentiating, we get S_T = S₀·e^((μ − ½σ²)T + σW_T). We can interpret that solution by calculating the mean and median:

E[S_T]      = S₀ · e^(μT)
Median(S_T) = S₀ · e^((μ − ½σ²)T)

Notice that they differ by a term: −½σ². It is the difference in growth rate between the mean and median price, which gives us a skewed distribution for returns (log normal).

Randomness hedging to Black Scholes PDE

If we combine an option with the right number of shares, we can eliminate stock price randomness. Once the portfolio is riskless, no-arbitrage tells us how much it should earn, which gives us the Black Scholes PDE:

∂V/∂t + ½σ²S²(∂²V/∂S²) + rS(∂V/∂S) − rV = 0

To understand that, we can first use the simpler binomial model as an example. Suppose:

S = 100,  S_u = 120,  S_d = 90,  V_u = 30,  V_d = 5

S = stock price
V = option value

And we have a portfolio like so:

Π = V − ΔS

If we don’t have the stock hedge ΔS and only hold the option V, our portfolio value will be:

Π_u = V_u = 30
Π_d = V_d = 5

In the up and down states. However, the final portfolio is uncertain, as the up and down states have different values. That’s why the ΔS term is there. We can solve for Δ using S_u and S_d

Π_u = V_u − ΔS_u
Π_d = V_d − ΔS_d
Π_u = Π_d

Δ = (V_u − V_d) / (S_u − S_d) = 5/6

Verified by:

Π_u = 30 − (5/6)(120) = 30 − 100 = −70
Π_d =  5 − (5/6)(90)  =  5 − 75  = −70
Π_u = Π_d

For the calculus (continuous) version, it is a similar concept. We start from the stock price model from geometric Brownian motion:

dS = μSdt + σSdW

As well as the differential of option value relative to the stock price:

dV = (∂V/∂t)dt + (∂V/∂S)dS

Above is the chain rule version without Itô’s lemma. Once it is applied, we get:

dV = (V_t + μSV_S + ½σ²S²V_SS)dt + σSV_S dW

V_t  = ∂V/∂t
V_S  = ∂V/∂S
V_SS = ∂²V/∂S²

We then construct a portfolio similar to the binomial tree model:

Π = V − ΔS

Over an interval, the calculus version (instantaneous change) is:

dΠ = dV − ΔdS

Substitute the equations from earlier and group:

dΠ = (V_t + μSV_S + ½σ²S²V_SS − ΔμS)dt + σS(V_S − Δ)dW

In this scenario, instead of a certain (hedged) portfolio value remaining identical across two states, the change in the portfolio depends only on time. There will be no random influence (dW) now. To do so, choose Δ so that:

Δ = V_S  ⟹  dΠ = (V_t + ½σ²S²V_SS)dt

A delta hedged portfolio is constructed. Under Black Scholes assumptions, a riskless portfolio must earn the risk free rate r for there to be no arbitrage opportunities. Thus:

dΠ = rΠdt,   Π = V − ΔS,   dΠ = r(V − SV_S)dt

We can combine the two equations:

dΠ = (V_t + ½σ²S²V_SS)dt
dΠ = r(V − SV_S)dt

which gives us the Black Scholes PDE:

V_t + ½σ²S²V_SS = rV − rSV_S

∂V/∂t + ½σ²S²(∂²V/∂S²) + rS(∂V/∂S) − rV = 0

Thank you for reading Issue 33 of Balaena Quant Insights. We look forward to sharing more of our research and perspectives with you in the next issue.