Option Pricing: PDEs and the Finite Element Method

Derivative prices can often be written as discounted conditional expectations of future payoffs. Feynman–Kac turns that expectation into a parabolic pricing equation. This article develops the finite element method for that equation, first in an abstract setting, then for Black–Scholes, and finally for seven important extensions.

1. From weak formulations to stability and error estimates

Let $J=(0,T]$, let $G=(x_L,x_R)$ be a one-dimensional spatial domain, and let $u=u(t,x)$ be unknown. The functions $\alpha(x)>0$, $\beta(x)$, and $\gamma(x)$ are the diffusion, convection, and reaction coefficients. The right-hand side $f=f(t,x)$ is known and does not depend on $u$. Consider

with an initial condition and suitable boundary conditions.

FDM starts from this strong equation and replaces derivatives by difference quotients. On a uniform grid, for example,

FEM takes a different route. It first rewrites the PDE in weak form and then looks for

in a finite-dimensional function space $V_h\subset V$¹. The basis functions $\phi_j$ describe the spatial shape. Only the coefficients $U_j(t)$ evolve in time. The Galerkin condition makes the discrete residual vanish against every $v_h\in V_h$, which preserves the energy structure of the weak problem and produces the matrix equations below.

1.1 When to use FEM instead of FDM?

For a smooth one-dimensional problem with constant coefficients and a regular domain, FDM is often the simplest and cheapest method. Three recurring features of derivatives pricing make a direct discretization of the strong equation more delicate.

Nonsmooth payoffs. A vanilla payoff has a kink at the strike; a digital payoff is discontinuous. A barrier boundary condition may also be incompatible with the payoff at $\tau=0$. The strong equation requires the classical derivatives appearing in the PDE, and those derivatives need not exist near maturity.

The usual truncation-error argument for FDM has the same issue. A stencil is the set of nearby nodes used to approximate a derivative. The centered approximation of $u_{xx}(x_i)$ uses $x_{i-1},x_i,x_{i+1}$. Proving a local error of $O(h^2)$ by Taylor expansion normally assumes $u\in C^4$ in a neighborhood covered by the stencil. A kink or jump violates that assumption at $\tau=0$.

Weak solution theory allows $u_0\in H=L^2(G)$ and still gives existence and uniqueness under the conditions stated below. This does not mean that FEM can solve problems that FDM cannot. It means that the weak framework remains available without first restoring classical smoothness. FDM may instead use payoff smoothing, local refinement, or a specialized error argument. Low regularity still matters for FEM: it increases the interpolation error and may reduce the observed convergence order.

Irregular boundaries and constraints. The exercise boundary of an American option is unknown before the problem is solved. Barriers, multi-asset domains, and other constraints may not align with a Cartesian grid. FDM can handle these cases through embedded boundaries, ghost points, penalty methods, or projections, but the local stencil and its analysis usually have to be adjusted.

FEM works element by element. Local matrices are computed on each element and assembled through shared degrees of freedom. This makes unstructured meshes and local refinement natural. Replacing an equation by a variational inequality also allows the same spatial space to be used for the LCP generated by an American claim.

More complicated operators. Correlation in a multi-asset model creates mixed derivatives. Jump diffusions and Lévy models introduce nonlocal integrals: the equation at $x$ depends on values far from $x$. A local finite-difference stencil is no longer sufficient and numerical quadrature adds another error component.

FEM must evaluate the same nonlocal integral, so it does not remove the computational cost. Its benefit is structural: the integral can be added to the bilinear form, after which continuity, a Gårding inequality, stability, and convergence can be discussed in the same language. The price is algebraic: a sparse banded stiffness matrix may become dense or block structured.

The reason to use FEM is therefore not that it dominates FDM in every setting. It provides a reusable formulation for problems with weak data, unusual geometry, constraints, or nonlocal operators.

1.2 From the strong equation to the weak equation

Changing from a strong to a weak formulation does not change the pricing model. It changes the sense in which the equation is required to hold. Multiply the strong equation by a test function $v$ and integrate over $G$:

Integrating the diffusion term by parts gives

