Special Issue in honor of John W. Neuberger Electronic Journal of Differential Equations, Special Issue 02 (2023), pp. 161–173. ISSN: 1072-6691. URL: https://ejde.math.txstate.edu or https://ejde.math.unt.edu INVERSE VOLATILITY PROBLEM FOR CURRENCY OPTIONS IAN KNOWLES, SUNDAR TAMANG Abstract. In transactions associated with future-oriented financial instru- ments, such as options, a huge amount of data is available buried inside of which is the market’s best guess as to what the future holds. We consider here the possibility of extracting future foreign exchange volatility informa- tion from foreign exchange option data with the aid of a new computational inverse algorithm using minimization of a convex functional. 1. Introduction Financial entities such as stocks, bonds and even foreign exchange (FX) rates, are often modeled by the stochastic differential equation dSt St = m(St, t)dt+ σ(St, t)dBt, (1.1) where, for each time t, St(ω) is a random variable representing the price of the financial entity for the sample path ω, m is the drift, which relates to the “trend” of the entity, σ is the volatility (“wobble”), and Bt(ω) is the Brownian motion stochastic process used to model this log-normal randomness. Financial derivatives are contracts that derive their value from such an underly- ing entity. In particular, a European call option on a stock, is a derivative financial instrument that, when purchased, gives the holder the right but not the obligation to buy a designated number of stock shares at a pre-agreed price (the strike price) on a specified date, and a European put option provides the right to sell under the same circumstances. In their seminal work [2] Black and Scholes showed that, for dividend-free stocks, the arbitrage-free price v(S, t;K,T ) of a European option contract satisfies the deterministic partial differential equation ∂v ∂t + 1 2 σ2(S, t)S2 ∂ 2v ∂S2 + µs ∂v ∂S − rv = 0 (1.2) in time t and the value S of the underlying asset, where σ(S, t) is the volatility, K is the strike price, T is the expiration time, µ is the risk-neutral drift, and r is the short-term risk-free interest rate. 2020 Mathematics Subject Classification. 35K10, 65N21. Key words and phrases. Inverse volatility; European FX options; convex functional. ©2023 This work is licensed under a CC BY 4.0 license. Published March 27, 2023. 161 162 I. KNOWLES, S. TAMANG EJDE/SI/02 Analogously, in the case of FX options, a European call option on a foreign exchange rate (currency) pair, is a derivative financial instrument that, when pur- chased, gives the holder the right but not the obligation to buy a designated amount of money, denominated in one currency, in another currency at a pre-agreed ex- change rate on a specified date, and a European put option provides the right to sell under the same circumstances. Assume that S is the price of the deliverable currency (domestic units per foreign unit), K is the exercise price of the option (domestic units per foreign unit), t is the current time, u(S, t;K,T ) is the current price of an FX call option (domestic units per foreign unit), rD is the domestic (riskless) interest rate, rF is the foreign (risk- less) interest rate, and σ is volatility of the spot currency price. The corresponding Black Scholes partial differential equation ∂u ∂t + 1 2 σ2(S, t)S2 ∂ 2u ∂S2 + (rD − rF )S ∂u ∂S − rDu = 0 (1.3) for this case was provided by Garman and Kohlhagen in their classic 1983 paper [8], where for consistency we have replaced the time to maturity T in [8, equ. 6] with the current time t, and we have allowed the volatility to vary with S and t. We note also the final condition u(S, T ) = { S −K, if S > K 0, if S ≤ K, and the boundary condition u(0, t) = 0, 0 ≤ t ≤ T hold for this call option. Now, financial data specifying the market price of an option is readily available in quantity at various strike values K around the current price of the underlying (the “spot” price), and for values of the maturity T up to around twelve months or so into the future. Given that the computer projections currently used by financial analysts are at best only valid for a rather limited time into the future, one is led quite naturally to the so-called inverse volatility problem: determine a market- inspired estimate of the future volatility function σ(S, t) from a knowledge of current market prices of options with different strikes and future maturities. In the case of stocks the solution of this problem generally goes as follows. The value v of an option contract depends crucially on the exercise (strike) price K, and the maturity date T of the contract. In 1994 Bruno Dupire [7] noticed that the function v(S, t;K,T ) satisfies the “dual” Black-Scholes equation ∂v ∂T − 1 2 K2σ(K,T ) ∂2v ∂K2 + µK ∂v ∂K − (µ− r)v = 0, (1.4) known also as the Dupire equation. A number of approaches have been proposed subsequent to the appearance of [7]. Minimization methods using regularized least- squares fitting have been proposed in [1, 3, 16]; the possible presence of spurious local minima is always an issue here. An integral equation approach is presented in [4, 5], where convergence problems are possible given the underlying ill-posedness, and in [6, 10] linearization of the inherently non-linear inverse problem is discussed. In [12] the volatility is recovered by the minimization of a convex functional, thereby avoiding the spurious local minima problem that plagues generic least square meth- ods. For FX options, observing that (1.3) may be obtained from (1.2) by replacing µ with rD − rF and r with rD. On using the same replacements in (1.4) we then EJDE-2023/SI/02 INVERSE VOLATILITY PROBLEM 163 obtain the Dupire-type equation ∂u ∂T − 1 2 K2σ(K,T ) ∂2u ∂K2 + (rD − rF )K ∂u ∂K + rFu = 0, (1.5) for the price u(S, t;K,T ) of a European FX call option. If, at a certain time t0 (when the spot price of the underlying is S0 = S(t0)), the prices u(S0, t0,K, T ) are known for all strikes K and maturities T then, as was noted first in [7], the volatility σ(K,T ) is uniquely determined in principle from the equation (1.5). But such a formula for σ is of little use in practice, as the market data for u is not only noisy (which would make the estimation of these derivatives highly ill-posed), but even worse, the data is typically both discrete and sparse in the variables T,K. We present here a new variational algorithm for recovering, via the Dupire equa- tion (1.5), the volatility σ(K,T ) from a knowledge of European FX option prices at various strikes and maturities, that is, from a knowledge of u(S0, t0;K,T ) for some time t0 ≤ T . Here, as the variable K represents prices around the current spot price and T ≥ t0 represents future times, we have in effect recovered σ(S, t) for future times t and FX prices S in a neighborhood of the spot price S0. The method used is an adaption, and improvement in the FX case, of the vari- ational approach used in [12] involving the minimization of convex functionals as we discuss below. We begin by outlining the approach used in [12] to set the stage for our new method. 2. Reconstruction of FX volatility Let T0 < T1 < · · · < Tn be maturity times, and for each maturity Ti, 0 ≤ i ≤ n, let Ki,1 < Ki,2 < · · · < Ki,mi be the known associated strike prices. From the known price data u(S0, t0;Ki,j , Ti), 0 ≤ i ≤ n and 1 ≤ j ≤ mi obtained at some time t0, we compute the function u(S0, t0;K,T ) by linear interpolation, and for convenience, we make use of the notation u(S0, t0;K,T ) = u(K,T ). Next, we assume that the volatility is piecewise constant in time, so that, for 1 ≤ i ≤ n, σ = σi(K) over the i-th time sub-interval [Ti−1, Ti]. Fixing i, set σi(K) = σ(K) and wλ(K) = ∫ Ti Ti−1 e−λTu(K,T ) dT, (2.1) where λ > 0 is a parameter. For each such fixed i, 1 ≤ i ≤ n, we now Laplace transform the Dupire equation (1.5) over [Ti−1, Ti] to obtain 0 = ∫ Ti Ti−1 e−λTuT dT − 1 2 K2σ2 ∫ Ti Ti−1 e−λTuKK dT︸ ︷︷ ︸ w′′λ + (rD − rF )K ∫ Ti Ti−1 e−λTuK dT︸ ︷︷ ︸ w′λ +rFwλ, 164 I. KNOWLES, S. TAMANG EJDE/SI/02 where the primes indicate differentiation with respect to K. On integrating the first term by parts we obtain [e−λTu]TiTi−1 + λ ∫ Ti Ti−1 e−λTu dT︸ ︷︷ ︸ wλ −1 2 K2σ2w′′λ + (rD − rF )Kw′λ + rFwλ = 0, and rearranging terms gives −1 2 K2σ2w′′λ + (rD − rF )Kw′λ + (λ+ rF )wλ = −u(K,Ti)e −λTi +u(K,Ti−1)e−λTi−1 . Next, dividing by 1 2K 2σ2 throughout, we see that −(w′′λ − 2(rD − rF ) Kσ2 w′λ) + λ+ rF 1 2K 2σ2 wλ = −u(K,Ti)e −λTi + u(K,Ti−1)e−λTi−1 1 2K 2σ2 . Finally, on multiplying by the integrating factor p(K) = e −2(rD−rF ) ∫K a dk kσ2(k) , (2.2) we now have an equation in Sturm-Liouville form: − (p(K)w′λ)′ + (λ+ rF )q(K)wλ = β(K,λ)q(K), (2.3) where q(K) = ( 2 K2σ2(K) ) p(K), (2.4) β(K,λ) = −u(K,Ti)e −λTi + u(K,Ti−1)e−λTi−1 . (2.5) The idea in [12] is that if one can recover the functions p(K) and q(K) from (2.3), one can then find the volatility σ(K) from the formula σ(K) = √ 2p(K) K2q(K) . (2.6) So, we first focus attention on a variational approach to the recovery of the pair of positive coefficient functions p, q defined on the generic interval a ≤ K ≤ b. It is assumed that we are given the functions wλ(K) for K in [a, b] and all λ > 0. For positive functions r and s also defined on [a, b], let c = (r, s). Define wc,λ(K) to be the solution to the boundary value problem Lc,λwc,λ = −(r(K)w′c,λ)′ + (λ+ rF )s(K)wc,λ = β(K,λ)q(K), wc,λ(a) = wλ(a), wc,λ(b) = wλ(b). (2.7) Let D be the set of all positive function pairs c = (r, s) such that the boundary value problem (2.7) is disconjugate on [a, b], i.e. every non-trivial solution has at most one zero on [a, b]. It is known [9, Theorem 6.1, p. 351] that (2.7) is disconjugate if and only if the boundary value problem (2.7) can always be solved uniquely. It is also known (c.f. [15, Proposition 2.1]) that this set is open and convex in L[a, b]×L[a, b] and L2[a, b]× L2[a, b]. For each λ > 0 define the functional Gλ on the convex set D by Gλ(c) = ∫ b a r(K)(w′2λ − w′2c,λ) + (λ+ rF )s(K)(w2 λ − w2 c,λ) − 2βs(K)(wλ − wc,λ) dK. (2.8) EJDE-2023/SI/02 INVERSE VOLATILITY PROBLEM 165 2.1. Properties of the functional Gλ. The main properties of the functional Gλ are summarized in the following theorem, which is [12, Theorem 3.1]. Theorem 2.1. (a) For any c = (r, s) in D, Gλ(c) = ∫ b a r(w′λ − w′λ,c)2 + (λ+ rF )s(wλ − wc,λ)2. (2.9) (b) Gλ(c) ≥ 0 for all c = (r, s) in D, and Gλ(c)=0 if and only if wλ = wc,λ. (c) The first Gâteaux derivative of Gλ is given by G′λ(r, s)[h1, h2] = ∫ b a (w′2λ − w′2c,λ)︸ ︷︷ ︸ L2 gradient in r h1 + [(λ+ rF )(w2 λ − w2 c,λ)− 2β(wλ − wc,λ)]︸ ︷︷ ︸ L2 gradient in s h2. (2.10) (d) The second Gâteaux derivative of Gλ is given by G′′λ(c)[h, k] = 2(L−1c,λ(e(h)), e(k)), (2.11) where h = (h1, h2) , k = (k1, k2), e(h) = −(h1w ′ λ,c) ′ + (λ+ rF )h2wc,λ − βh2, and (· , · ) denotes the usual inner product in L2[a, b]. As Lc,λ is a positive operator on W 1 0 [a, b], we have from Theorem 2.1(d) that G′′λ(c) ≥ 0 for all c in the convex set D. By [18, Corollary 42.8] the functional Gλ is therefore convex on D. We know from Theorem 2.1(b) that Gλ has a global minimum (zero) at c = (r, s) if and only if wλ = wc,λ. Choose N ≥ 3 positive distinct real numbers λj , 1 ≤ j ≤ N , so that 0 < (λj + rF )T < 2, T ∈ [Ti−1, Ti]. In [12] the convex functional G was defined on the domain D (defined above) by G(c) = N∑ j=1 Gλj (c). (2.12) From the uniqueness theorem [11, Theorem 3.5] we know that, under certain (computer-verifiable) conditions on the nature of the flows of certain associated vec- tor fields (which amount here to an admissibility restriction on the data u(K,T )), the condition that wλ = wc,λ for at least three distinct values of λ implies that c = (r, s) = (p, q). By [18, Proposition 42.6(1)] we know that if the convex func- tional G has a stationary point at (p, q) then it must have a global minimum there, and from the foregoing (assuming admissible data) that stationary point must uniquely occur at (p, q). So, the desired function pair (p, q) now appears as the unique global minimum of a convex functional with a unique stationary point. In practical numerics this is an important consideration, as many (if not most) least-square-type minimization methods suffer greatly from the minimization pro- cess getting stuck in spurious local minima. That this cannot happen here is one of the significant advantages of this approach. As it happens, in the FX option case the relevant functional G takes on quite tiny starting values (≈ 10−4 for the data in our example) and thus even though we have a well-posed minimization (so errors arising from ill-posedness are not an issue here [13]), this still constitutes a difficult inverse problem, because we are minimizing to zero and there is not much wiggle room. 166 I. KNOWLES, S. TAMANG EJDE/SI/02 Our solution in this paper is to pursue the volatility σ(K) directly, and rather than first minimizing the two-variable functional G(r, s) to obtain (p, q) followed by computing σ from (2.6), we now construct from G a new convex functional H(ρ) with a unique global minimum at ρ = σ. The philosophy here is that a one variable minimization is expected to be better behaved numerically than a two variable minimization. In the sequel we show that this is indeed the case for the function H(ρ) defined below. 3. The functional H(ρ) For a generic volatility ρ(K) define pρ(K) = e −2(rD−rF ) ∫K a dτ τρ2(τ) , (3.1) qρ(K) = 2pρ(K) K2ρ2(K) . (3.2) Note that from (2.2) and (2.4) we have p = pσ and q = qσ. For λ > 0 we define the non-negative functional Hλ(ρ) by c = (pρ, qρ) and Hλ(ρ) = Gλ(c) = ∫ b a pρ.(w ′ λ − w′λ,c)2 + (λ+ rF )qρ.(wλ − wc,λ)2. (3.3) Theorem 3.1. (a) Hλ(ρ) ≥ 0 for all ρ, and Hλ(ρ) = 0 if and only if wλ = wc,λ, where c = (pρ, qρ). (b) The first Gâteaux derivative of Hλ is given by H ′λ(ρ)[h] = ∫ b a (w′2λ − w′2c,λ)h1 + [(λ+ rF )(w2 λ − w2 c,λ)− 2β(wλ − wc,λ)]h2, (3.4) where h1(K) = 4(rD − rF )e −2(rD−rF ) ∫K a dk kρ2(k) ∫ K a h(τ)dτ τρ3(τ) , and h2 = − 4h K2ρ3 e −2(rD−rF ) ∫K a dτ τρ2(τ) + 8 K2ρ2 e −2(rD−rF ) ∫K a dτ τρ2(τ) (rD − rF ) ∫ K a h τρ3 dτ. (c) The second Gâteaux derivative of Hλ is given by H ′′λ (ρ)[h, k] = 2(L−1c,λ(e(h1, h2)), e(k1, k2)), (3.5) where e(h1, h2) = −(h1w ′ λ,c) ′ + h2[(λ+ rF )wc,λ − β], and (·, ·) denotes the usual inner product in L2[a, b] and h1 and h2 are calculated above from h and ρ. EJDE-2023/SI/02 INVERSE VOLATILITY PROBLEM 167 Proof. (a) This follows from Theorem 2.1(a). (b) We have H ′λ(ρ)[h] = lim ε→0 Hλ(ρ+ εh)−Hλ(ρ) ε = lim ε→0 Gλ(pρ+εh, qρ+εh)−Gλ(pρ, qρ) ε = lim ε→0 Gλ(pρ + εh∗1, qρ + εh∗2)−Gλ(pρ, qρ) ε = G′λ(pρ, qρ)[h1, h2]. where we define h∗1 and h∗2 by pρ+εh = pρ + εh∗1 and qρ+εh = qρ + εh∗2 and note that h1 = lim ε→0 h∗1 = lim ε→0 1 ε [pρ+εh − pρ] = lim ε→0 1 ε [ e −2(rD−rF ) ∫K a dτ τ(ρ(τ)+εh(τ)2 − e−2(rD−rF ) ∫K a dτ τρ2(τ) ] = lim ε→0 1 ε e −2(rD−rF ) ∫K a dτ τρ2(τ) { e −2(rD−rF ) ∫K a 1 τ [ 1 (ρ+εh)2 − 1 ρ2 ] dτ − 1 } = lim ε→0 1 ε e −2(rD−rF ) ∫K a dτ τρ2(τ) { e 2(rD−rF )ε ∫K a 1 τ 2ρh+εh2 ρ2(ρ+εh)2 dτ − 1 } = lim ε→0 e −2(rD−rF ) ∫K a dτ τρ2(τ) 2(rD − rF ) ∫ K a 1 τ 2ρh+ εh2 ρ2(ρ+ εh)2 dτ × { e2(rD−rF )ε ∫K a 1 τ 2ρh+εh2 ρ2(ρ+ε)2 dτ − 1 2(rD − rF )ε ∫K a 1 τ 2ρh+εh2 ρ2(ρ+εh)2 dτ } = 4e −2(rD−rF ) ∫K a dτ τρ2(τ) (rD − rF ) ∫ K a h τρ3 dτ, using that lim z→0 ez − 1 z = 1; and in similar fashion h2 = lim ε→0 h∗2 = lim ε→0 1 ε [qρ+εh − qρ] = lim ε→0 2 εK2 [ pρ+εh { 1 (ρ+ εh)2 − 1 ρ2 } + 1 ρ2 {pρ+εh − pρ} ] = − 4h K2ρ3 pρ + 2 K2ρ2 lim ε→0 1 ε [pρ+ε − pρ] = − 4h K2ρ3 e −2(rD−rF ) ∫K a dτ τρ2(τ) + 8 K2ρ2 e −2(rD−rF ) ∫K a dτ τρ2(τ) (rD − rF ) ∫ K a h τρ3 dτ. (c) This follows from Theorem 2.1(d). � We next calculate the L2[a, b]-gradient of the functional Hλ. 168 I. KNOWLES, S. TAMANG EJDE/SI/02 Theorem 3.2. Further to (3.4) we have H ′λ(ρ)[h] = ∫ b a [ 1 Kρ3(K) ∫ b K m(τ) dτ − n(K) ] h(K) dK, (3.6) where m(K) = 4(rD − rF )e −2(rD−rF ) ∫K a dτ τρ2(τ) × [ [(λ+ rF )(w2 λ − w2 c,λ)− 2β(wλ − wc,λ)] 2 K2ρ2 + (w′2λ − w′2c,λ) ] , (3.7) and n(K) = [(λ+ rF )(w2 λ − w2 c,λ)− 2β(wλ − wc,λ)] 4 K2ρ3 e −2(rD−rF ) ∫K a dτ τρ2(τ) , (3.8) and finally, g(K) = 1 Kρ3(K) ∫ b K m(τ) dτ − n(K) is the L2[a, b]-gradient of Hλ(ρ). Proof. From Theorem 3.1(b) we have H ′λ(ρ)[h] = ∫ b a { (w′2λ − w′2c,λ)4(rD − rF )e −2(rD−rF ) ∫K a dτ τρ2(τ) ∫ K a h(τ)dτ τρ3(τ) − h(K)[(λ+ rF )(w2 λ − w2 c,λ)− 2β(wλ − wc,λ)] 4 K2ρ3 e −2(rD−rF ) ∫K a dτ τρ2(τ) + [(λ+ rF )(w2 λ − w2 c,λ)− 2β(wλ − wc,λ)] 8(rD − rF ) K2ρ2 e −2(rD−rF ) ∫K a dτ τρ2(τ) × ∫ K a h τρ3 dτ, } dK = ∫ b a { 4(rD − rF )e −2(rD−rF ) ∫K a dk kρ2(k) ∫ K a h(τ)dτ τρ3(τ) × [ [(λ+ rF )(w2 λ − w2 c,λ)− 2β(wλ − wc,λ)] 2 K2ρ2 + (w′2λ − w′2c,λ) ] − h[(λ+ rF )(w2 λ − w2 c,λ)− 2β(wλ − wc,λ)] 4 K2ρ3 e −2(rD−rF ) ∫K a dτ τρ2(τ) } dK. Collecting the terms containing ∫K a h(τ)dτ τρ3(τ) , the above expression becomes∫ b a { m(K) ∫ K a h(τ)dτ τρ3(τ) − h(K)n(K) } dK, (3.9) where m(K) = 4(rD − rF )e −2(rD−rF ) ∫K a dk kρ2(k) × [ [(λ+ rF )(w2 λ − w2 c,λ)− 2β(wλ − wc,λ)] 2 K2ρ2 + (w′2λ − w′2c,λ) ] , (3.10) n(K) = [(λ+ rF )(w2 λ − w2 c,λ)− 2β(wλ − wc,λ)] 4 K2ρ3 e −2(rD−rF ) ∫K a dτ τρ2(τ) . (3.11) EJDE-2023/SI/02 INVERSE VOLATILITY PROBLEM 169 We now integrate the first term in the integrand in (3.9) by parts to isolate the h factor. Specifically,∫ b a m(K) ∫ K a h(τ)dτ τρ3(τ) dK = ∫ b a h(K) 1 Kρ3(K) [ ∫ b K m(τ) dτ ] dK, and hence from (3.9) H ′λ(ρ)[h] = ∫ b a h(K) { 1 Kρ3(K) ∫ b K m(τ) dτ − n(K) } dK, (3.12) which completes the proof. � Finally, analogously to (2.12), we set DH = {ρ : (pρ, qρ) ∈ D}, and define the convex functional H(ρ) on the domain DH by H(ρ) = N∑ j=1 Hλj (ρ). (3.13) 4. A recovery algorithm H(ρ) is a nonnegative convex functional since it is the sum of nonnegative convex functionals, and it also has a unique stationary point at σ. The idea here is that by using H rather than just one of the Hλ, in addition to gaining favourable uniqueness properties, we are blending additional time-based data into the inverse problem, and this is intended to improve the well-posedness of the problem, given that ill- posedness is basically due to a paucity of information. We note in passing from [13] that this inverse recovery is conditionally well-posed in the weak-L2 sense, so the recoveries are expected to be, and indeed were, quite stable. We minimize this functional for N = 20 using the steepest descent method to recover the coefficient σ(T,K). The L2-direction of steepest ascent for H at ρ is ∇L2H(ρ) = N∑ i=1 { 1 Kρ3(K) ∫ b K mi(τ) dτ − ni(K) } , where, for 1 ≤ i ≤ N , mi(K) = 4(rD − rF )e −2(rD−rF ) ∫K a dτ τρ2(τ) [[ (λi + rF )(w2 λi − w 2 c,λi) − 2β(wλi − wc,λi) ] 2 K2ρ(K)2 + (w′2λi − w ′2 c,λi) ] , (4.1) ni(K) = [(λi + rF )(w2 λi − w 2 c,λi)− 2β(wλi − wc,λi)] × 4 K2ρ(K)3 e −2(rD−rF ) ∫K a dτ τρ(K)2(τ) , (4.2) Instead of using this L2-gradient directly, it is numerically preferable to use the corresponding Neuberger gradient (see [17]) given that the L2-gradient has numerical problems that are extensively discussed in [15]. In particular, the L2- gradient is zero on the boundary of [a,b] given that wλ and wc,λ are equal there, and thus the algorithm is unable to properly recover σ. The Neuberger gradient smooths the L2-gradient and preserves boundary data during the descent, an important property not shared by other descent techniques. Our Neuberger gradient g = 170 I. KNOWLES, S. TAMANG EJDE/SI/02 ∇H1H can be found from an L2-gradient ∇L2H by solving the Dirichlet-Neumann boundary value problem −g′′ + g = ∇L2H, g(a) = g′(b) = 0, (4.3) where [a, b] represents the generic K-interval. The choice of the Neumann condition at b allows for some flexibility in the shape of the iterates. Forcing g(b) = 0 appears to be too confining and inhibits the descent. We have option prices for discrete sets of strikes and maturities. We generated the function v(K,T ) by linearly interpolating the option price in both strike and maturity. The function v(K,T ) was mollified (c.f. [14, §6]) so that it could be dif- ferentiated, and the derivative vK(K,T ) was found using central differences. For 20 fixed values of λ the functions v(K,T ) and v′K(K,T ) were Laplace transformed us- ing (2.1) to wλ(K) and w′λ(K) respectively. The functions pσ0 (K) and qσ0 (K) were initialized using (2.2) and (2.4) with the initial σ chosen to be the implied volatil- ity. We performed a series of descents in σ using the aforementioned Neuberger steepest descent algorithm such that the functional could not be minimized any further. The line minimization of f(α) = H(σ0−α∇H1H(σ0)) in α was done using the well known Brent minimization technique, by adapting the one-variable code in the Numerical Recipes in C function brent(). To avoid possible catastrophic cancellation in the Simpson rule formula used in the calculation of the integrals in the formula (2.8) for the functional Hλ, we used the alternate formula (3.3) instead. Below is the steepest descent algorithm used to get one descent step in σ: (i) Initialize σ = σ0, for example using a known implied volatility. (ii) Set c0 = (pσ0 , qσ0). (iii) For all i, find wλi,c0 and w′λi,c0 by solving (2.7) with λ = λi. (iv) Find the L2 gradient of H at σ0, ∇L2H(σ0). (v) Find the Neuberger gradient, ∇H1H(σ0). (vi) Find α = αmin that gives the lowest value of f(α) = H(σ0 − α∇H1H(σ0)). (vii) Compute σnew = σ0 − αmin∇H1H(σ0). One may then iterate by overwriting σ0 with σnew and repeating steps (ii) through (vii) until H fails to descend using the descent direction −∇H1H(σnew). 5. Results The most widely traded of the FX options is the EURUSD option. To recover the coefficient function σ(T,K) we wrote a C program, which used, as the initial guess in the descent process, the “implied volatility” obtained directly from the standard option price formula of Garman and Kohlhagen [8] by substituting the known option price and solving for the implied volatility σ as an unknown. The implied and recovered volatility surfaces and the convergence graph are shown in Figures 1–3 6. Conclusions We have presented herein a new inverse algorithm for the computation of foreign exchange volatility from foreign exchange option data. The major difficulty that we encountered in this particular inverse problem was that the typical values for our previously effective ([12]) minimizing functional G(r, s) (see (2.12)) were quite tiny (of the order 10−4) and as we were minimizing to zero the two-variable functional EJDE-2023/SI/02 INVERSE VOLATILITY PROBLEM 171 Figure 1. Implied volatility surface for EURUSD options Figure 2. Recovered EURUSD option volatility surface σ(K,T ) G was too insensitive to handle this task. So we reformulated to a new one-variable functional H(ρ) defined by (3.13) and this proved quite effective. It seems likely that this new approach would also improve the volatility recovery not only in [12] but elsewhere, such as in groundwater flow-parameter recovery, where the smallness of the functional values is not an issue. References [1] M. Avellaneda, C. Friedman, L. Holmes, L. Sampieri; Calibrating volatility surfaces via relative entropy minimization. Appl. Math. Finance, 4:37–64, 1997. [2] Fischer Black, Myron Scholes; The pricing of options and corporate liabilities. Journal of Political Economy, 81:637–654, 1973. 172 I. KNOWLES, S. TAMANG EJDE/SI/02 Figure 3. Minimization of the functional H [3] J. N. Bodurtha, M. Jermakyan; Non-parameric estimation of an implied volatility surface. J. Computational Finance, 2:29–61, 1999. [4] Ilia Bouchouev, Victor Isakov; The inverse problem of option pricing. Inverse Problems, 13:L11–L17, 1999. [5] Ilia Bouchouev, Victor Isakov; Uniqueness, stability and numerical methods for the inverse problem that arises in financial markets. Inverse Problems, 15:R95–R116, 1999. [6] Ilia Bouchouev, Victor Isakov; Recovery of volatility coefficient by linearization. Quant. Fi- nance, 2:257–263, 2002. [7] B. Dupire; Pricing with a smile. RISK, 7:18–20, 1994. [8] Mark B. Garman, Steven W. Kohlhagen; Foreign currency option values. Journal of Inter- national Money and Finance, 2(3):231 – 237, 1983. [9] Philip Hartman; Ordinary differential equations. S. M. Hartman, Baltimore, Md., 1973. Cor- rected reprint. [10] Victor Isakov; The inverse problem of option pricing. Preprint, 2004. [11] Ian Knowles; Uniqueness for an elliptic inverse problem. SIAM J. Appl. Math., 59(4):1356– 1370, 1999. [12] Ian Knowles, Li Feng, Ajay Mahato; The inverse volatility problem for European options. In Proceedings of the Variational and Topological Methods: Theory, Applications, Numerical Simulations, and Open Problems, volume 21 of Electron. J. Differ. Equ. Conf., pages 183– 195. Texas State Univ., San Marcos, TX, 2014. [13] Ian Knowles, Mary A. LaRussa; Conditional well-posedness for an elliptic inverse problem. SIAM J. Appl. Math., 71:952–971, 2011. [14] Ian Knowles, Tuan A. Le, Aimin Yan; On the recovery of multiple flow parameters from transient head data. J. Comp. Appl. Math., 169:1–15, 2004. [15] Ian Knowles, Robert Wallace; A variational method for numerical differentiation. Numerische Mathematik, 70:91–110, 1995. [16] R. Lagnado, S. Osher; A technique for calibrating derivation of the security pricing models: numerical solution of the inverse problem. J. Computational Finance, 1:13–25, 1997. [17] J. W. Neuberger; Sobolev gradients and differential equations, volume 1670 of Lecture Notes in Mathematics. Springer-Verlag, Berlin, second edition, 2010. [18] Eberhard Zeidler; Nonlinear functional analysis and its applications. III. Springer-Verlag, New York, 1985. Variational methods and optimization, Translated from the German by Leo F. Boron. EJDE-2023/SI/02 INVERSE VOLATILITY PROBLEM 173 Ian Knowles Department of Mathematics, University of Alabama at Birmingham, Birmingham AL 35294, USA Email address: iknowles@uab.edu Sundar Tamang Department of Math and Computer Science, Western New Mexico University, Silver City, NM 88062, USA Email address: sundar.tamang@wnmu.edu 1. Introduction 2. Reconstruction of FX volatility 2.1. Properties of the functional G 3. The functional H() 4. A recovery algorithm 5. Results 6. Conclusions References