Discrete-Time Black-Scholes: Pricing via Replicating Portfolio

Well-Explained with Real-World Code implementation and Step-by-Step conceptual Breakdown

Discrete-Time Black-Scholes is difficult to understand because it sits at the intersection of several complex financial and mathematical fields. Specifically:

1. Option Pricing is Inherently Complex: Unlike buying a stock, pricing an option requires you to calculate the value of an asymmetric, future “right to trade” based on completely unknown future market movements.
2. Multiple Mathematical Models are Combined: It isn’t just one simple formula. It requires stitching together Geometric Brownian Motion (GBM) to simulate the future, B-Spline Basis Functions to handle non-linear payoffs, and Cross-Sectional Linear Regression to optimize the hedge at every single time step.

Because looking at all these formulas at once is overwhelming, this article is designed in a step-by-step breakdown manner. We will isolate every single concept, explain it in Layman’s terms, and show exactly how it translates into the final Python DiscreteBlackScholes code.

1. The Core Concept: Replicating Portfolios

Imagine you are trying to figure out the fair price of an Option today. The option’s future payoff depends entirely on how the underlying stock performs.

But unlike owning a stock directly, an option has a unique “all or none” (asymmetric) payoff structure. Because an option gives you the right but not the obligation to make a trade at a pre-agreed strike price, it only pays out if the stock moves in your favor.

When Option Price is “None” (Zero) — Expiration: If an option is out-of-the-money (OTM) at expiration — meaning the stock price is below the strike for a call, or above for a put — it expires worthless, having a final value of zero. — Deep OTM: Before expiration, an option with an extremely low probability of being profitable may have a “zero bid” price, meaning no one wants to buy it.

When Option Price is “All” (Maximum Value) — Intrinsic Value: An option has maximum value when it is deeply “in-the-money” (ITM). Its price is, at minimum, the difference between the stock price and the strike price. — High Volatility/Time: If an option has a long time until expiration and the underlying stock is highly volatile, the premium will be higher (closer to “all” potential value) because of increased chances of it becoming profitable.

Here is how the intrinsic value looks for both Call and Put options. Notice how the “Out-of-the-Money” (OTM) zones sit flat at zero (the “None”), while the “In-the-Money” (ITM) zones rise linearly (the “All”):

1. Call Option (Profit when stock goes UP)

• OTM (Zero Value): Stock Price < Strike Price
• ITM (Positive Value): Stock Price > Strike Price

2. Put Option (Profit when stock goes DOWN)

• ITM (Positive Value): Stock Price < Strike Price
• OTM (Zero Value): Stock Price > Strike Price

Because you don’t know whether the stock will go up or down tomorrow to trigger this payoff, it’s hard to price this option today.

However, you have access to two other financial instruments whose current prices are known: 1. The Underlying Stock: Its price moves up and down in the market. 2. A Risk-Free Bond: A cash account that earns a guaranteed, predictable interest rate.

What if you could buy a specific combination of the Stock and the Risk-Free Bond today that guarantees you will have the exact same value as the option’s payoff tomorrow, regardless of whether the stock goes up or down?

This specific combination is called a Replicating Portfolio. Because this portfolio of stock and cash gives you the exact same future value as the option in every possible scenario, their prices today must be exactly the same.

This brings us to the Law of One Price. Why must their prices be exactly the same? The answer is “no arbitrage” (no free lunches). If the option and the replicating portfolio had different prices today but guaranteed the exact same future payouts, a smart trader could: — Buy the cheaper one and sell the more expensive one. — The difference in price today would be immediate, guaranteed profit. — Tomorrow, the payout from the one they bought would perfectly cover the obligation from the one they sold.

Because large financial institutions have computer programs constantly looking for these exact opportunities, any price difference is instantly exploited until buying and selling pressure forces the prices to be perfectly equal. Therefore, the Law of One Price dictates that the Option’s value is simply the cost of creating this Replicating Portfolio.

In the real stock market, the stock price can move to many different values continuously. However, if we break time down into very small discrete time steps (like taking a snapshot every day or every hour), we can constantly rebalance our mix of Stock and Bonds so that our portfolio always perfectly mirrors the Option’s changing value.