The second derivative of $u$ has been replaced by one weak derivative of $u$ and one derivative of $v$. A homogeneous Dirichlet condition removes the boundary term by requiring $v=0$ on $\partial G$; a Neumann condition remains as a boundary flux.

For homogeneous Dirichlet data, set

and define

Equivalently, $A:V\to V^*$ is defined through

Given $u_0\in H$ and $f\in L^2(J;V^*)$, the weak problem is to find

such that

with $u(0)=u_0$.

The embeddings

form a Gelfand triple. They are continuous and dense. A function in $V$ can be regarded as an element of $H$, and the Riesz representation theorem identifies $H$ with $H^*$. Each $h\in H$ therefore defines a continuous functional on $V$ by

The three spaces are not equal. The embeddings let the same object be read either as a regular function or as a weaker functional. This makes $a(u,v)$ meaningful for $u(t)\in V$, permits $\dot u(t)\in V^*$ to act on $v\in V$, and yields the energy identity

A sufficiently regular strong solution satisfies the weak equation. A weak solution need not possess all classical derivatives required by the strong equation.

1.3 Well-posedness

Writing a weak problem does not by itself guarantee that a solution exists. The standard parabolic theorem uses two estimates.

The first is continuity:

The estimate implies, for fixed $v$,

which says that $A:V\to V^*$ is bounded.

The second estimate is the Gårding inequality:

Diffusion supplies positive control of the gradient. Convection and reaction terms need not be positive, but they may be absorbed into the lower-order $H$ term. Under these assumptions, the weak problem has a unique solution and

This estimate gives continuous dependence on the initial condition and forcing term. In particular, $u_0$ only has to belong to $L^2(G)$. A discontinuous digital payoff can therefore generate a unique weak solution, although optimal error estimates may require additional regularity or positive smoothing time.

1.4 Galerkin semidiscretization and assembly

Choose

Every $v_h\in V_h$ is a linear combination of the basis functions. If they are linearly independent, the coefficient vector is unique. Substitute

into the weak equation and test successively with $\phi_i$. The result is

where

On a one-dimensional element of length $h$, continuous piecewise-linear basis functions give

and a constant diffusion coefficient gives

Adding local entries to their global degrees of freedom produces sparse matrices. A convection contribution generally makes $A$ nonsymmetric; this reflects the first-order drift and is not an assembly error.

1.5 The theta scheme and stability

Let $k=t_{m+1}-t_m$ and define

The $\theta$-scheme is

Stability comparison. Here $z=k\lambda\ge0$.

Accuracy and nonsmooth payoffs. The error orders assume a sufficiently smooth solution.

The stated total-error orders require enough spatial and temporal regularity. Stability may remain valid even when the formal convergence order does not.

For $0\le\theta<1/2$, stability requires

For linear finite elements, $\lambda_A\asymp Ch^{-2}$, so

For $\theta\ge1/2$, no such mesh-ratio condition is needed. This follows by considering a generalized eigenmode $Av=\lambda Mv$. With $z=k\lambda$ its one-step amplification factor is

Explicit Euler is stable only for $0\le z\le2$. Backward Euler satisfies $R_1(z)\to0$ and therefore damps high frequencies. Crank–Nicolson is stable for every $z\ge0$, but $R_{1/2}(z)\to-1$; unresolved modes may alternate in sign and decay very slowly.

"Unconditionally stable" does not mean accurate for arbitrarily large $k$. Explicit and implicit Euler have global temporal error $O(k)$ for a smooth solution. Crank–Nicolson has error $O(k^2)$. A kink or jump at $\tau=0$ may destroy the time derivatives needed for the second-order argument and cause order reduction.

1.6 Convergence and error estimates

Split the error at $t_m$ as

The approximation properties of $V_h$ control $\eta^m$. The term $\xi_h^m$ satisfies a discrete parabolic equation with a residual and is controlled by stability. For a sufficiently smooth solution and continuous piecewise-linear elements,

while

These orders are not properties of an algorithm name alone. Payoff singularities, degenerate diffusion, free boundaries, and jumps may violate the required regularity. That is the first point to check before reusing a Black–Scholes solver for another model.

