Black-Scholes Theory and Diffusion Processes on the Cotangent Bundle of the Affine Group

The Black-Scholes partial differential equation (PDE) from mathematical finance has been analysed extensively and it is well known that the equation can be reduced to a heat equation on Euclidean space by a logarithmic transformation of variables. However, an alternative interpretation is proposed in this paper by reframing the PDE as evolving on a Lie group. This equation can be transformed into a diffusion process and solved using mean and covariance propagation techniques developed previously in the context of solving Fokker–Planck equations on Lie groups. An extension of the Black-Scholes theory with coupled asset dynamics produces a diffusion equation on the affine group, which is not a unimodular group. In this paper, we show that the cotangent bundle of a Lie group endowed with a semidirect product group operation, constructed in this paper for the case of groups with trivial centers, is always unimodular and considering PDEs as diffusion processes on the unimodular cotangent bundle group allows a direct application of previously developed mean and covariance propagation techniques, thereby offering an alternative means of solution of the PDEs. Ultimately these results, provided here in the context of PDEs in mathematical finance may be applied to PDEs arising in a variety of different fields and inform new methods of solution.


Introduction
The Nobel Prize-winning Black-Scholes Equation [1] is arguably the most well-known partial differential equation in mathematical finance. The equation rests on a parsimonious option pricing model that has been fairly successful in informing banks and portfolio managers of the construction of risk-free hedges. Additionally, the model has since provided the framework for the import of a variety of tools, such as the theory of stochastic processes, from physics and mathematics. In this paper, we offer a new Lie group-theoretic interpretation of the Black-Scholes equation by reformulating the original equation and extensions of it as diffusion processes on Lie groups.
Group-theoretic approaches have been used extensively in the analysis of symmetries of partial differential equations (PDEs) in mathematical finance [2][3][4]. One of the central questions there is to identify the group of transformations of variables that can be applied to the equations that preserve the structure of the equations while reducing it to a simpler form for analysis and solution. For instance, the one-asset [2] and in general, the multi-asset Black Scholes equation [5] can be reduced to a heat equation through a logarithmic transformation of variables.
In this paper, instead of analysing the symmetry properties of PDEs, we reformulate the PDE as a diffusion equation on the Lie group. Functions of real parameters are upgraded to functions that take group elements as their arguments. That is, we start with a linear (parabolic) PDE of the form, where L is the linear partial differential operator and q is a vector of coordinates that may be chosen based on how the differential equation is derived. In the case of the Black-Scholes equation and related asset models considered in this paper, L and q are originally defined in an N-dimensional Euclidean space. By matching the derivatives in L to Lie derivatives of correctly chosen groups, one can rewrite (1) as the following linear PDE over the Lie group G, where q parameterizes the group element g ∈ G andf (g(q), t) = f (q, t). The operatorL is now a differential operator consisting of Lie directional derivatives. Such differential equations arise in bevel-tip needle steering [6,7], error propagation [8,9], DNA statistical mechanics [10], multi-robot localization [11], stochastic kinematic cart models in SLAM [12] and in image contour completion and enhancement [13][14][15]. The benefit of reframing differential equations this way would be that variable-coefficient PDEs of N variables in R N may reduce to constant-coefficient PDEs on the Lie group G of dimension of at least N, which can be analysed using techniques developed in [6][7][8][9][10][11][12]. These applications motivated the development of mean and covariance propagation techniques to approximate the solution of diffusion equations in the special Euclidean group SE(N) and more generally, in unimodular groups [16,17]. Other numerical techniques to solve differential equations on Lie groups, through a generalisation of Euler and Runge-Kutta schemes have been developed in [18][19][20].
In this paper, we extend the regime of applicability of mean and covariance propagation techniques by first considering the groups GL + (1) and GL + (1) × GL + (1) that arise in the one-asset and two-asset Black-Scholes equations, respectively, and then the affine group, A f f + (1) that arises in a coupled-asset dynamics extension of the Black-Scholes theory. The application of mean and covariance propagation to A f f + (1) becomes the main subject of the paper. There exist other methods to approximately solve PDEs in mathematical finance, such as finite difference methods [21,22], finite element methods [23] and the Adomian decomposition method [24,25]. In this paper, we also make a comparison between the mean and covariance propagation technique and a standard finite difference scheme used to solve the governing equations.
The Lie group A f f + (1) is not unimodular, thereby precluding the direct application of the theory of mean and covariance propagation used in the solution of diffusion equations on the group. However, we show that the cotangent bundle group of a Lie group (with a trivial center) is unimodular. Thus, we analyse diffusion processes on the affine group by matching the diffusion on the affine group with a degenerate diffusion on the cotangent bundle to which mean and covariance propagation techniques can be applied. This paper considers three types of asset dynamics models: one-asset, two-asset and a coupled asset model. Each of these models give rise to a PDE describing the evolution of the option price as a function of the asset prices, and are introduced in Section 2. These PDEs in Euclidean space are reframed as diffusion processes on Lie groups in Section 3. The theory of mean and covariance propagation in the solution of diffusion processes on Lie groups is reviewed in Section 4. The one-asset and two-asset Black-Scholes models, which can be trivially solved using the logarithmic transformation of variables, are used as examples to illustrate the techniques, and thereby sets the ground for the non-trivial case of a coupled asset model. Since the coupled asset model leads to a diffusion equation over a non-unimodular group, Section 5 proves that a related structure, the cotangent bundle group of a Lie group with trivial center, is always unimodular. Section 6 describes the mean and covariance propagation over the unimodular cotangent bundle of the affine group, thereby solving the option price dynamics for the coupled asset model. Section 7 provides some numerical results for the mean and covariance propagated solution in comparison with finite difference methods. Finally, Section 8 seeks to demonstrate backward compatibility and completes the analysis by applying mean and covariance propagation on the cotangent bundle group to the solution of the one-asset Black-Scholes equation.

Asset Dynamics Models
In this section, we review the derivation of the well-known one-asset Black-Scholes equation, the two-asset Black-Scholes equation that evolve with correlated Wiener process, and follow a similar derivation to develop a coupled asset model.