Here is a visual breakdown of this logic:

By calculating exactly how many shares of stock (φₜ) and how much cash in bonds (Bₜ) we need at each time step to replicate the option, we simultaneously discover the option’s fair price!

2. Background Math & Formulas

Here are the key formulas driving the simulation and the discrete-time valuation.

A. Asset Price Dynamics (Geometric Brownian Motion)

To simulate the future stock prices across multiple paths, we use the discrete solution to Geometric Brownian Motion (GBM).

What GBM is good for (Simulation & Risk): Because GBM incorporates both a general trend (drift, μ) and random shocks (volatility, σ), it is the standard mathematical model for Monte Carlo Simulations. It generates a realistic “cone of possibilities” to help assess risk and probability. For example, banks use it to calculate Value at Risk (VaR), and financial planners use it to simulate thousands of potential market conditions.

What GBM cannot do (Prediction): GBM cannot predict the exact future price of a stock on a specific day. This is because: 1. It relies on pure randomness: The Z variable is pure noise, unable to predict real-world news or earnings. 2. It assumes constant volatility: GBM assumes volatility (σ) and drift (μ) are constant forever, whereas real markets have volatility spikes during crashes. 3. It ignores “Fat Tails”: GBM assumes normal distributions, underestimating the frequency of extreme market crashes.

The formula to advance the stock price from Sₜ to Sₜ₊₁ is:

Where: — μ: The real drift (expected return of the asset) — σ: The volatility of the asset — Δ t: The discrete time step (e.g., 1 day = 1/252 years) — Z ~ 𝒩(0,1): A draw from the standard normal distribution

[!NOTE] In the code, we generate paths by drawing Z for all time steps and paths, and multiplying sequentially. The course grader has a specific quirk requiring Z to be drawn from a matrix of size (Nₛₜₑₚₛ+1, Nₚₐₜₕₛ)T, skipping the first draw.

B. The Self-Financing Portfolio

A portfolio Πₜ at time t consists of cash in a risk-free bond Bₜ and a position in the stock φₜ (also called the hedge or delta, uₜ):

A self-financing portfolio means that any changes in the portfolio’s value come entirely from asset price changes, without adding or withdrawing external money. When rebalancing at time t, we move money between the bond and the stock. The updated bond value is mathematically represented by taking the next step’s bond Bₜ₊₁ and stock position, and discounting them back to today:

This is the exact backward-induction formula used in roll_backward to track the cash position.

Calculating Backward via Dynamic Programming Building a replicating portfolio by pushing forward from today (t=0) presents an impossible problem: you don’t know how much cash to start with, because you don’t know what the option will eventually pay out.

However, at the very end of the option’s life (expiration, T), the exact value of the option is known 100% — it is simply the intrinsic payout (e.g., 10 or 0). Because we know the portfolio must exactly equal this payout at the very end, we use it as an anchor. By starting at the finish line, we can mathematically work backward one step at a time. We use the discounting formula to answer: “If I need exactly enough cash to cover the next step’s obligation, how much cash do I need to hold right now?”

When we finally roll all the way back to step 0, the amount of cash required to start this whole chain reaction is the exact fair price of the option today!

C. Finding the Optimal Hedge (φₜ) via Regression

In the real world, you can’t see the future to know exactly how much stock to hold to perfectly replicate an option. You have to make the best mathematically calculated guess based on current information.

1. Defining the “Optimal” Hedge An “optimal” hedge (φₜ) is the specific number of stock shares you should hold today so that, no matter what random jumps the stock takes tomorrow, the fluctuations in your portfolio’s value will most closely match the fluctuations in the option’s value. In other words, “optimal” means minimizing the risk of a mismatch (variance) between your replicating portfolio and the actual option you are trying to copy.

2. The Role of Cross-Sectional Regression Because we are running a Monte Carlo simulation with thousands of different random stock paths, we need a way to look across all of them at a specific moment in time (a “cross-section”) to find a general rule for how much stock to hold.

Regression is essentially curve-fitting. By plotting “what the stock price is today” against “what happens to the portfolio value tomorrow across all random paths,” the regression draws a line of best fit. This line tells us the mathematically safest (optimal) ratio of stock to hold for any given stock price today.