1.7 A fully worked Black–Scholes example

The preceding derivation reduces FEM to

This example starts from a standard European put and computes the basis functions, element integrals, assembled matrices, and first time step explicitly. The hand calculation deliberately uses a very coarse mesh; the interactive cell then refines the problem and compares it with the Black–Scholes formula.

Step 1: transform the Black–Scholes equation

The strong equation for a European put is

Set

Since $V=Ku$ and $\partial x/\partial S=1/S$,

Dividing by $K$ gives

with

Truncate the real line to $G=(x_L,x_R)$. For a put, use

These are nonhomogeneous Dirichlet conditions. Their contribution to the right-hand side must be retained.

Step 2: basis functions and weak formulation

Let $x_0\lt x_1\lt\cdots\lt x_{N+1}$ be a uniform grid with spacing $h$. The piecewise-linear hat functions satisfy

On $K_e=[x_{e-1},x_e]$, introduce

Then

After multiplying by $v\in H_0^1(G)$ and integrating the second-order term by parts,

where

On a bounded interval this form is continuous on $H_0^1(G)\times H_0^1(G)$. Cauchy–Schwarz, Young, and Poincaré control the convection term and give a Gårding inequality. After subtracting a boundary lift, every $u_0\in L^2(G)$ therefore generates a unique weak solution

Step 3: compute the element matrices

The element mass matrix is

The required integrals are

Hence

The diffusion matrix is

For example,

The convection matrix is

For instance,

The complete element matrix is therefore

For $N$ interior nodes, assembly gives

while

and

$C_{II}$ is skew-symmetric. Unless $b=0$, the Black–Scholes stiffness matrix is nonsymmetric.

Step 4: insert numbers and take one time step

Use

Then $\alpha=0.02$, $b=0.03$, and the nodes are

The element matrices are

There are three interior degrees of freedom. The assembled blocks are

and

with the boundary-coupling blocks

Initially,

Take one backward-Euler step with $k=0.25$. At $\tau_1=0.25$,

Partitioning interior and boundary degrees of freedom gives

The right-hand side is

and the $3\times3$ solve yields

Since $x=0$ corresponds to $S=K=100$, the coarse-grid price is

The exact price is $3.3728$. This large error is useful: the diffusion scale over a quarter year is about $\sigma\sqrt{\tau}=0.1$, while the hand mesh has $h=0.5$. It cannot resolve the payoff kink. The small negative value at the rightmost interior node is another symptom of under-resolution. Correct assembly and adequate resolution are different questions.

The following cell repeats the assembly on $[-5,2]$ with $h=0.05$ and advances the solution for one year with Crank–Nicolson.

After refinement, the FEM price at $S=K$ is about $5.5534$, compared with the exact value $5.5735$. The absolute error is approximately $2.01\times10^{-2}$.

2. Solving the Black–Scholes equation in practice

Section 1 used Black–Scholes to make the finite-element algebra concrete. This section returns to the same model from the risk-neutral dynamics, connecting the SDE, generator, pricing PDE, localization, numerical treatment of nonsmooth payoffs, and Greeks. Once the resulting bilinear form is continuous and satisfies a Gårding inequality, the Galerkin, $\theta$-scheme, stability, and error machinery from Section 1 can be reused.

2.1 Risk-neutral dynamics and the generator

Assume that a non-dividend-paying stock follows

under the risk-neutral measure. Itô’s formula gives

The infinitesimal generator measures the instantaneous conditional change of a test function:

The drift term in Itô’s formula gives

The second-order term comes from quadratic variation. Correlated multidimensional noise produces mixed derivatives; jumps produce an integral term. Both appear in Section 3.

2.2 Feynman–Kac and the pricing PDE

For payoff $g(S_T)$,

Applying Itô’s formula to $e^{-rt}V(t,S_t)$ and setting its drift to zero gives

or

The conditional expectation and the PDE are probabilistic and analytic representations of the same price.

2.3 Time reversal, log-moneyness, and the weak form

The calendar-time PDE has terminal data at $t=T$. Introduce $\tau=T-t$ so that the payoff becomes initial data at $\tau=0$, and set