One-Asset Black-Scholes Equation
The one-asset Black-Scholes equation is derived from the following Itô stochastic differential equation [26] governing the dynamics of an asset value a, da/a = µdt + σdW, where dW is the increment of a Wiener process, corresponding to random draws from a Gaussian probability distribution [27,28]. The increment of the Wiener process satisfies the following relations: and in general for an N-dimensional Wiener increment dW we have, where in both cases · · · denotes the expected value with respect to the underlying probability distribution function and δ ij is the Kronecker-delta (δ ij = 1 if and only if i = j, otherwise δ ij = 0). The µ and σ in Equation (3) respectively describe a drift and volatility for the asset price evolution. The volatility is defined such that σ 2 a 2 is the variance of the random price fluctuation. Equation (3) is a stochastic differential equation whose solution is a geometric Brownian motion. We emphasise that the equation in (3) is to be interpreted as an Itô stochastic differential equation. Stratonovich stochastic differential equations also appear in this paper (see Equation (110) for instance) and will be distinguished from Itô equations with a , i.e., for a stochastically varying quantity x, a general one-dimensional Stratonovich stochastic differential equation with Wiener noise would be, for a drift a(x, t) and diffusion coefficient b(x, t). Nevertheless the properties of dW in (4) and of dW in (5) hold for both Itô and Stratonovich equations. The option price V = V(a, t) is a function of time t and the underlying asset price a. Using Itô's Lemma [29], we have for dV correct to order dt, In the Black-Scholes model, there also exists another asset, the bond price, which evolves deterministically (for a non-stochastic interest rate) as dB = rBdt where r is the risk-free interest rate [26]. A portfolio Π is constructed as a linear combination of the two assets and can be written as, This form of Π is due to the fact that the portfolio is assumed to be self-financing with no external money flows [30]. Banks choose ∆ in order to obtain a zero-risk portfolio [31], i.e., so that d(V + Π) would have no Wiener noise terms. Assets without stochasticity would evolve at the risk-free interest rate r and therefore we obtain, and using dΠ = ∆dV + βdB, a hedge of the form ∆ = −a∂V/∂a would cancel out the (σa∂V/∂a)dW term in (6) and therefore make the dynamics of V + Π risk-free. Substituting this hedge into (8) and using the form of dV in (6) we obtain the one-asset Black-Scholes equation as,

∂V ∂t
The parameters β and µ do not feature in this equation.

Two-Asset Black-Scholes Equation
In the two-asset model, one has asset values a and b that evolve with correlated Wiener processes-this is a specific case of the multi-asset model in [5]. That is, where dW 1 dW 2 = ρdt and · · · represents the expectation. In the two-asset problem, the option price is V = V(a, b, t). Applying Itô's Lemma to dV now gives: The portfolio is now Π = ∆ 1 a + ∆ 2 b + βB where ∆ 1 = −a∂V/∂a and ∆ 2 = −b∂V/∂b; ∆ 1 and ∆ 2 are both chosen such that V + Π experiences risk-free dynamics-very similar to the procedure in the derivation of the one-asset Black-Scholes equation in the previous section. The two-asset Black-Scholes equation [32] will then be, where r is the risk-free interest rate. Like the one-asset Black-Scholes equation, µ 1 , µ 2 and β do not feature in the two-asset equation.

Option Price Evolution with Coupled Assets
The evolution equations of the two assets a and b in the two-asset Black-Scholes model (10) featured correlated Wiener processes. In this section, we imagine a different form of coupling of the form: The evolution of the value of asset a is independent of the evolution of b; a is therefore the 'leading' asset value. On the other hand, the drift and variance in the stochastic differential equation for db are dependent on the current values of the leading asset. Hence b is the 'trailing' asset value. One may imagine this form of dependency when a represents a raw material and b represents a finished good that makes use of this raw material.
Both da and db are forced by the same Wiener process dW; this implies that σ 2 /σ 1 = (µ 2 − r)/(µ 1 − r) for the risk-free interest rate r. This can be derived by using the risk-neutral measure with the knowledge that all assets evolve with a risk-free interest rate of r in this measure [31]. The option price is V = V(a, b, t) and applying Itô's Lemma to dV now gives, We now assume a general hedge of the form Π = ∆ 1 a + ∆ 2 b + βB. Since there is only one Wiener process dW, removing uncertainty would not be a sufficient condition to determine both ∆ 1 and ∆ 2 . Thus, for simplicity, we let ∆ 2 = 0 which leads to the relationship ∆ 1 = −∂V/∂a − (σ 2 /σ 1 )∂V/∂b. The governing partial differential equation for V(a, b, t) with two coupled assets would be, where the leading asset/trailing asset coupling between the variables breaks the symmetry between the assets a and b that existed in (12).

Reframing Partial Differential Equations as Diffusion Processes on Lie Groups
The governing PDEs for the three types of asset dynamics models in (9), (12) and (15) can be re-expressed in terms of Lie derivatives of GL(1) + , GL(1) + × GL(1) + and A f f + (1), respectively, by matching the group parameters to the asset variables. A review of the concept of Lie derivatives is presented in Section 3.1.