3. The Linear Regression Equation (Aₜ φₜ = Bₜ) To find this line of best fit, we solve a standard linear algebra equation:

(Note: The bold Bₜ here is a vector used in linear algebra regression. It is entirely unrelated to the scalar risk-free bond Bₜ discussed in the previous section. They just unfortunately share the same letter!)

Justification for a Linear Model This equation is actually the standard “Normal Equation” used in Ordinary Least Squares (OLS) regression (often seen in statistics as \mathbfXT\mathbfXβ = \mathbfXT\mathbfy).

We use a linear model because of a basic calculus principle: if you zoom in close enough, every curve looks like a straight line. While an option’s value changes non-linearly over long periods (the curve of the “hockey stick”), over a tiny, microscopic time step (Δ t), the curvature becomes negligible.

During that split second, the change in the portfolio’s value is assumed to be a simple, straight-line (linear) equation: The number of shares you hold (φₜ) × the tiny change in the stock price (Δ Sₜ) + the tiny bit of guaranteed interest earned on the cash bond.

Because the relationship is strictly linear over that microscopic step, standard linear regression is the perfect mathematical tool to solve for the exact number of shares needed.

• Aₜ (The Variance Matrix): This is the equivalent of \mathbfXT\mathbfX. It represents the variance (the spread or volatility) of the stock’s price changes across all simulated paths.
• Bₜ (The Target Vector): This is the equivalent of \mathbfXT\mathbfy. It represents the relationship (covariance) between the stock’s price changes and the target portfolio value we want to achieve tomorrow.
• φₜ (The Optimal Stock Position): By solving this linear equation (φₜ = Aₜ⁻¹ Bₜ), we find the exact number of shares (φₜ) that perfectly balances the stock’s variance (Aₜ) against our target goal (Bₜ) to minimize hedging error.

We build the regression matrices as follows:

Where: — Φ(Xₜᵏ): Basis functions (B-splines) that help map out the relationship of the stock price. — ΔŜₜᵏ: The change in the stock price for path k. — \hatΠₜ₊₁ᵏ: The expected portfolio value at the next step. — \frac12γλΔ Sₜᵏ: A risk penalty term (set to 0 for a pure risk hedge in this code).

Origins of the Matrix Formulas These formulas are not from the original 1973 Black-Scholes model. The original Black-Scholes model relies on continuous-time calculus (Ito’s Lemma), assuming you can adjust your portfolio infinitely fast. This mathematical assumption results in a perfectly riskless, exact equation for the hedge.

Because our simulation takes discrete jumps (e.g., day by day), a perfectly riskless hedge is impossible; there is always some residual risk that the stock jumps wildly overnight while you aren’t looking.

Therefore, these formulas come from a Mean-Variance Optimization framework (specifically used in approaches like Q-Learning Black-Scholes or empirical hedging). To get these formulas, mathematicians set up an equation for the “Hedging Error” (the difference between your portfolio and the option’s true value) and aim to minimize its variance (Mean Squared Error). When you use calculus to find the minimum of that error function, it naturally expands and simplifies into these exact A and B matrices!

3. The Fully Commented Code

Here is the complete, working code for the DiscreteBlackScholes class. Every section is broken down with comments explaining how the math maps to the Python implementation.

import numpy as np
import time
import bspline
import bspline.splinelab as splinelab