This maps $s>0$ to $x\in\mathbb R$, moves the strike to $x=0$, makes the diffusion coefficient constant, and normalizes the price by $K$. Since

the transformed equation is

where

The normalized put and call payoffs are

On $H^1(\mathbb R)$, integration by parts gives

Hölder’s inequality gives continuity, while

is a Gårding inequality. The well-posedness theorem from Section 1 therefore applies.

2.4 Localization

The natural log-price domain is $\mathbb R$, but a computer needs a bounded interval. Choose

With zero Dirichlet data, the localized problem has the probabilistic interpretation of a knock-out claim that dies when the process exits $G$. For a payoff of polynomial growth and a fixed interior point $x$,

The interior error decays exponentially as $R$ grows, although a zero boundary is not the correct far-field behavior of a vanilla claim. In implementations one often uses put or call asymptotics instead. Error should always be measured away from the artificial boundary.

Localization is part of the numerical approximation, not a harmless preprocessing step. At minimum,

Refining $h$ and $k$ while keeping $R$ too small eventually produces an error plateau.

2.5 Discretization and nonsmooth payoffs

On $G$, linear FEM gives

The time step remains

The same PDE can be discretized by FDM. On a uniform one-dimensional grid, both methods produce tridiagonal systems. The useful comparison is not a claim that one is universally better, but a record of the boundary treatment and error sources, checked against the same closed-form price.

Vanilla payoffs have a kink at the strike and a digital payoff is discontinuous. Well-posedness remains valid, but the optimal error theorem from Section 1 requires more regularity. Common remedies are:

• initialize with an $L^2$ projection rather than nodal interpolation;
• use a few backward-Euler half-steps before switching to Crank–Nicolson, known as Rannacher smoothing;
• grade the time mesh near $\tau=0$.

A graded mesh may be written as

It resolves the initial smoothing layer with short steps and increases the step size after diffusion has regularized the solution. The grading exponent must match the initial singularity; overly weak grading misses the layer, while excessive grading increases conditioning and cost. Rannacher smoothing addresses the same problem by damping high-frequency components before Crank–Nicolson is used.

2.6 Greeks

Delta and Gamma are derivatives with respect to the state:

In log-price coordinates,

They may be obtained by differentiating the finite element function or by postprocessing the final price vector. Differentiation amplifies mesh noise, and Gamma requires two more spatial derivatives than the price. Convergence of the price alone does not imply that Gamma is accurate.

Vega is a parameter derivative. If

then $\widetilde u=D_\eta u[\delta\eta]$ satisfies

The discrete sensitivity equation reuses the left-hand matrix from the price solve and adds a source term generated by the price. This is more systematic than bump-and-revalue and is better suited to error analysis.

The next cell uses Crank–Nicolson FDM to keep the browser code short. FEM uses the same time stepping and the same error diagnostics.

3. Seven extensions of the Black–Scholes framework

Black–Scholes is the template. Changing the model or contract usually changes the state variables, generator, energy space, boundary conditions, or one block of the discrete algebra; the overall route remains recognizable.

Each subsection follows the same checklist: specify a contract and parameters, derive the strong and weak forms, state the conditions for existence and uniqueness, identify the discrete system, and finish with an executable numerical example.

The common well-posedness tool is the following. Let

be a Gelfand triple. If $a(t;\cdot,\cdot)$ is measurable and there are constants $C,\alpha>0$ and $\lambda\ge0$, independent of $t$, such that

and

then $u_0\in H$ and $f\in L^2(0,T;V^*)$ determine a unique solution

For each model below, the main question is how to obtain these two estimates, or which replacement theorem is needed when the state space is weighted, constrained, or fractional. A numerical plot checks an implementation for one parameter set; it does not prove well-posedness or convergence.

3.1 CEV: degenerate diffusion at the boundary

The CEV model replaces $\sigma S_t$ by

The second-order coefficient is $\tfrac12\sigma^2s^{2\rho}$ and vanishes at $s=0$. An ordinary $H^1$ norm no longer matches the energy of the operator. A natural norm is

For some values of $\rho$, the measure also needs a weight $s^{2\mu}$, leading to a space $W_{\rho,\mu}$.

