1. LEARNING OBJECTIVES

By the end of this lesson, you will be able to:

• Define Stochastic Differential Equations (SDEs) and distinguish between drift and diffusion terms.

• Classify SDEs as linear, non-linear, autonomous, and non-autonomous.

• Solve linear SDEs using the integrating factor method and Ito’s Lemma.

• Derive the explicit solution to the Ornstein-Uhlenbeck (OU) process.

• Derive the explicit solution to Geometric Brownian Motion (GBM).

• Solve the Cox-Ingersoll-Ross (CIR) process using the square root transformation.

• Understand the Vasicek model for interest rates and its properties.

• Apply the Feynman-Kac formula to solve PDEs using SDEs.

• Understand the relationship between SDEs and partial differential equations.

2. STOCHASTIC DIFFERENTIAL EQUATIONS – DEFINITION AND CLASSIFICATION

2.1 Definition

A Stochastic Differential Equation (SDE) is an equation of the form:

dX_t = μ(t, X_t) dt + σ(t, X_t) dW_t

Where:

• μ(t, X_t) is the drift coefficient (deterministic part).

• σ(t, X_t) is the diffusion coefficient (stochastic part).

• W_t is a Brownian motion.

2.2 Classification of SDEs

By Linearity:

• Linear SDE: The coefficients μ and σ are linear in X_t:

  dX_t = (a(t) X_t + b(t)) dt + (c(t) X_t + d(t)) dW_t

• Non-Linear SDE: The coefficients are non-linear in X_t (e.g., X_t², sin(X_t), e^{X_t}).

By Time Dependence:

• Autonomous (Time-Homogeneous): μ and σ do not depend explicitly on t:

  dX_t = μ(X_t) dt + σ(X_t) dW_t

• Non-Autonomous (Time-Inhomogeneous): μ and σ depend on t:

  dX_t = μ(t, X_t) dt + σ(t, X_t) dW_t

2.3 Existence and Uniqueness Conditions

For an SDE to have a unique solution, the coefficients must satisfy:

1. Lipschitz Condition: There exists a constant K such that for all x, y:

   |μ(t, x) – μ(t, y)| + |σ(t, x) – σ(t, y)| ≤ K |x – y|

2. Linear Growth Condition: There exists a constant K such that:

   |μ(t, x)| + |σ(t, x)| ≤ K (1 + |x|)

These conditions ensure that the SDE has a unique strong solution.

3. SOLVING LINEAR SDES – THE INTEGRATING FACTOR METHOD

3.1 The General Linear SDE

The general linear SDE is:

dX_t = (a(t) X_t + b(t)) dt + (c(t) X_t + d(t)) dW_t

3.2 The Homogeneous Solution

First, solve the homogeneous equation (b(t) = 0, d(t) = 0):

dX_t = a(t) X_t dt + c(t) X_t dW_t

This is a linear SDE with multiplicative noise. The solution is:

X_t^h = X_0 * exp( ∫_0^t (a(s) – (1/2) c(s)²) ds + ∫_0^t c(s) dW_s )

3.3 The Particular Solution

Using the variation of constants method, the general solution is:

X_t = Φ_t * ( X_0 + ∫_0^t Φ_s^{-1} (b(s) – c(s) d(s)) ds + ∫_0^t Φ_s^{-1} d(s) dW_s )

Where Φ_t is the fundamental solution:

Φ_t = exp( ∫_0^t (a(s) – (1/2) c(s)²) ds + ∫_0^t c(s) dW_s )

3.4 Application to GBM

For GBM, a(t) = μ, b(t) = 0, c(t) = σ, d(t) = 0.

The fundamental solution is:

Φ_t = exp( (μ – (1/2) σ²) t + σ W_t )

The general solution is:

X_t = X_0 * exp( (μ – (1/2) σ²) t + σ W_t )

This matches the known solution for GBM.

4. THE ORNSTEIN-UHLENBECK (OU) PROCESS

4.1 Definition

The OU process is a mean-reverting process defined by the SDE:

dX_t = θ (μ – X_t) dt + σ dW_t

Where:

• θ > 0 is the speed of mean reversion.

• μ is the long-term mean level.

• σ is the volatility.

4.2 Properties of the OU Process

1. Mean Reversion: X_t tends to drift toward μ.

2. Stationary Distribution: As t → ∞, X_t ~ N(μ, σ² / (2θ)).

3. Autocorrelation: For s ≤ t:

   Corr(X_s, X_t) = e^{-θ(t – s)}

4.3 Solving the OU Process

Rewrite the SDE as:

dX_t + θ X_t dt = θ μ dt + σ dW_t

Multiply by the integrating factor e^{θt}:

e^{θt} dX_t + θ e^{θt} X_t dt = θ μ e^{θt} dt + σ e^{θt} dW_t

The left side is d(e^{θt} X_t):

d(e^{θt} X_t) = θ μ e^{θt} dt + σ e^{θt} dW_t

Integrating from 0 to t:

e^{θt} X_t = X_0 + θ μ ∫_0^t e^{θs} ds + σ ∫_0^t e^{θs} dW_s

Solve for X_t:

X_t = e^{-θt} X_0 + μ (1 – e^{-θt}) + σ e^{-θt} ∫_0^t e^{θs} dW_s

4.4 Mean and Variance of the OU Process

Mean:

E[X_t] = e^{-θt} X_0 + μ (1 – e^{-θt})

Variance:

Var(X_t) = (σ² / (2θ)) (1 – e^{-2θt})

4.5 Financial Application – Interest Rate Modelling

The OU process is used as the Vasicek model for interest rates:

dr_t = θ (μ – r_t) dt + σ dW_t

Where r_t is the short rate. The Vasicek model allows for positive and negative interest rates.

5. THE COX-INGERSOLL-ROSS (CIR) PROCESS

5.1 Definition

The CIR process is a mean-reverting process with a square-root diffusion term:

dX_t = θ (μ – X_t) dt + σ √X_t dW_t

Where θ > 0, μ > 0, σ > 0.

5.2 Properties of the CIR Process

1. Non-Negativity: If X_0 > 0 and 2θμ > σ² (the Feller condition), then X_t > 0 for all t.

2. Stationary Distribution: As t → ∞, X_t ~ Gamma(2θμ/σ², 2θ/σ²).

3. Mean Reversion: X_t tends to drift toward μ.

5.3 Solving the CIR Process

The CIR process does not have a closed-form explicit solution like GBM or OU. However, we can solve it approximately or use the square root transformation.

Let Y_t = √X_t. Applying Ito’s Lemma:

dY_t = ( (1/(2√X_t)) θ (μ – X_t) – (1/8) σ² X_t^{-3/2} ) dt + (σ/2) dW_t
dY_t = ( (1/(2X_t^{1/2})) θ (μ – X_t) – (1/8) σ² X_t^{-3/2} ) dt + (σ/2) dW_t

This is not linear, but the CIR process is widely used in finance for interest rates and volatility modelling.

5.4 Financial Application – Interest Rate Modelling

The CIR process is used as an alternative to the Vasicek model:

dr_t = θ (μ – r_t) dt + σ √r_t dW_t

The CIR model ensures positive interest rates (if the Feller condition is satisfied).

6. THE VASICEK MODEL

6.1 Definition

The Vasicek model is an OU process for the short rate:

dr_t = θ (μ – r_t) dt + σ dW_t

6.2 Solution to the Vasicek Model

From the OU process solution:

r_t = e^{-θt} r_0 + μ (1 – e^{-θt}) + σ e^{-θt} ∫_0^t e^{θs} dW_s

6.3 Distribution of r_t

The short rate is normally distributed:

r_t ~ N( e^{-θt} r_0 + μ (1 – e^{-θt}), (σ²/(2θ)) (1 – e^{-2θt}) )

6.4 Bond Pricing in the Vasicek Model

The price of a zero-coupon bond in the Vasicek model is:

P(t, T) = A(t, T) * e^{-B(t, T) r_t}

Where:

B(t, T) = (1 – e^{-θ(T-t)}) / θ
A(t, T) = exp( (B(t, T) – (T-t)) (θ² μ – σ²/2) / θ² – σ² B(t, T)² / (4θ) )

6.5 Yield Curve

The continuously compounded yield is:

Y(t, T) = – (ln(P(t, T))) / (T – t)

The Vasicek model generates upward-sloping, downward-sloping, and humped yield curves.

6.6 Limitations of the Vasicek Model

1. Negative Interest Rates: The normal distribution allows for negative interest rates.

2. Constant Volatility: The volatility does not depend on the level of interest rates.

3. No Jumps: The model does not capture sudden changes in interest rates.

7. THE FELLER CONDITION AND THE SQUARE ROOT PROCESS

7.1 The Feller Condition

For the CIR process, the Feller condition is:

2θμ > σ²

If this condition holds, the process is strictly positive (X_t > 0 for all t). If the condition fails, the process can hit zero.

7.2 The Square Root Process

The CIR process is also known as the square root process. It is the foundation of the Heston model for stochastic volatility.

7.3 The Heston Model

The Heston model extends the Black-Scholes model by allowing volatility to be stochastic:

dS_t = r S_t dt + √v_t S_t dW_t^{(1)}
dv_t = θ (μ – v_t) dt + σ √v_t dW_t^{(2)}
dW_t^{(1)} dW_t^{(2)} = ρ dt

Where:

• v_t is the variance (volatility squared).

• θ is the speed of mean reversion.

• μ is the long-term variance.

• σ is the volatility of variance (vol-of-vol).

• ρ is the correlation between asset returns and variance.

8. THE FEYNMAN-KAC FORMULA

8.1 The Theorem

The Feynman-Kac formula relates SDEs to PDEs. Consider the SDE:

dX_t = μ(t, X_t) dt + σ(t, X_t) dW_t

Let f(t, x) be defined as:

f(t, x) = E[ g(X_T) | X_t = x ]

Then f satisfies the PDE:

∂f/∂t + μ(t, x) ∂f/∂x + (1/2) σ(t, x)² ∂²f/∂x² = 0

With the terminal condition:

f(T, x) = g(x)

8.2 Proof Sketch

Using Ito’s Lemma on f(t, X_t):

df = (∂f/∂t + μ ∂f/∂x + (1/2) σ² ∂²f/∂x²) dt + σ ∂f/∂x dW_t

Taking expectations:

E[df] = (∂f/∂t + μ ∂f/∂x + (1/2) σ² ∂²f/∂x²) dt

Since f is a martingale (E[df] = 0), we get the PDE.

8.3 Financial Application – Option Pricing

The Feynman-Kac formula is the foundation of the risk-neutral pricing framework. The price of a European option V(t, S_t) is:

V(t, S_t) = E^Q[ e^{-r(T-t)} g(S_T) | S_t = S ]

This satisfies the Black-Scholes PDE.

9. SUMMARY OF KEY SDES IN FINANCE

10. PRACTICAL IMPLEMENTATION

A. Simulating the OU Process:

import numpy as np
import matplotlib.pyplot as plt

def simulate_OU(theta, mu, sigma, X0, T, n_steps, n_paths):
    dt = T / n_steps
    dW = np.random.normal(0, np.sqrt(dt), (n_paths, n_steps))
    X = np.zeros((n_paths, n_steps + 1))
    X[:, 0] = X0
    for i in range(n_steps):
        X[:, i+1] = X[:, i] + theta * (mu - X[:, i]) * dt + sigma * dW[:, i]
    return X

# Parameters
theta = 2.0
mu = 0.05
sigma = 0.02
X0 = 0.03
T = 1.0
n_steps = 252
n_paths = 100

# Simulate
X = simulate_OU(theta, mu, sigma, X0, T, n_steps, n_paths)
t = np.linspace(0, T, n_steps + 1)