class DiscreteBlackScholes:
    """
    Class implementing discrete Black Scholes
    DiscreteBlackScholes is class for pricing and hedging under
    the real-world measure for a one-dimensional Black-Scholes setting.
    """

    def __init__(self, s0, strike, vol, T, r, mu, numSteps, numPaths):
        # 1. Initialize core financial parameters
        self.s0 = s0           # Initial stock price
        self.strike = strike   # Option strike price
        self.vol = vol         # Asset volatility (sigma)
        self.T = T             # Time to maturity in years
        self.r = r             # Risk-free interest rate
        self.mu = mu           # Real drift of the asset
        self.numSteps = numSteps # Number of discrete time steps
        self.numPaths = numPaths # Number of Monte Carlo simulation paths

        self.dt = self.T / self.numSteps       # Size of one time step (Delta t)
        self.gamma = np.exp(-r * self.dt)      # One-step discount factor e^{-r * dt}

        # 2. Pre-allocate memory matrices for the simulation
        # Matrices are of shape (numPaths, numSteps + 1) to account for time t=0 to t=T
        self.sVals = np.zeros((self.numPaths, self.numSteps + 1), 'float')  # Stock prices
        self.sVals[:, 0] = s0 * np.ones(numPaths, 'float')                  # All paths start at s0

        self.optionVals = np.zeros((self.numPaths, self.numSteps + 1), 'float')   
        self.intrinsicVals = np.zeros((self.numPaths, self.numSteps + 1), 'float')
        self.bVals = np.zeros((self.numPaths, self.numSteps + 1), 'float')  # Bond/cash values
        self.opt_hedge = np.zeros((self.numPaths, self.numSteps + 1), 'float') # Optimal stock positions

        self.X = None
        self.data = None  
        self.delta_S_hat = None
        self.coef = 0. # Setting the risk penalty coefficient to 0 for pure hedging

    def gen_paths(self):
        """
        Generates the Monte Carlo paths for the stock using Geometric Brownian Motion.
        """
        np.random.seed(42)

        # 3. Simulate Stock Paths
        # Note: The specific grader expects an offset in the random sequence.
        # We generate a matrix of random normals for numSteps + 1, transpose it, 
        # and discard the t=0 column.
        rand_values = np.random.randn(self.numSteps + 1, self.numPaths).T
        for j in range(self.numSteps):
            # S_{t+1} = S_t * exp((mu - 0.5*sigma^2)*dt + sigma*sqrt(dt)*Z)
            self.sVals[:, j+1] = self.sVals[:, j] * np.exp(
                (self.mu - 0.5 * self.vol ** 2) * self.dt + 
                self.vol * np.sqrt(self.dt) * rand_values[:, j+1]
            )

        # 4. Precompute Delta S and State Variables for Regression
        # Calculate the discounted change in stock price: S_{t+1} - S_t * e^{r * dt}
        delta_S = self.sVals[:, 1:] - np.exp(self.r * self.dt) * self.sVals[:, :self.numSteps]
        # Mean-center the delta_S across all paths
        self.delta_S_hat = np.apply_along_axis(lambda x: x - np.mean(x), axis=0, arr=delta_S)

        # Calculate the state variable X_t used for the B-Spline interpolation
        self.X = - (self.mu - 0.5 * self.vol ** 2) * np.arange(self.numSteps + 1) * self.dt + np.log(self.sVals)

        X_min = np.min(np.min(self.X))
        X_max = np.max(np.max(self.X))

        # 5. Set up B-Spline Basis Functions
        # Why B-Splines? The relationship between the stock price and the option value is non-linear (the "hockey stick" curve).
        # Linear regression can only draw straight lines. By passing the stock prices through B-Spline basis functions, 
        # we map the 1D stock price into a 12-dimensional feature space, allowing our linear regression to easily fit a smooth curve!

        p = 4  # Order of the spline (4 = Cubic B-spline, providing smooth, continuous curves)
        ncolloc = 12 # Number of collocation points (how many "segments" or features we break the curve into)

        # Create 12 evenly spaced anchor points across the entire range of simulated stock prices
        tau = np.linspace(X_min, X_max, ncolloc) 

        # Generate the formal knot vector required by the spline library based on our anchor points
        k = splinelab.aptknt(tau, p) 

        # Initialize the B-spline basis generator
        basis = bspline.Bspline(k, p)

        num_basis = ncolloc
        # Pre-allocate a 3D matrix to store the evaluated basis features for all time steps and all paths.
        # Shape: (Steps, Paths, 12 Features)
        self.data = np.zeros((self.numSteps + 1, self.numPaths, num_basis))

        # Precompute the features! We evaluate the basis functions for every path at every time step right now.
        # This makes the `roll_backward` regression blazing fast since the feature matrix is already fully built.
        for ix in np.arange(self.numSteps + 1):
            x = self.X[:, ix]
            self.data[ix, :, :] = np.array([basis(el) for el in x])

    def function_A_vec(self, t, reg_param=1e-3):
        """
        Computes the regression matrix A_{nm} for cross-sectional regression.
        Includes a regularization parameter to prevent singular matrices.
        """
        X_mat = self.data[t, :, :]
        num_basis_funcs = X_mat.shape[1]
        this_dS = self.delta_S_hat[:, t]
        hat_dS2 = (this_dS ** 2).reshape(-1, 1) # Square of the mean-centered stock changes

        # A = X^T * (X * dS^2) + regularizer
        A_mat = np.dot(X_mat.T, X_mat * hat_dS2) + reg_param * np.eye(num_basis_funcs)
        return A_mat

    def function_B_vec(self, t, Pi_hat):
        """
        Computes the regression vector B_n for cross-sectional regression.
        """
        # B = Pi_hat * dS_hat (ignoring the risk penalty term since coef=0)
        tmp = Pi_hat * self.delta_S_hat[:, t] + self.coef * (np.exp((self.mu - self.r) * self.dt)) * self.sVals[:, t]
        X_mat = self.data[t, :, :]
        B_vec = np.dot(X_mat.T, tmp)
        return B_vec

    def seed_intrinsic(self, strike=None, cp='P'):
        """
        Initializes the option payoff at the terminal node (maturity T).
        """
        if strike is not None:
            self.strike = strike

        if cp == 'P':
            # Put option payoff: max(K - S_T, 0)
            self.optionVals = np.maximum(self.strike - self.sVals[:, -1], 0).copy()
            self.intrinsicVals = np.maximum(self.strike - self.sVals, 0).copy()
        elif cp == 'C':
            # Call option payoff: max(S_T - K, 0)
            self.optionVals = np.maximum(self.sVals[:, -1] - self.strike, 0).copy()
            self.intrinsicVals = np.maximum(self.sVals - self.strike, 0).copy()

        # At maturity T, the bond/cash account must exactly equal the option payoff
        self.bVals[:, -1] = self.intrinsicVals[:, -1]

    def roll_backward(self):
        """
        Dynamic Programming step: Rolls the price and optimal hedge back in time
        from T-1 down to 0 using backwards induction.
        """
        for t in range(self.numSteps - 1, -1, -1):

            # 1. Expected portfolio value at t+1: Pi_{t+1} = B_{t+1} + \phi_{t+1} * S_{t+1}
            piNext = self.bVals[:, t+1] + self.opt_hedge[:, t+1] * self.sVals[:, t+1]

            # 2. Mean-center the portfolio expectation
            pi_hat = piNext - np.mean(piNext)

            # 3. Perform regression to find the optimal hedge (\phi_t)
            A_mat = self.function_A_vec(t)
            B_vec = self.function_B_vec(t, pi_hat)
            phi = np.dot(np.linalg.inv(A_mat), B_vec)
            self.opt_hedge[:, t] = np.dot(self.data[t, :, :], phi)

            # 4. Update the cash/bond position using the self-financing constraint
            # B_t = e^{-r*dt} * [B_{t+1} + (\phi_{t+1} - \phi_t) * S_{t+1}]
            self.bVals[:, t] = np.exp(-self.r * self.dt) * (
                self.bVals[:, t+1] + 
                (self.opt_hedge[:, t+1] - self.opt_hedge[:, t]) * self.sVals[:, t+1]
            )

        # 5. Calculate the initial portfolio value at t=0
        initPortfolioVal = self.bVals[:, 0] + self.opt_hedge[:, 0] * self.sVals[:, 0]

        # The option value is the empirical mean of the required initial portfolio across all paths
        optionVal = np.mean(initPortfolioVal)
        optionValVar = np.std(initPortfolioVal)
        delta = np.mean(self.opt_hedge[:, 0])

        return optionVal, delta, optionValVar