For a European put with $K=1$, the strong equation in time to maturity is

On $G=(0,R)$, use $u(\tau,0)=e^{-r\tau}$ and $u(\tau,R)=0$. Set

Since

the unweighted weak form is

where

With measure $s^{2\mu}ds$, integration by parts adds the derivative of the weight. The first-order coefficient becomes

This is the weighted weak form; it is not obtained by merely multiplying the unweighted final formula by a weight.

Existence and uniqueness. On $[\varepsilon,R]$ with $\varepsilon>0$,

so ordinary $H^1$ continuity and a Gårding inequality apply. If the domain includes $s=0$, uniform ellipticity is lost. In $W_{\rho,\mu}$, suitable weight parameters make the mass, drift, and trace terms finite and recover a weighted Gårding inequality. The same parabolic theorem then gives a unique weak solution in the weighted Gelfand triple. Whether a boundary condition must be imposed at zero also depends on the boundary classification of the CEV process.

For piecewise-linear basis functions,

and

The coefficients vary with $s$, so the element matrices are evaluated by quadrature rather than copied from one constant matrix.

The first cell shows the degeneracy directly.

Now take

Crank–Nicolson FEM on $60$, $120$, and $240$ elements gives $0.093404$, $0.093589$, and $0.093635$ at $S_0=1$. The successive differences decrease by roughly a factor of four. This is a mesh-consistency check rather than a comparison with a closed-form CEV price.

3.2 American options: early exercise as an obstacle

CEV changes the coefficients of the pricing operator. American exercise changes the problem class itself: the price must satisfy both the evolution equation and an inequality constraint.

An American claim is an optimal stopping problem:

The value must satisfy $V\ge g$. In the continuation region $V>g$, the Black–Scholes equation holds; in the stopping region $V=g$. With

the strong complementarity system is

In forward time, define

Then

If $u>g$, the PDE residual is zero. If $u=g$, the residual is nonnegative.

After the log-price transform and a boundary lift, define the closed convex set

The weak formulation is to find $u(\tau)\in K$ such that

Existence and uniqueness. $K$ is nonempty, closed, and convex. The localized Black–Scholes form is continuous and satisfies a Gårding inequality. An exponential shift makes it coercive, after which the standard parabolic variational-inequality theorem gives existence. Testing two solutions against one another and using monotonicity and Grönwall’s inequality gives uniqueness.

Implicit time stepping and FEM or FDM in space produce an LCP:

Each degree of freedom either satisfies the pricing equation in the continuation region or is fixed to the exercise payoff. For backward Euler, $B=M+kA$. The mass matrix and Gårding lower bound make its symmetric part positive for sufficiently small $k$; a standard monotone discretization gives a $P$-matrix and a unique finite-dimensional LCP solution.

Projected SOR first performs an SOR update,

and then projects onto the obstacle:

PSOR is an iterative solver for the LCP; it is not the early-exercise condition itself.

For $K=100$, $r=5\%$, $\sigma=25\%$, and $T=1$, using $160$ spatial intervals on $[0,300]$ and $160$ backward-Euler steps gives an exercise boundary near $S^*=75.0$. Interpolating to $S=100$ gives an American value of about $7.9667$, a European value of $7.4648$, and an early-exercise premium of $0.5018$. The most expensive time step uses $21$ PSOR iterations.

A credible LCP calculation should also report the equation residual, obstacle violation, and complementarity product. The iteration tolerance must be smaller than the target discretization error.

3.3 Path dependence: restoring the Markov state

The American obstacle adds a constraint without adding a state variable. Path dependence usually does the reverse: it restores Markovian structure by enlarging the state.

Path dependence does not automatically rule out PDE methods. The question is whether the relevant history can be summarized by a finite-dimensional state. If $Y_t$ satisfies an SDE or ODE driven by the current $(S_t,Y_t)$, the augmented process is Markov and Feynman–Kac still gives a finite-dimensional PDE. The cost is an extra dimension or a jump condition at monitoring dates.

Barrier options. A continuously monitored barrier is summarized by the first hitting time