Preliminary Definitions
Let G be an N-dimensional matrix Lie group with Lie algebra G. Then, let an element g ∈ G be parameterized as g = g(q) where q ∈ R N , using the notation in [29,33]. The 'right' Jacobian J R (g(q)) of the group is defined [34] as the following matrix, and the 'left' Jacobian J L (g(q)) is the matrix, where the square brackets reinforce that we are dealing with a matrix. This is not to be confused with the 'right' and 'left' Jacobian determinants that arise in the volume forms of the group, which would be the determinants of the matrices in (16) and (17), respectively. Note that the 'right' Jacobian has the g −1 appearing on the right whereas the 'left' Jacobian has the term on the left. However, the 'right' Jacobian is left invariant, i.e., J R (h • g(q)) = J R (g(q)) and the 'left' Jacobian is right invariant, i.e., J L (g(q) • h) = J L (g(q)), assuming that q parameterizes the whole group G and that these shifts are permitted in the function domain. Finally, the ∨ operator is defined as a bijection mapping G to R N , and vectorizes the matrix element in G. The inverse of the ∨ operator is a ∧ that maps R N to G.
We make an additional remark regarding the parameterization of the group with q. This is to say that the whole group (except for a set of measure zero) is parameterized by one coordinate chart. For instance, for groups such as SO(3) or SE(3) where one may use Euler angles to parameterize the rotations, there exists a set of measure zero corresponding to the set of Euler angles where the Jacobian matrices for the parameterization becomes singular. For (α, β, γ) denoting the ZXZ or ZYZ Euler angles, singularities occur at β = 0 and β = π. Additionally, since a rotation by (α, β, γ) describes the same rotation as that by (α + 2π, β, γ + 2π), the open coordinate chart will be (α, β, γ) ∈ (−π, π) × (0, π) × (−π, π), which has a one-to-one correspondence with a subset of rotations that excludes the rotations at β = 0 and π, and the rotations at α = π and γ = π. The closure of this coordinate chart will however establish a many-to-one map with the group and parameterize all group elements. A similar issue exists in the case of using the Iwasawa decomposition to parameterize SL(2, R) with (θ, λ, ξ) where θ parameterizes the 2D rotation, and (λ, ξ) parameterize the 2 × 2 upper triangular matrix. Here, the coordinate chart would be (θ, λ, ξ) ∈ (−π, π) × R + × R; yet again, the closure of this chart will establish a many-to-one map with SL(2, R) and parameterize all group elements. For other cases, such as using a vector drawn from R N 2 to parameterize GL(N) as well as for the affine group A f f + (1) and the cotangent bundle group of A f f + (1) considered in this paper, the coordinate chart used parameterizes the entire group.
The operator of the adjoint representation of group G at g ∈ G is given by Ad(g) where Ad(g)X = gXg −1 for any X ∈ G. When expressed as a matrix, we denote the operator as [Ad(g)]. If E i for i = 1, · · · , N forms an N-dimensional orthonormal basis for G, [Ad(g)] is given as, Orthonormality of the Lie algebra basis is defined here with respect to an inner product of the to define the Lie bracket, we can also represent the "little ad" operator, where ad( where X ∈ G and Ad(exp(X)) = exp(ad(X)). For a differentiable function on the group f : G → R, one can construct the right and left Lie directional derivatives as, In parametric form, wheref (q) = f (g(q)),J R,L (q) = J R,L (g(q)) and ∇ q = [∂/∂q 1 , · · · , ∂/∂q N ] T is the gradient operator for R N . Here, '−T' represents the inverse transpose operation, i.e., the inverse of the transpose of the matrix. The Lie directional derivative operators are E R = [E R 1 , · · · , E R N ] T and E L = [E L 1 , · · · , E L N ] T . In the sequel, we drop the tildes, as it will be clear from the arguments whether f (g(q)) orf (q) is considered, and likewise for the Jacobians.

One-Asset Black-Scholes as a Diffusion on Gl + (1)
The set of positive real numbers equipped with the multiplication operation forms a commutative (Abelian) group. This is a subgroup of the general linear group of one dimension and is represented by GL + (1). The basis of the Lie algebra is the number 1 and for a group element a ∈ GL + (1), the corresponding element in the Lie algebra is log a. Using (21) for a differentiable function f : where both left and right Lie derivatives yield the same result since the group is Abelian. Using this relationship, we can rewrite the one-asset Black-Scholes Equation (9) as,

∂V ∂t
using the shorthand E 2 V = E(EV). The variable-coefficient Black-Scholes equation has thus been transformed to a constant-coefficient equation on GL + (1). Additionally, using time reversal t = −t to convert the backward parabolic equation to a forward parabolic equation and setting V(a, t ) = u(a, t )e −rt we have, which is a diffusion equation with drift r − σ 2 /2 and diffusivity σ 2 . Here, u = u(a, t ) is interpreted as a function over GL + (1) × R ≥0 . Henceforth the term initial condition will be with respect to t , which due to time reversal, refers to a final condition with respect to t. We note that (9) by itself is not a diffusion process. Instead, we obtain a diffusion equation in (24) only after the time reversal t = −t and making the exponential transformation V = u(a, t )e −rt . Therefore, solving for u(a, t ) indirectly solves for V(a, t ) and solves the Black-Scholes equations. In subsequent sections, we apply similar transformations, while noting that solutions of such diffusion equations for u indirectly provides a solution for V.
3.3. Two-Asset Black-Scholes as a Diffusion on Gl + (1) × Gl + (1) An element g ∈ GL + (1) × GL + (1) parameterized as g = g(a, b) can be represented by a diagonal matrix as, The Lie algebra of GL + (1) × GL + (1) can be represented by the orthonormal basis, The group is Abelian and much like GL + (1), the left and right Lie derivatives coincide for a differentiable function f : GL + (1) × GL + (1) → R as, where f (g(a, b)) =f (a, b). This allows us to write the two-asset Black-Scholes Equation (12) in terms of Lie derivatives as, where t = −t and V(a, b, t ) = u(g(a, b), t )e −rt . Hence, the solution of (28) allows one to construct the solution to the two-asset Black-Scholes Equation (12).

Option Price Evolution with Coupled Assets as a Diffusion on A f f + (1)
The affine group of the positive real line consists of all (a, b) ∈ R + × R that transforms the scalar The group action is a matrix multiplication of g with x = [x, 1] T (where x is expressed in homogeneous coordinates as x). The Lie algebra of A f f + (1) is two dimensional and spanned by the orthonormal basis elements E 1 and E 2 , An important feature of A f f + (1) is that the determinant of [Ad(g)] is a, which is generally not equal to 1. This implies that the group is not unimodular. For non-unimodular groups, there exist distinct left-invariant and right-invariant Haar measures [35,36]. Using (16) and (17), the left and right Jacobians J L and J R of A f f + (1) expressed as matrices are, and since det J L = det J R the Lie group is not unimodular. The left and right Lie directional derivatives can be evaluated using (21) for a function f : where f (g(a, b)) =f (a, b). Therefore, we can now rewrite (15) in terms of Lie derivatives of A f f + (1) as where t = −t and V(a, b, t ) = u(g(a, b), t )e −rt . Hence, the solution of (34) allows one to construct the solution to the coupled asset model in Equation (15).

Mean and Covariance Propagation for Unimodular Lie Groups
Diffusion equations on unimodular Lie groups, such as (24) and (28) can be solved approximately using mean and covariance propagation techniques developed in [6][7][8][9][10][11][12]16]. This technique would not be applicable directly on the affine group in the solution of Equation (34) since the group is not unimodular; instead, the later sections will show how the technique can be modified by converting (34) to a diffusion over the cotangent bundle of the affine group, which is unimodular. This section will review the theory of covariance propagation for a general N-dimensional unimodular matrix Lie group G with Lie algebra G and apply the technique to solve (24) and (28).
Consider a general diffusion equation for u(g, t) on G in terms of right Lie derivatives as, where the drift vector m(t) = [m 1 (t), · · · , m N (t)] T and diffusivity matrix D(t) are independent of g but can be time-dependent in general. Since the equation is a Fokker-Planck equation over a Lie group [17,33], one can interpret u(g, t) as a probability density (see Appendix B). This analogy also assumes that u(g, t) is square-integrable in G. Additionally, G u(g, 0)dg = 1, which also holds true for all values of time. Here, dg is the Haar measure [35,36].
Given an initial condition u 0 = u(g, 0), we aim to propagate the solution to the next time-step (t = δt) to obtain u(g, δt). In general, knowing u(g, t) we seek to obtain the solution at u(g, t + δt). The solution at t + δt can be obtained by a group convolution [33] from the solution at t as, where u f (g, δt) is the fundamental solution/Green's function that propagates the solution over δt.
We recognise that over a sufficiently small δt the fundamental solution evolves in G. Since G can be bijectively mapped to the Euclidean space R N , an approximate expression for u f (g, δt) would be a multivariate Gaussian that is a function of the N-dimensional q used to parameterize the group element g ∈ G. This Gaussian would evolve in R N with a drift of mδt. This drift can be injected on the group G by an evolution of the form µ(δt) = exp(m ∧ δt) where ∧ maps an element of R N to G [12,17]. The covariance of this Gaussian would be Dδt, which would equal the group-theoretic covariance Σ(δt) when restricted to small δt. The definitions of the group-theoretic mean, µ(t), and covariance, Σ(t), would be provided later in (39) and (40), respectively, but here it suffices to note that the group-theoretic mean and covariance match the mean and covariance defined in Euclidean space R N (see Appendix C) for small δt, thereby allowing us to construct u f (g, δt) in terms of the group parameters [17,33] as, where the log maps the group element to G and ∨ maps the element of the Lie algebra to R N . This form naturally extends to a Gaussian over the group in (47). If u 0 = u(g, 0) = δ G (gg −1 0 ), which is a group Dirac delta distribution centred about g 0 , the distribution at δt remains tightly focused (the rigorous definition of such a distribution is provided in [9]) and u(g, δt) = u f (g, δt) from (36). The propagated solution would be, which is a specialized case of propagation of probability density using transition probabilities in a continuous-time Markovian process. We can define the group-theoretic mean µ(t) and covariance Σ(t) using the following equations [17,33]: where [log(µ(t) −1 • g) ∨ ] defines a distance from the mean. Note that G [log g] ∨ u(g, t)dg is usually only approximately equal to log(µ(t)) ∨ [8]. Similar definitions of mean and covariance are constructed in [37][38][39] and alternative definitions of an algebraic covariance have also been discussed in [17,40,41].
It is possible to obtain expressions for Σ(t) and µ(t) without explicitly solving for u(g, t) using the concepts of mean and covariance propagation developed in [6][7][8][9][10][11][12]16] by using the relation in (38). For X, Y, Z ∈ G, we have the Baker-Campbell-Hausdorff (BCH) formula: where [8,11], Here, [X, Y] = XY − YX, is the Lie bracket in G. Since x = log(e X ) ∨ and making use of the fact that [X, Y] ∨ = [ad(X)]y we have, where the expansion in (41) (39) and (40), they give rise to products of covariances or higher moments, which are assumed to be negligible. Using the definitions of mean and covariance in (39) and (40) for a convolved probability density u(g, t + δt) = u(g, t) * u(g, δt), we obtain the following expressions for the mean and covariance of the convolved distribution correct to second order: where F captures the second order propagation in the covariance. Its exact form is derived for SE (3) in [8] and more generally for unimodular groups in [11,33]. In this section, we concentrate on first order propagation since it admits a closed-form solution. Successive compositions of µ(δt) using (43) provides the following closed form solution for the mean, which is obtained by integrating the exponent over time [12] (assuming that the initial condition is the identity element of the group): where m(τ) is defined in (35). Setting the initial condition to be the identity element does not lead to a loss of generality since the solution at any other initial condition can be obtained by convolving the fundamental solution (37) with a group Dirac delta function at the initial condition. Using Equation (41) we see that a successive composition becomes an integration in the exponent when the Lie bracket [ t 0 m(τ) ∧ dτ, m(δt) ∧ δt] = 0 [7]; this holds in all examples considered in this paper since m ∧ is time-independent. Approximating Equation (44) to first order and recursively applying Equation (44) by discretising a domain from 0 to t into n segments with step-size δt and taking the limits n → ∞ and δt → 0, we obtain the following integral [12]: Then, an approximate solution to (35) can be constructed as, where u 0 (g) is the initial condition. While initial conditions such as that for European put options [22] can also be considered, we focus on a Dirac delta initial condition centred at the group identity, u 0 (g) = δ G (g) since the solution for any other initial condition can be constructed by a convolution using (47). Moreover, the initial condition should be such that u 0 (g) is square-integrable in G.