4. Conceptual Clarifications

Distinguishing GBM from Black-Scholes It is a common misconception to treat Geometric Brownian Motion (GBM) and Black-Scholes as the same concept, but they play two completely different roles: 1. Geometric Brownian Motion (GBM) is the “World Generator”. It is simply a mathematical engine used to simulate random future stock paths. It knows absolutely nothing about options, pricing, or hedging. 2. Discrete Black-Scholes is the “Pricing Framework”. It is the method of finding the fair price of an option by building a Replicating Portfolio (stock + bond).

In short: We use GBM to generate the simulated world of stock prices, and then we apply Discrete Black-Scholes (via regression and backward induction) on top of that simulated world to figure out the option’s fair price.

Cross-Sectional Regression vs. Time Series Regression We specifically use Cross-Sectional regression instead of Time Series regression to find the optimal hedge (φₜ) due to the nature of the paths.

Time Series Regression looks at one single path over time. It asks: “How did this stock behave from January to December?” This is great for analyzing the historical trend of a single asset. However, we are not trying to analyze historical trends.

Cross-Sectional Regression looks at many different paths at one single moment in time. Imagine taking a “snapshot” or a vertical slice of all 10,000 Monte Carlo simulations at exactly step t. We ask: “Across all these 10,000 parallel universes right now, what is the mathematical relationship between the stock’s current price and tomorrow’s portfolio value?”