For a knock-out claim, solve the Black–Scholes equation in the surviving region and impose an absorbing boundary at the barrier. With discrete monitoring, solve the PDE between monitoring dates and apply

Without a rebate, in/out parity provides a useful check:

Asian options. Define the accumulated price

Since $dA_t=S_tdt$, $(S_t,A_t)$ is Markov. Under Black–Scholes its generator is

For a continuous arithmetic-average call, the calendar-time equation is

with $V(T,s,A)=(A/T-K)^+$. In time to maturity the strong equation is

On a truncated rectangle $G$, let $a(s)=\sigma^2s^2/2$. Integrate only the $s$-diffusion by parts. The weak form is

where

Existence and uniqueness. There is no $w_Av_A$ term, so this form does not control the full two-dimensional $H^1$ norm. The natural setting is a Kolmogorov transport-diffusion space: a weighted first derivative in $s$, the transport derivative $s u_A$ in the dual space, and data on the inflow boundary. Away from $s=0$, with bounded coefficients and compatible inflow data, the closed transport-diffusion operator generates a unique evolution. A common alternative is to add $-\varepsilon_Au_{AA}$, solve a uniformly parabolic problem for each $\varepsilon_A>0$, and pass to the limit $\varepsilon_A\downarrow0$.

With tensor-product basis functions $\Phi_{ij}(s,A)=\phi_i(s)\psi_j(A)$,

where

The next cell supplies an independent benchmark for a future two-dimensional FEM implementation. With $S_0=K=100$, $r=5\%$, $\sigma=20\%$, $T=1$, and $128$ monitoring dates, $60{,}000$ fixed-seed paths give about $5.8209$ with Monte Carlo standard error $0.0329$.

Compound options. If $C(t,s)$ is the price of the underlying option and a compound call pays $K_c$ at $t_c$ to acquire it, first compute $C(t_c,s)$ and then use

as the terminal condition for a second PDE.

Swing options. Let $q$ be the number of exercise rights remaining. At an exercise date,

Between exercise dates, solve one linear pricing PDE for each $q$; at an exercise date, apply the maximum operation. Continuous exercise produces a coupled system of variational inequalities.

3.4 Interest-rate derivatives: discounting as a state variable

State augmentation also clarifies interest-rate models: the short rate is both a state variable and the quantity that drives discounting.

Let the short rate be a one-dimensional Markov diffusion under the risk-neutral measure:

The money-market account is

For a terminal payoff $g(r_T)$,

and Feynman–Kac gives

With $\tau=T-t$,

On $G=(r_{\min},r_{\max})$, let $a(r)=\sigma(r)^2/2$. Since

the diffusion, convection, and reaction coefficients in divergence form are $a$, $a'-b$, and $r$. Boundary data should follow from asymptotics or boundary accessibility, not from the truncation alone.

For a Vasicek model,

the strong zero-coupon equation, written with $y$ as the short-rate coordinate, is

The weak form is

where

Existence and uniqueness. On a bounded interval the drift and reaction coefficients are bounded, and $\sigma>0$ gives uniform diffusion. The form is continuous. Young’s inequality controls the drift, while the negative part of $y$ is absorbed into the Gårding shift. A unique weak solution follows even when $y_{\min}<0$.

For piecewise-linear elements,

The Vasicek zero-coupon price is also known in closed form:

where

This gives an external error benchmark. With $\kappa=0.6$, $\bar r=0.04$, $\sigma=0.02$, $T=5$, $160$ elements on $[-0.05,0.15]$, and $500$ Crank–Nicolson steps, the FEM value at $r_0=3\%$ is $0.833031769$ versus the exact $0.833032037$. The maximum nodal error is about $2.93\times10^{-7}$.

The CIR model replaces the constant diffusion by $\sigma\sqrt r$. It preserves nonnegativity but degenerates at zero, requiring the same weighted-space and boundary classification ideas as CEV.

3.5 Multiple assets: the multidimensional template

Once the state is vector-valued, the same generator-to-weak-form derivation produces the multidimensional template used for several assets or factors.

Let $W_t\in\mathbb R^m$ be Brownian motion and let the risk-neutral state $X_t\in\mathbb R^d$ satisfy