Mean and Covariance Propagation for Diffusion Processes on Gl + (1)
In this subsection, we apply the propagation technique to (23). An initial condition of u 0 (a) = δ G (a) = aδ(a − 1) is set (here δ(x) is the Euclidean Dirac delta distribution defined for x ∈ R + ). Using (45) and (46), we have µ(t ) = exp(−(r − σ 2 /2)t ) and Σ(t ) = σ 2 t . The expressions are in terms of reversed time t = −t. Substituting these expressions into (47) and noting that this is the fundamental solution u f (a, t ), we have, and the corresponding solution for the one-asset Black-Scholes Equation (9) The solution obtained from covariance propagation exactly matches the standard analytical solution obtained after applying a logarithmic transformation x = log a to (9), as described in Appendix A.

Mean and Covariance Propagation for Diffusion Processes on Gl
We now apply the propagation technique to (28) with an initial condition of Then, the mean propagation equation would be, and Σ(t ) = Dt . Substituting these expressions in (47), we obtain the fundamental solution u f (g(a, b), t ) as, for l 1 (t ) = log a + (r − σ 2 1 /2)t and l 2 (t ) = log b + (r − σ 2 2 /2)t . Propagation yields a solution that is equivalent to the standard solution obtained after applying the logarithmic transformations x = log a and y = log b to (12) as described in Appendix A. In summary, mean and covariance propagation applied to the one-asset and two-asset Black-Scholes models yields results that exactly match analytical solutions.

Mean and Covariance Propagation Requires Unimodularity of the Lie Group
Both GL + (1) and GL + (1) × GL + (1) are unimodular Lie groups. That is, the Haar measure dg for a group G is bi-invariant and is the only 'natural' measure (upto an arbitrary scaling by a constant) for which the following holds, for a fixed g 0 ∈ G. However A f f + (1) is not a unimodular group, and one can define a right-invariant as well as a left-invariant Haar measure (but no bi-invariant measure exists). The theory of mean and covariance propagation in this section implicitly relies on the group being unimodular; this ensures that there is a unique bi-invariant Haar measure with respect to which a probability density function can be defined. If the group was not unimodular, a probability density defined with respect to one measure may not be a density with respect to the other. Therefore, to apply mean and covariance propagation to obtain the solution for (34), one approach would be to map the diffusion on the group to a diffusion on a related but unimodular group. In the next section, we show that that cotangent bundle of a group, when equipped with semidirect product group operation, is unimodular. This property will be exploited thereon by reinterpreting (34) as a diffusion on the cotangent bundle of the affine group.

Unimodularity of the Cotangent Bundle Group
Let G be an N dimensional matrix Lie group and the Lie algebra G is the tangent space of the group at identity with N basis vectors, represented by E i for i = 1, · · · , N. The tangent bundle TG can be constructed and equipped with a group operation as, Similarly, the cotangent bundle, which is the dual space to the tangent bundle is TG * , and can be equipped with the group operation as, where there exists a bijection between G and R N and similarly between G * and R N . Note that the two expressions indicate that the tangent and cotangent bundles have been endowed with a semidirect product. For an element g ∈ G and X ∈ G, the corresponding element in the tangent bundle would be (g, X) and the group operation in the tangent bundle will be, (g 1 , X 1 ) (g 2 , X 2 ) = (g 1 • g 2 , Ad(g 1 )X 2 + X 1 ) where • is the group operation in G and Ad(g 1 ) is the adjoint representation of g 1 ∈ G where Ad(g) is defined such that Ad(g)X = gXg −1 for X ∈ G. If we were to represent Ad(g) as matrix, we have [Ad(g)] and the coadjoint representation would be [Ad(g)] −T . This is because for X ∈ G and Y ∈ G * , we can construct a bijection such that X ∨ = x and Y ∨ = y where x, y ∈ R N . Note that the ∨ operator for elements from G and elements from G * are different but its usage will be clear from the object to which it is applied to. The inner product is defined to be (X, Y) = x T y = k, k ∈ R. If we use the adjoint representation of the group, a typical element of the tangent space would be Ad(g)X and where (Ad(g) −T Y) ∨ = [Ad(g)] −T y for Y ∈ G * and Y ∨ = y ∈ R N . Tangent and cotangent bundles, endowed with these operations have been used before in [42]. A matrix representation can be constructed for both the tangent and cotangent bundles incorporating the semidirect product as, and, The dimension of the cotangent and tangent bundle group would be 2N provided that the adjoint representation Ad(G) is N-dimensional. This may not hold in certain cases, such as for GL(1), where Ad(G) is zero dimensional. More generally, these constructions will result in a 2N-dimensional semidirect product group if the group G has a trivial center. For the affine group A f f + (1), whose cotangent bundle group will be considered in the paper, this is true. However for G = GL(N), the center is one-dimensional (consisting of scalar multiples of the identity matrix I N ) and Ad(G) would be N 2 − 1 dimensional. We postpone the discussion of constructing the cotangent and tangent bundle groups for groups with non-trivial centers to a future paper.