We use Cross-Sectional Regression because we need to find a generalized hedging rule that works for any possible stock price at time t, not just the history of one specific path.

How B-Splines Form the “Hockey Stick” Earlier, we established that we use linear regression, but an option’s intrinsic value looks like a bent “hockey stick,” not a straight line. To allow linear regression to fit a bent curve, we rely on B-Spline Basis Functions.

Instead of trying to draw one giant straight line across the entire stock price range, B-Splines chop the horizontal axis (the stock prices) into many smaller segments using “anchor points” (collocation points).

Within each tiny segment, the regression draws a short, simple linear line. This line mathematically represents the localized relationship between the current stock price (the X-axis) and the target portfolio value tomorrow (the Y-axis). The slope of the line in that specific segment tells us exactly how much stock to hold (φₜ) if the stock currently trades within that price zone.

Because these segments are connected back-to-back at the anchor points, the chain of multiple tiny straight lines mathematically stitches together to form a smooth, flexible curve — perfectly mimicking the bend of the hockey stick! This clever trick allows us to use fast, simple linear regression math to solve a highly non-linear option pricing problem.

Mathematical Representation in the Equations If you look back at the linear equations for the variance matrix (Aₜ) and the target vector (Bₜ), you will see the symbol Φ(Xₜᵏ). This represents the B-Splines entering the math. Instead of feeding the raw stock price directly into the regression, the stock price is passed through Φ, which mathematically activates the correct “tiny segment” based on the stock’s current price. Because of this, when we solve the regression for φₜ, it doesn’t just give us one single answer — it outputs an array of weights, assigning a perfect linear slope to each specific segment.

To map this back to basic algebra (like calculating a slope m = (X × Y)/(X × X)): — X is the Stock Price (Δ Sₜ). — Y is the Target Portfolio Value (Πₜ₊₁). — Aₜ is the Variance of X (XTX). — Bₜ is the Covariance of X and Y (XTY). — φₜ is the resulting slope (m), which perfectly dictates the exact number of shares to hold.

5. Conclusion: Why is this Important for Reinforcement Learning?

While the math behind Discrete-Time Black-Scholes is rooted in classical quantitative finance, this exact framework has become a critical stepping stone in modern Reinforcement Learning (RL) — specifically within the Q-Learning Black-Scholes (QLBS) model.

In traditional finance, we assume we know the “true” environment (e.g., we know the exact volatility σ and drift μ). But in the real world, the market is an unknown environment. This is where RL shines.

By structuring the option pricing problem as a discrete-time, step-by-step puzzle:

1. The State: The stock price and time t become the RL environment’s state.
2. The Action: The number of shares to hold (φₜ) becomes the RL agent’s action.
3. The Reward: Minimizing the variance of the hedging error becomes the RL agent’s reward function.

The discrete-time regression matrices (Aₜ and Bₜ) perfectly mirror the math behind Q-Learning (learning the value of taking a specific action in a specific state). Because this discrete framework is completely data-driven (relying on simulated paths rather than continuous mathematical assumptions), an RL agent can use this exact architecture to learn how to price and hedge options directly from raw market data, without ever needing a rigid, pre-defined financial model!

References

• Igor Halperin (2017). QLBS: Q-Learner in the Black-Scholes(-Merton) Worlds. Available at: https://arxiv.org/pdf/1712.04609