Multidimensional Itô calculus gives the generator

With killing rate $c(x)$ and terminal payoff $g(X_T)$,

satisfies

Set $\tau=T-t$, $D=Q/2$, and

The forward equation is

On a bounded domain $G\subset\mathbb R^d$, with homogeneous Dirichlet data for simplicity, choose $V=H_0^1(G)$. The weak form is

where

If there is a position-independent $\lambda>0$ such that

then diffusion controls the full gradient. Cauchy–Schwarz, Young, and Poincaré control the remaining terms, giving continuity and a Gårding inequality.

In a multi-asset Black–Scholes model with volatilities $\sigma_i$ and correlations $\rho_{ij}$, log-price coordinates give

Off-diagonal entries of $Q$ generate mixed derivatives. The route from SDE to weak form remains unchanged:

On a product domain, use

In two dimensions,

For constant or separable coefficients, the mass matrix is $M_1\otimes M_2$; pure second derivatives give $S_1\otimes M_2$ and $M_1\otimes S_2$; mixed derivatives give Kronecker products of first-derivative matrices. General domains and nonseparable coefficients still admit ordinary multidimensional assembly, but not a short Kronecker representation.

For a two-asset arithmetic basket call, set $x_i=\log(S_i/K)$ and $b_i=r-\sigma_i^2/2$. The strong equation is

with

The constant diffusion matrix is

and

Existence and uniqueness. For $\sigma_1,\sigma_2>0$ and $|\rho|<1$,

$D$ is positive definite, its smallest eigenvalue controls $\|\nabla v\|_{L^2}$, and the weak problem has a unique solution. At $|\rho|=1$, this argument fails because the diffusion is lower rank.

If both directions use matrices $M,S,C$, with $C_{ij}=\int\phi_j'\phi_i$, then

This formula shows exactly where the mixed derivative enters the global matrix.

For $S_{1,0}=S_{2,0}=K=100$, $r=5\%$, $\sigma_1=20\%$, $\sigma_2=30\%$, $T=1$, and $\rho=0.4$, $100{,}000$ fixed-seed paths give a basket-call benchmark of about $10.8480$ with standard error $0.0497$.

If the error remains $O(h^2)$ while the number of degrees of freedom is $N\asymp h^{-d}$, then

This deterioration with dimension is the curse of dimensionality. Sparse grids, low-rank tensor methods, and Monte Carlo can all be understood as responses to this scaling.

3.6 Heston: a degenerate multidimensional extension

Heston is a two-dimensional diffusion of the preceding type, with variance as the second state and a new degeneracy at the boundary $v=0$.

The Heston model introduces stochastic variance:

In $(x,v)=(\log S,v)$, the covariance matrix is

The model contains a mixed derivative, but the entire covariance matrix is multiplied by $v$. At $v=0$ its smallest eigenvalue vanishes and uniform ellipticity is lost.

For a European call, let $u(\tau,x,v)=V(T-\tau,e^x,v)$. The strong equation is

with

Define

Since

and

the weak form is

where

where

Existence and uniqueness. For $|\rho|<1$ and $v>0$, $D(v)$ is positive definite, but its smallest eigenvalue is proportional to $v$. The natural energy controls

with a weight $\omega$ compatible with the CIR variance process. Suitable weights and boundary conditions make $a_H$ continuous and give a weighted Gårding inequality, hence a unique weak solution in a weighted Gelfand triple. The Feller condition

makes $v=0$ inaccessible from a strictly positive starting value, but it does not restore uniform ellipticity because $D(0)=0$.

For finite element basis functions $\Phi_j$,

The matrix has the same structural blocks as the two-asset problem, but the coefficients vary with $v$ and elements near $v=0$ must be interpreted in the weighted energy.

The first cell displays the tensor-product cost and the collapse of the smallest diffusion eigenvalue.

For an actual price, take

Full-truncation Euler keeps the coefficients nonnegative in the variance update. At $\rho=-0.5$, $60{,}000$ paths and $250$ steps give a call benchmark of about $9.1027$ with standard error $0.0489$. This is an external validation value for a two-dimensional FEM calculation; its own time-discretization bias must be checked separately.