Properties of the Adjoint Operator
The proofs of unimodularity will depend on the properties of the adjoint operator in the tangent and cotangent bundle group. We present a few properties that will be useful in constructing the proofs in the later subsections.
we can convert the matrix representation of the operator into its coordinate free form through the ∧ operator. For x ∈ R N where x ∧ = X ∈ G, we have, The right-hand side can be simplified as, maps the Lie algebra of the tangent bundle to R 2N whereas the ∨ operator acting on E i maps G to R N . The overloading of the ∨ symbol should not lead to confusion since the object on which it acts will determine its interpretation.

Lie Algebra of TG and TG *
The basis vectors spanning the Lie algebra of the tangent bundle can be obtained by differentiating an element representation (57) with respect to each of the N variables parameterizing it, at identity.
The Lie algebra basis can be obtained for the parameters parameterizing the group G as well as for the parameters parameterizing G as, Note that e ∧ j−N = E j−N where E j−N is a basis of G. ForẼ i defined to be a basis element of the Lie algebra of (TG, ), the ∨ operator is defined such thatẼ ∨ i = e i where this time i = 1, · · · , 2N and the ∨ operator maps from this Lie algebra to R 2N . A similar definition can also be created for the Lie algebra of the cotangent bundle, (TG * , ). In this case, the basis of the Lie algebra can be deduced as, Here, E i ∈ G but e * i ∈ (G * ) ∨ , the latter of which is related to the cotangent space at identity. In general for an element of the Lie algebra of G G, the vee operator is defined such that, where X ∈ G and y ∈ G ∨ . A similar definition exists for the cotangent bundle, where X ∈ G and y ∈ (G * ) ∨ . The specific ∨ used will be clear from the context in terms of which argument it is applied to. The Lie algebra for the (co-)tangent bundle group defined this way is only 2N-dimensional for groups with trivial centers. Proof. This can be proven by considering the adjoint representation of an element within the tangent bundle. The adjoint, Ad(h)Ẽ i where h ∈ G G, is given through matrix notation as,