# Plot
plt.figure(figsize=(10, 6))
for i in range(min(10, n_paths)):
    plt.plot(t, X[i, :], alpha=0.6)
plt.xlabel('Time')
plt.ylabel('X_t')
plt.title('Simulated OU Process Paths')
plt.axhline(y=mu, color='r', linestyle='--', label='Long-term mean μ')
plt.legend()
plt.show()

# Verify stationary distribution
print(f"Mean of X_T: {np.mean(X[:, -1]):.4f}")
print(f"Expected Mean: {mu:.4f}")
print(f"Variance of X_T: {np.var(X[:, -1]):.4f}")
print(f"Expected Variance: {sigma**2 / (2 * theta):.4f}")

B. Simulating the CIR Process:

def simulate_CIR(theta, mu, sigma, X0, T, n_steps, n_paths):
    dt = T / n_steps
    dW = np.random.normal(0, np.sqrt(dt), (n_paths, n_steps))
    X = np.zeros((n_paths, n_steps + 1))
    X[:, 0] = X0
    for i in range(n_steps):
        # Ensure non-negativity
        X[:, i+1] = X[:, i] + theta * (mu - X[:, i]) * dt + sigma * np.sqrt(np.maximum(X[:, i], 0)) * dW[:, i]
        X[:, i+1] = np.maximum(X[:, i+1], 0)  # Reflect at zero
    return X

# Parameters
theta = 2.0
mu = 0.05
sigma = 0.10
X0 = 0.03
T = 1.0
n_steps = 252
n_paths = 100

# Simulate
X = simulate_CIR(theta, mu, sigma, X0, T, n_steps, n_paths)

# Plot
plt.figure(figsize=(10, 6))
for i in range(min(10, n_paths)):
    plt.plot(t, X[i, :], alpha=0.6)
plt.xlabel('Time')
plt.ylabel('X_t')
plt.title('Simulated CIR Process Paths')
plt.axhline(y=mu, color='r', linestyle='--', label='Long-term mean μ')
plt.legend()
plt.show()

# Check Feller condition
feller_condition = 2 * theta * mu / (sigma**2)
print(f"Feller Condition (2θμ/σ²): {feller_condition:.4f}")
print(f"Condition holds: {feller_condition > 1}")

C. Vasicek Bond Pricing:

def vasicek_bond_price(r, t, T, theta, mu, sigma):
    tau = T - t
    B = (1 - np.exp(-theta * tau)) / theta
    A = np.exp((B - tau) * (theta**2 * mu - sigma**2/2) / theta**2 - sigma**2 * B**2 / (4 * theta))
    return A * np.exp(-B * r)

# Parameters
r = 0.03
t = 0
T = 5
theta = 2.0
mu = 0.05
sigma = 0.02

P = vasicek_bond_price(r, t, T, theta, mu, sigma)
Y = -np.log(P) / T
print(f"Bond Price: {P:.4f}")
print(f"Yield: {Y:.4f}")

D. Feynman-Kac Verification:

def feynman_kac_solution(x, t, T, mu, sigma, g):
    # This is the solution to the PDE using the Feynman-Kac formula
    # For a simple case: μ = 0, σ = 1, g(x) = x^2
    # The solution is f(t, x) = x^2 + (T - t)
    return x**2 + (T - t)

# Verify for a simple SDE: dX_t = dW_t
def g(x):
    return x**2

x = 1.0
t = 0.5
T = 1.0
mu = 0.0
sigma = 1.0

f_value = feynman_kac_solution(x, t, T, mu, sigma, g)
print(f"Feynman-Kac solution: {f_value:.4f}")

# Monte Carlo verification
n_paths = 10000
n_steps = 1000
dt = (T - t) / n_steps
W = np.random.normal(0, np.sqrt(dt), (n_paths, n_steps))
X_T = x + np.sum(W, axis=1)
expected = np.mean(g(X_T))
print(f"Monte Carlo expectation: {expected:.4f}")