Delta and Gamma are still state derivatives. Sensitivity to $\rho$ or another model parameter can be obtained from a parameter-derivative equation. Under the required regularity, prices and sensitivities can achieve the same convergence order in Heston and two-dimensional Black–Scholes discretizations.

3.7 Lévy models: a nonlocal generator

Heston remains a local diffusion: its differential operator only couples nearby points after FEM assembly. Lévy jumps introduce a nonlocal term and therefore change both the analysis and the matrix structure.

Continuous diffusions cannot independently control skew, heavy tails, and jumps. For a Lévy log-return process, the generator contains

where $\nu$ is the Lévy measure. The pricing equation becomes a PIDE.

Local differential operators couple nearby basis functions. A jump couples $x$ to every reachable $x+z$, so the stiffness matrix is generally dense. The energy space may also cease to be $H^1$: the high-frequency growth of the symbol determines an appropriate fractional Sobolev space $\widetilde H^s(G)$. If the resulting form is continuous and satisfies a Gårding inequality, the well-posedness and Galerkin framework from Section 1 remains valid.

As a concrete finite-activity example, consider Merton jump diffusion. Let $N_t$ have intensity $\lambda$, let one log jump satisfy

and define

The risk-neutral continuous log drift is

For a European call, the strong PIDE is

with $u(0,x)=(e^x-K)^+$. Because $\lambda<\infty$, the small-jump compensation term is unnecessary here; its risk-neutral correction is already contained in $b$.

The weak form is

where

Existence and uniqueness. If $\sigma>0$, local diffusion controls the $H^1$ seminorm. For finite jump activity, translation is norm preserving on $L^2(\mathbb R)$, so the jump term is a bounded $L^2$ perturbation. Continuity and a Gårding inequality follow, and the localized weak solution is unique. If $\sigma=0$ and the Lévy measure has infinite activity, this $H^1$ argument is no longer appropriate; the symbol determines a fractional energy space.

The nonlocal finite element block is

Even basis functions with disjoint local supports can be connected by a jump, which explains the dense matrix.

Localization now produces at least four error sources:

plus any algebraic approximation used to compress the dense matrix. The next cell shows how a Gaussian jump kernel spreads a localized function.

Merton’s model also has an independent Poisson-mixture price. Conditional on $N_T=n$, the terminal log-price is normal with an additional mean $n\mu_J$ and variance $n\delta_J^2$. Weight the corresponding lognormal call prices by

For $S_0=K=100$, $r=3\%$, $\sigma=18\%$, $T=1$, $\lambda=0.6$, $\mu_J=-0.1$, and $\delta_J=0.25$, the mixture price is $11.5354$. Setting $\lambda=0$ gives the Black–Scholes benchmark $8.6407$.

Infinite-activity models such as Variance Gamma require more careful quadrature and fractional-space analysis, but the computational chain is unchanged: identify the order of the generator, choose the energy space, prove continuity and a Gårding inequality, localize, and assemble the nonlocal matrix. The same sensitivity framework covers diffusion, Heston, and Lévy parameters.

3.8 One computational pipeline

All seven extensions fit the same diagram:

flowchart LR
    A[State process and payoff] --> B[Generator]
    B --> C[PDE / PIDE / obstacle problem]
    C --> D[Weak form and energy space]
    D --> E[Localization and boundaries]
    E --> F[Spatial discretization]
    F --> G[Theta time stepping]
    G --> H[Price, Greeks, and error]

• CEV changes the diffusion weight and energy space.
• American exercise replaces an equation by an obstacle constraint.
• Path dependence is handled by boundaries, state augmentation, or dynamic recursion.
• Short-rate models make the discount term state dependent.
• Multi-asset models replace scalar diffusion by a covariance matrix and tensor-product spaces.
• Heston adds boundary degeneracy to the multidimensional framework.
• Lévy models replace a local second-order operator by a nonlocal integral operator.

The reusable object is not one particular matrix. It is the sequence from stochastic dynamics to generator, weak formulation, well-posedness, discretization, and error control.