Unimodularity of TG and TG
Solving this by substituting the two different forms ofẼ i obtained from (61) and using the relationship that [Ad If G has a non-trivial center Z(G), then for a k ∈ Z(G), one can decompose g = h(q 1 , · · · , q m ) • k(q m+1 , · · · , q N ) such that h and k commute. This is true since one can partition G such that any g ∈ G can be decomposed into g = h • k where h ∈ F G/Z(G) , k ∈ Z(G), and F G/Z(G) ⊂ G is the 'fundamental domain' [33] associated with the quotient group G/Z(G). The dimension of Ad(G) will then be dim(G) − dim(Z(G)) = m. Therefore, the current construction of using Ad(G) to represent the group G will no longer be faithful when dim(Z(G)) = 0 implying that since the theorem and proof are restricted to the special construction Ad(G) G, they only hold in their current form for the case where G has a trivial center and Ad(G) G is 2N dimensional. Theorem 2. The cotangent bundle group (TG * , ) = G G * of an N-dimensional Lie group with a trivial center is always a 2N-dimensional unimodular Lie group independent of the unimodularity of G.
Proof. We construct a proof in a very similar way as that in Theorem 1. Here, we obtain the adjoint representation of the cotangent bundle in the following form for h ∈ G G * , Since det [Ad(h)] = 1, this ensures that the cotangent bundle group is unimodular, independent of the unimodularity of G. The restriction to trivial centers follows from the fact that Ad(G) −T is only N-dimensional for groups with trivial centers and therefore Ad −T (G) G will only have the same number of dimensions as the cotangent bundle with such a restriction.

Option Price Evolution with Coupled Assets as a Diffusion Process on the Cotangent Bundle of the Affine Group
Now we apply the results from the previous section to construct the cotangent bundle of the affine group. The motivation is to re-express Equation (34) in terms of the Lie derivatives of the cotangent bundle of A f f + (1). An element h from the cotangent bundle of the affine group, h ∈ Ad(A f f + (1)) −T R 2 , can be expressed in the following form using (31) and (58), where (a, b, x, y) ∈ R + × R 3 . Using (62), the orthonormal basis of the Lie algebra of this group, expressed in matrices, is , and hence we can construct a bijection from the Lie algebra to R 4 such thatẼ ∨ i = e i for i = 1, 2, 3, 4. We argue that although the number of dimensions is doubled when we consider the cotangent bundle, this has the benefit of making the group unimodular. To see this, consider the left and right Jacobians of the cotangent bundle group: Then, we observe that detJ R = detJ L = 1/a, implying unimodularity. The adjoint of the affine cotangent bundle [Ad(h)] is, where det [Ad(h)] = 1, as expected. For completeness, [ad( The left and right Lie derivatives are, where f (h(a, b, x, y)) =f (a, b, x, y). We now observe from (33) and (73) thatẼ R 1 = E R 1 andẼ R 2 = E R 2 , which allows us to rewrite (34) in terms of the cotangent bundle group Lie derivatives as, which is a degenerate diffusion over the affine cotangent bundle group. We denote the 'master' function as, u = u(h(a, b, x, y), t), andũ =ũ(a, b, t) as the marginalised version of this 'master' function marginalised over the variables x and y.

Mean and Covariance Propagation for Diffusion Processes in the Cotangent Bundle of A f f + (1)
Since Ad(A f f + (1)) −T R 2 is unimodular, the theory of mean and covariance propagation can be applied directly from Section 4 to solve (74). The initial condition is u 0 (h(a, b, x, y)) = δ G (h) = aδ(a − 1)δ(b − 1)δ(x)δ(y) where µ(0) would be the 2 × 2 identity matrix and the mean propagation equation would be, for k 1 = r − σ 2 1 /2 and which is a singular matrix, implying that this is a degenerate diffusion over the cotangent bundle group. The covariance propagation equation is then, such that, where κ = k 2 /k 1 , I 1 = (e k 1 t − 1)/k 1 and I 2 = (e 2k 1 t − 1)/(2k 1 ). Noticeably Σ(t ) is singular, but we treat this issue by writing, where I 2 is the 2×2 identity matrix. Then, u f (h(a, b, x, y), t ) can be written as: For the affine cotangent bundle, it is possible to obtain a closed-form expression for log(µ(t ) −1 • h) ∨ . Any µ(t ) ∈ Ad(A f f + (1)) −T R 2 can be expressed in the following form, using the expression of h in (68) and where a and b are functions of t . Then, one has, where, The derivations of these expressions are provided in Appendix D. For convenience of notation, we express the first two components of log(µ(t ) −1 • h) ∨ as log(µ(t ) −1 • h) ∨ and the third and fourth components with F so that we have, Substituting this form into (80) and taking the limit → 0, we obtain, where δ(· · · ) is the Dirac delta function. Additionally, we re-express δ(F 1 )δ(F 2 ) in terms of δ(x)δ(y) as, which can be derived by calculating the Jacobian determinant of the system of Equations in (83) and (84). We note that the variables x and y are extraneous to the original problem and we aim to marginalise the distribution u f (h(a, b, x, y), t ) from (80) toũ f (a, b, t ) as, Using (87) and (86), we have, which is the fundamental solution to (74).

Normalisation of Probability Distribution Functions on the Affine Cotangent Bundle Group
For the affine cotangent bundle group, we can define a general Gaussian as, for h ∈ A f f + (1) R 2 , where a general normalisation factor c(Σ), which is a function of the covariance Σ, is used. In the previous section, we used a factor of the form, c(Σ ) = |det Σ −1 | 1/2 /(2π) N/2 for N = 4 in the affine cotangent bundle group (prior to marginalisation; see Equation (80)). However, this is only correct for small covariances and assuming that the determinant of the group Jacobian (70) is sufficiently close to 1. The goal of this section is to derive a higher order correction to this normalisation factor for the affine cotangent bundle group. To do so, we first convert to exponential coordinates so that for an arbitrary element g ∈ A f f + (1) R 2 , where q = [α, β, γ, φ] T is the exponential coordinate parameterization of the group. Then, the right Jacobian of the group defined in (16) for this coordinate system would be (1/a)∂a/∂α 0 0 0 (1/a)∂b/∂α (1/a)∂b/∂β 0 0 where the use of (a, b, x, y) makes contact with the earlier parameterization in (68); expressing α, β, γ and φ as functions of (a, b, x, y) we have, where since the cotangent bundle group is unimodular, we have dropped the R subscript noting that det J R = det J L independent of the parameterization. Expanding det J exp to small α we have, We can also write this relation as, for Now if we consider the integral of (90) over the cotangent bundle group, and make the substitution g = µ −1 • h, we can rewrite (100) as, by making use of the unimodularity of the cotangent bundle and noting that log(g(q)) ∨ = q for a parameterization in terms of the exponential coordinates q. Thus, where we use the relationship for the determinant of the Jacobian in (99). Equating (102) to the number 1 and noting that the integral evaluates to (2π) In the context of the degenerate diffusions on the affine cotangent bundle group that arise in the coupled asset model, we have, We also note that the sign of the non-zero element ofW is negative. Whereas, all terms ofΣ −1 are positive. This suggests that at a sufficiently large covariance, it is possible that c(Σ) = 0 as a consequence of the current approximation. While this was not observed for the small covariances used in this paper, this suggests an opportunity to use a different method of approximating the integral, by using the following general relation for q ∈ R N [29,43]: where tr(· · · ) is the trace operator. Nevertheless, in this paper we do not pursue this approximation and instead use the result in (104) assuming that the eigenvalues ofΣ −1 are sufficiently large (which is the case for the range of parameters used in the numerical simulations in the subsequent sections).

Numerical Results for Option Price Evolution with Coupled Assets
We solve the PDE in (74) by four methods: (1) Finite difference method (implicit and explicit), (2) first order propagation, (3) second order propagation and (4) Euler-Maruyama integration of the underlying stochastic differential equations. The theory of second order covariance propagation has not yet been introduced for A f f + (1) R 2 , and will also be described in this section. All simulations were performed using MATLAB R2019a on a 2.7 GHz Dual-Core Intel Core i5 processor. The CPU time to run 10 time steps of second order propagation, explicit finite difference method and implicit finite difference method was 18.94 s, 0.99 s and 11.55 s (assuming that the finite difference matrices are constructed before-hand), respectively; however, if the matrix logarithms are evaluated analytically rather than numerically-which is possible for the affine cotangent bundle group (82)-the CPU time for first order propagation reduces to 0.30 s for 10 time steps.

Finite Difference Method
The 2D finite difference scheme was implemented on a rectangular a-b grid with second-order accuracy. Explicit and implicit schemes were constructed, using the forward and backward Euler scheme, respectively. A time-step of 10 −6 units was used for the explicit scheme and a time-step of 10 −3 units for the implicit scheme. A smaller time-step was used for the explicit method to ensure stability. A grid spacing of approximately 4.92 × 10 −3 units along a and 5.24 × 10 −3 units along b was used. The simulation domain was chosen to best capture the distribution while ensuring that the boundaries B were sufficiently far from the mode of the distribution. This was to ensure that the Dirichlet boundary condition ofũ(a, b)| B = 0 could be set. A Dirac delta initial condition was approximated using a circular Gaussian with a standard deviation equal to twice the grid spacing along a.

Euler-Maruyama Integration of Underlying Stochastic Differential Equations
It is possible to convert a Fokker-Planck equation over an N-dimensional Lie group to a stochastic differential equation in the Lie algebra. A Fokker-Planck equation for a probability density of the form u(g, t) in the group G with a drift of m(t) and diffusivity D(t) is, whereẼ R i is a right directional Lie derivative in G. From [17,29], we recognise that the stochastic differential equation described by this Fokker-Planck equation on the Lie algebra G would be of the form, where dW is a vector of increments of N uncorrelated Wiener processes, W 1 , · · · , W N , corresponding to random draws from a Gaussian with zero mean and variance dt, and BB T = D. Equation (108) can be interpreted either as a Stratonovich or Itô equation since the diffusion term B is independent ofx. Sincex ∈ G ∨ , this process occurs in the Lie algebra of G. This process can then be injected on to the group [17] as, thereby defining a stochastic process on G. If we were to parameterize the group elements with a vector q ∈ R N such that dx = [J R (q)]dq where [J R (q)] is the right Jacobian matrix, we obtain a Stratonovich stochastic differential equation in the parameter space R N as, where emphasises that this is a Stratonovich equation. For the cotangent bundle group of the affine group, dq = [da, db, dx, dy] T and the form of [J R ] is provided in (70). Matching (107) with (74), we obtain, and, Using the form of [J R ] from (70) and the expressions for m and B from (111) and (112), we have the following Stratonovich stochastic differential equations for a, b, x, y: where dx = 0 and dy = 0 highlight the degenerate nature of the stochastic differential equation. However, to implement an Euler-Maruyama integration in parameter space, we require the stochastic differential equations to be expressed as Itô equations; the distinction between the Itô and Stratonovich form of a stochastic differential equation is especially important here since the diffusivity is now a function of a. The Itô version of the equations is thus, where the terms 1 2 (B 2 11 + B 2 12 ) and 1 2 (B 21 B 11 + B 22 B 12 ) correct for the drift. Equivalently, instead of solving (114), it is also possible to obtain an evolution in parameter space by projecting the stochastic evolution of g(a(t), b(t), x(t), y(t)) in the cotangent bundle (109) on to the space of parameters (a, b, x, y). This is because for any g(t) ∈ A f f + (1) R 2 we have, which determines a unique point in a stochastic trajectory in parameter space (a(t), b(t), x(t), y(t)).
Due to the degenerate nature of the diffusion process, there is no evolution in the parameters (x, y). This method of obtaining the stochastic process in parameter space, although equivalent in principle to that obtained by integrating (114), is henceforth referred to as an Itô-Gangolli method since it makes use of the McKean-Gangolli injection [17] and numerically solves the stochastic differential equations in the Lie algebra of the group (108) by an Euler-Maruyama integration. The group-theoretic covariance and mean of the probability density function u(g, t) were deduced from the ensemble generated by (109) on the group using the methods in [8,33]. That is, the discrete version of (39) can be obtained by setting the sampled probability density to be u s (g, t) = 1 N ∑ N s i=1 δ G (gg −1 i ) for a total of N s samples. Here, δ G (g) is the group Dirac delta function. A similar substitution in (40) gives the sample covariance. Thus, the sampled mean µ s (t) and covariance Σ s (t) are, In the context of the Euler-Maruyama integration, N s is the total number of sample paths and the averaging is performed at each time slice t. The value of g i (t) is obtained from the evolution process in G described in (109). One would expect that for large N s , µ s (t) and Σ s (t) will approximate the corresponding mean and covariance of the solution of (74) to the extent that the exact solution remains a Gaussian on the cotangent bundle group. Hence, the results from the Euler-Maruyama integration of a large number of sample paths can serve as a baseline truth to compare results against. Note that (116) needs to be solved iteratively by beginning with the following guess solution, and the mean at the (j + 1) th iteration can be computed using, where we observe that 1 ] quantifies an error that goes to zero once the mean converges. This converged value is substituted as µ s (t) in (117) to obtain the sampled covariance. In the simulations, N s = 30,000 and a time-step of 10 −3 units was used. After generating the ensemble of sample paths, a continuous probability distribution was created in (a, b) space by kernel density estimation using a Gaussian kernel. This kernel density estimated probability distribution was used as a baseline to compare against finite difference and propagation solutions (shown later in the contour plots of Figures 3 and Figures 6).

Second Order Propagation for
The mean propagates by (43) where µ(δt) = exp(m ∧ δt) for m = −[k 1 , k 2 ] T using the definitions in (75). Using (72) and (42) and substituting it in the definition of covariance in (40), we have the following expression for second-order covariance propagation [11,33]: where, Much like the case for first order propagation, only the first 2 × 2 subspace of Σ(t ) is non-zero, corresponding toΣ(t ). A time-step of 10 −3 units was used in second order propagation. Figure 1 depicts the convergence of the error in Σ(t ) with respect to time-step. The baseline truth in this case was Σ s (t ) sampled from the ensemble of paths generated in the Euler-Maruyama integration (117).

Results
For the numerical simulations, we consider the PDE in (74) with two sets of parameters, S l : {σ 1 = 0.1, σ 2 = 0.05, r = 3, µ 1 = 1, µ 2 = 2} and S h : {σ 1 = 1, σ 2 = 0.5, r = 3, µ 1 = 1, µ 2 = 2}; the set of parameters S l describes a scenario where diffusion of asset price is low and S h describes a scenario where the diffusion is relatively large. Since the values of µ 1 , µ 2 and r are same in both cases, the values of σ 1 and σ 2 will be used to distinguish these two sets. The constraint σ 2 /σ 1 = (µ 2 − r)/(µ 1 − r) was used to fix the ratio σ 2 /σ 1 in both cases. Additionally, the initial condition is a Dirac delta distribution at the identity element of the cotangent bundle group of A f f + (1).
Firstly, we show the differences between first and second order group-theoretic covariance propagation for the large diffusion scenario (S h ) in Figure 2. The error in covariance is plotted as a relative deviation using the following error metric, where || · · · || denotes the Frobenius matrix norm, Σ s (t ) is the sampled covariance from (117) and Σ(t ) is either the first or second order propagated covariance. The relative error in mean was evaluated in a similar fashion: where µ s (t ) is the sampled mean from (116) and µ(t ) is either the first or second order propagated mean. For the set of parameters describing 'large diffusion' (parameter set S h ) the relative error in the mean was approximately 0.1 to 0.2% and first order and second order mean propagation were indistinguishable. Furthermore, in the case of 'small diffusion' (parameter set S l ), the error in covariance was approximately 0.4 to 1.2% and with minimal difference between first order and second order propagated results; finally, the relative error in mean for this set of parameters was in the order of 0.02% and again with no difference between first and second order propagation. Only approximate values are given since this range of error is within the variability of the sampled mean and covariance itself, which are used as the baseline. We now proceed to compare the results from first order and second order propagation with those from finite-difference methods, relative to the probability density function obtained from the Itô-Gangolli method used to indirectly solve the stochastic differential equations in (114). The probability distribution corresponding to the ensemble of points (a, b) ∈ R + × R in parameter space was considered as the ground truth for the following numerical studies. One can then estimate the mean and covariance of this ground truth at time t by, for N s sample paths and q i = [a i (t), b i (t)] T . Note that these are the expressions for the mean and covariance defined in Euclidean space (see Appendix C). It is important to emphasise that here we are not comparing on the basis of the group-theoretic mean and covariance but rather on the basis of a mean and covariance defined in R + × R 3 . Since the solution from propagation (followed by marginalisation) or finite difference methods would yield a probability density over the affine cotangent bundle group, it is important to convert these results to an equivalent probability density function on parameter space. That is, if f (g, t) is a probability distribution on the affine cotangent bundle, andf (q, t) =f (q, t) det J(q) would be a probability distribution in parameter space R + × R 3 . In the special case of degenerate diffusion for the coupled asset model,f (q, t) is a probability distribution in the Euclidean half-space R + × R. The mean and covariance of such a distribution can then be evaluated as, For q = [a, b] T ; we also have det J(q) = 1/a from (70). The mean and covariance defined this way is evaluated for the finite difference and propagation solution (based on (89) but using the higher-order normalisation factor from (104)) and compared against the sampled mean and covariance obtained from the ground truth in (124,125). 7.4.1. Small Diffusion: σ 1 = 0.1, σ 2 = 0.05 Contour plots were generated for first order propagation, second order propagation, explicit and implicit finite difference methods. The baseline probability density was smoothed by a kernel density estimation procedure using a Gaussian kernel, which was used as the ground truth for the contour plots to qualitatively assess the shape of the distribution. These are shown at 300 time steps into the simulation in Figure 3.
The poor performance of the finite difference solution relative to the propagation solution is also evident in Figure 4. In this figure, the relative error is measured with respect to the sampled covariance (125) and calculated using the same formula in (122), but where Σ(t ) is the covariance in terms of parameters (a, b) from (128).  The relatively poor performance of the finite difference solution is because the covariance is very small (but not zero) such that one observes spurious oscillations in the finite difference solution (see Figure 5). Negative values ofũ f (a, b, t) are artefacts of the discretisation. In simulating the evolution of distributions with low covariance, a finite difference solution requires a very fine mesh near the mode of the distribution to avoid such artefacts whereas a propagated solution requires no such discretisation in (a, b) space and automatically avoids these issues. Nevertheless, the relative error in mean, measured through a Euclidean norm, was approximately 0.2 to 0.4% for both propagation and finite difference simulations. Contour plots were generated for first order propagation, second order propagation, explicit and implicit finite difference methods. These are shown at 300 time steps into the simulation in Figure 6. Figures 7 and 8 show the mean and covariance evaluated using (127) and (128) compared relative to (124) and (125). The error in covariance was measured by a relative Frobenius norm and the relative error in mean was measured using a Euclidean norm defined as, for = µ − µ s and µ s is the sampled mean (124).   Mean and covariance propagation when applied to the coupled asset model tend to yield results with lower error when the covariances are small, in line with the original assumptions made to derive the technique. Additionally, covariance propagation is more suitable in dealing with Dirac delta initial conditions than a finite difference method. Another advantage of the propagation technique as opposed to a finite difference solution is that there is no grid involved: one effectively has a solution that is not discretized in the (a, b) domain and only requires a temporal discretization (for the second order propagation) or no discretization at all (for the first order propagation where the solution can be written in closed-form). Finally, a synergy of the two methods would be useful in spanning a broader range of covariances than either method can handle on its own.

Backward Compatibility of Propagation on the Cotangent Bundle with the One-Asset Black-Scholes Equation
In this final section, we consider reframing the one-asset Black-Scholes equation as a diffusion process on the affine cotangent bundle group. This is possible because the Lie derivative of GL + (1) (22) can also be represented byẼ r 1 from (73). Hence, we can write (24) in terms of the affine cotangent bundle group Lie derivatives as, By first order propagation, we know that the mean is (assuming an identity initial condition), The diffusion matrix is now, and, since [Ad(µ(t ) −1 )] and D are both diagonal. Similar results can also be obtained numerically using second order propagation for both mean and covariance. Following the discussion in Section 6.1, we construct a Gaussian on the cotangent bundle of the form, where using the form of µ(t ) in (81) we have, so that, F 1 = x + a a − a + aa log(a/a ) (a − a ) 2 by (137) and log(µ(t ) −1 • h) ∨ = log(µ(t ) −1 • h) ∨ · e 1 . Then we can show that, by constructing a Jacobian matrix J from the system of equations in (136)-(138) and using its determinant det J to normalise the Dirac delta function. Marginalising over the variables b, x, y, we obtain the solution for the one-asset Black-Scholes equation given by propagation on the affine cotangent bundle as, u f (h(a, b, x, y), t ) = 1 σ √ 2πt exp −1 2σ 2 t [(r − σ 2 2 )t + log a] 2 1 aa a − a log(a /a) 2 . (140) It is important to note that the solution in (140) differs from (48) by the Jacobian factor 1/det J and therefore is not an exact solution. However we see that for small values of = (a − a )/a , 1/det J ≈ 1 to O( 2 ). In this case, it is also possible to compare the propagated result with the analytical solution for the Black-Scholes equation. We make this comparison for S l : {r = 3, σ = 0.5} representing a small diffusion scenario and S h : {r = 3, σ = 1} representing a large diffusion scenario, at t = 0.30 and show the plots in Figures 9 and 10. Moreover, the higher-order normalisation factor in (104) is used to normalise the propagated result in (140) instead of 1/σ √ 2πt . We see a closer match with the analytical solution for the small diffusion scenario.

Conclusions
Reframing PDEs in R N as diffusion processes on Lie groups offers an alternative method to solve PDEs by using mean and covariance propagation techniques developed previously in the context of Fokker-Planck equations on Lie groups [6][7][8][9][10]12]. In the case of asset dynamics from mathematical finance, the method yields the exact solution for the one-asset and two-asset problems by matching the Lie derivatives of the one-asset and two-asset Black-Scholes equations with the Lie derivatives of GL + (1) and GL + (1) × GL + (1), respectively; this trivially reduces to the logarithmic coordinate transformation that converts these equations to heat equations.
While using the apparatus of mean/covariance propagation on Lie groups is undue for the one-asset and two-asset Black-Scholes equations, the matching is especially useful for the model of option price evolution under coupled asset dynamics introduced in the paper where the logarithmic coordinate transformation characteristic of the one-asset and two-asset Black-Scholes PDE can no longer be applied. Instead, we solve the equation by matching the derivatives with the Lie derivatives of A f f + (1). We provide proofs of the unimodularity of the cotangent bundle of a Lie group, and exploit this property to perform a mean/covariance propagation on the cotangent bundle of A f f + (1).
The mathematical apparatus developed can be applied to different PDEs in mathematical finance, such as those arising from stochastic volatility models or other forms of asset coupling in multi-asset models, and to linear convection-diffusion equations in transport theory, to name a few. Due to the unimodularity of the cotangent bundle group, the Lie group to which the derivatives are matched with need not be unimodular. Additionally, subsequent research can also be directed towards extending the cotangent bundle group construction presented here to groups with non-trivial centers, and to a general stability analysis of propagation schemes.
which is the Green's function for the one-asset Black-Scholes equation.

Appendix B. The Solution of a Diffusion Equation on a Unimodular Lie Group Is a Probability Density
Consider the diffusion equation for u(g, t) on an N-dimensional unimodular Lie group G of the form, ∂u ∂t If we let P = G u dg and integrating both sides of (A7) with respect to the bi-invariant Haar measure dg, we have, If u and E R i u are 'well-behaved' in that their absolute values and the square of the functions integrate to a finite value over the group then, by integration by parts, both integrals on the right-hand-side vanish and we are left with ∂P/∂t = 0 [17], which is the statement that the total probability (correct to a constant normalisation factor) is conserved. That is, if the initial condition u(g, 0) is a probability density, u(g, t) will also be a probability density.

Appendix C. Mean and Covariance in R N
For a probability density functionũ(q, t) where q ∈ R N , the mean µ(t) and covariance Σ(t) are defined as, R N (q − µ(t))ũ(q, t) dq= 0, (A9)