Recovery of time dependent volatility coefficient by linearization

We study the problem of reconstruction of special special time dependent local volatility from market prices of options with different strikes at two expiration times. For a general diffusion process we apply the linearization technique and we conclude that the option price can be obtained as the sum of the Black-Scholes formula and of an operator $W$ which is linear in perturbation of volatility. We further simplify the linearized inverse problem and obtain unique solvability result in basic functional spaces. By using the Laplace transform in time we simplify the kernels of integral operators for $W$ and we obtain uniqueness and stability results for for volatility under natural condition of smallness of the spacial interval where one prescribes the (market) data. We propose a numerical algorithm based on our analysis of the linearized problem.


Introduction. Basic results
The Black-Scholes formula [8] is an efficient method to price financial derivatives under the assumption that the stock price is log-normally distributed. However, prices of options with different strikes generated by the Black-Scholes formula differ from observed market prices [14]. One way to reconcile the differences is to replace the constant volatility by a more volatility that might depnd on time and the underlying asset. This approach is very popular in applications.
It is usually difficult to achieve a unique and stable fitting of the model to actual market prices. While in the constant volatility can be found very efficiently and uniquely from one option price, the space and time dependent volatility function must be restored from collection of simulteneous option quotes with different strikes and expiration times. Even though many numerical algorithms for this inverse problem have already been published [1], [4], [19], their convergence properties have not been satisfactory for practitioners. In particular, since these algorithms use minimization of a non convex regularized misfit functional, there are well known difficulties with avoiding local minima, and their convergence is not guaranteed. We refer to [5] for a survey of available theory and of numerical methods. In the time dependent case a relation between implied and local volatilities was discovered [2], [3] and a convergent numerical algorithm was designed. In many practical situations the interval ω * with strike prices is relatively small and values of the volatility coefficient outside ω * do not influence option prices inside ω * very much due to very fast diffusion. This suggest a possibility of linearization. One of these linearizations was suggested in [7] where linearized inverse problems are still complicated and no numerical algorithms were discussed. In addition, expiration times are sparse, so it is hard to expect a detailed recovery of general time dependence of volatility. In [6] we assumed that volatility coefficient only depends on stock price, considered a linearized version of the inverse problem, obtained a simple convenient representation of this linearization, and desinged and tested a very fast and reliable numerical algorithm based on this linearization.
In this paper we attempt to recover volatility which linearly depends on time while the additional data are given for two expiration times. So we are making use of realistic data to predict volatility in near future which seems to be most important for market decisions. Below we outline the results in more detail.
For any stock price, 0 < s < ∞, and time , 0 < t < T , the price u for an option expiring at time T satisfies the following partial differential equation Here, σ(s, t) is the volatility coefficient that satisfies 0 < m < σ(s, t) < M < ∞ and is assumed to belong to the Hölder space C λ (ω * ), 0 < λ < 1 for some interval ω * ⊂ (0, +∞) and outside this interval, and µ and r are, respectively, the risk-neutral drift and the risk-free interest rate assumed to be constants. The backward in time parabolic equation (1.1) is augmented by the final condition specified by the pay off of the call option with the strike price K u(s, T ) = (s − K) It is known (e.g. see [5] for Hölder σ) that there is a unique solution u to (1.1),(1.2) which belongs to C 2 ((ω * )×(0, T ] and to C((0, ∞)×[0, T ]) and satisfies the bound |u(s, t)| < C(s+1). The inverse problem of option pricing seeks for σ given Here s * is market price of the stock at time t * , u * (K) denote market price of options with different strikes K for a given expiry T . We suppose that ω * contains s * and is small, and attempt to recover volatility in the same interval. In this paper we assume that I * = {T 1 , T 2 } for some 0 < T 1 < T 2 we are trying to recover the volatility coefficient where f * 0 , f * 1 are small C(ω * )-perturbation of constant σ 2 0 and f * 0 = f * 1 = 0 outside ω * . Constant σ 0 is the implied volatility defined as the (unique) solution to the Black-Scholes equation where the normal distribution function Then the option price can be calculated (up to a quadraticaly small error) as the sum of the Black-Scholes formula with volatility σ 0 and of a certain linear operator U ((f * 0 , f * 1 )). To obtain our results we will use that the option premium u(., .; K, T ) satisfies the equation dual to the Black-Scholes equation (1.1) with respect to the strike price K and expiry time T : The equation (1.6) was found by Dupire [10] and rigorously justified, for example, in [5].
In the section 2 we use the standard linearization procedure and derive a partial differential equation for the principal linear (with respect to (f * 0 , f * 1 ) part V of u in new logarithmic variables. After another substitution this equation is reduced to the heat equation with the right-hand side linear with respect to (f * 0 , f * 1 ). We will drop * in new variables. Due to difficulties with the linearized inverse problem in section 3 we propose its simplified version. This simplification is a good approximation in the most practically important case of small ω * . For the simplified inverse problem by using the Fourier transformation with respect to s we obtain uniqueness, stability, and existence results in Sobolev spaces when the data are given for all strikes at two expire times. Due to strong diffusion the data (1.3) are small outside ω * , so the original data can be extended onto (0, +∞) \ ω * as some C 2 -small function.
In section 4 we study the linearized inverse option pricing problem. We apply the Laplace transform to explicitly evaluate the time integrals in the standard integral representation of the solution of the direct option pricing problem. As a result we arrive in the following system of linear integral equations for (f 0 , f 1 ) (obtained from (f * 0 , f * 1 ) via a simple transformation): with the kernels given by the error function and the right-hand side F j (x) computed from the market data (1.3). Using this explicit formulae we prove that the system of integral equations (1.7) has not more than one solution under a certain constraint on ω which is determined by the range of strikes for which the data (1.3) are given. It follows also a that for small change in F j that can be viewed as a market bid-ask spread or any other noise in the market data (1.3), solutions changes remain small as well. In other words, the linearized version of the volatility reconstruction is stable. Numerical results obtained by other authors and by author of this paper using different methodologies ( [1], [4], [6], [11], [19]) exhibit similar reconstruction properties. In ( [1], [19]) they are looking for the volatility σ = σ(t, s) from the data given for several T . We believe that values of the volatility for t * < t < T ( especially in the near future) are of most interest. But in the time dependent case the available data are not sufficient to identify it. That is why we think one should assume special time dependence and utilize other analytic features of this inverse problem. In [6] we considered time independent σ and use the data only at T 1 .
Our current contribution is twofold. First, we demonstrated uniqueness, stability, and existence of solution for the approximate linearized inverse problem. Furthemore, we are able to explicitly describe (Theorem 4.1) the domain where the stable reconstruction is ensured for the linearized inverse problem. The closest result in this direction is obtained in [5], where the linearization was replaced by an approximate linearization, and in [6] where only f = 0 and I * = {T 1 }. As in [6], the answer is obtained in terms of the (more complicated) ratios of the spatial variable to the square root of time to maturity. It is very much consistent with trader's intuition.
Our approach suggests the actual reconstruction methodology. Algorithms proposed in [1], [7], [19] attempt to use standard methods (regularized best L 2 -fit) to numerical solution of inverse problems. They involve very heavy computations with poor convergence properties which is typical for inverse parabolic problems. The numerical algorithm in [4] is based on some features of the inverse option pricing problem. In fact, it is a second order correction of the Black-Scholes formula. Unfortunately, that this correction lacks a rigorous justification and numerically it is not very efficient. In ( [7]) they used perturbation methods which look reasonable when volatility is close to a constant, but the corresponding linear integral equation remains complicated as well as its proposed numerical inversion. In [6] we were able to transform the linearized inverse problem into a one-dimensional integral equation with an explicit kernel. This led to a very fast, simple, and stable numerical algorithm for the reconstruction of the local volatility. In this paper we transform the linearized inverse option pricing problem into the system of linear integral equations (1.7). We already started numerical solution of this system.

Linearization at constant volatility
The substitution transforms the equation (1.6) and the initial data (1.2) into while the additional (market) data become Here ω is the transformed interval ω * (ω * in y-variables (2.1)) and U j (y) = u * (s * e y , T j ).
Observe that τ j = T − t * j . The equations (2.2) and (2.3) for functions U (τ, y), a(y) form the so-called inverse parabolic problem with the final overdetermination. The known uniqueness conditions for this problem [15], section 6.2, [17], section 9.2, are not satisfied in our particular situation, and we are not aware of any uniqueness result.
To derive the linearized inverse problem we observe that due to the assumption (1.4), Here V 0 solves (2.2) with a = σ 0 and v is quadratically small with respect to (f 0 , f 1 ), while the principal linear term V satisfies the equations where V j (y) = U j (y) − V 0 (y, τ j ). One can completely justify this linearization by using standard theory of direct parabolic boundary value problems [13], [18], as it was done in [5], pp. R 103-104, for the inverse option pricing problem or in [17], section 4.5, for some elliptic inverse problems.
The new substitution V = e cy+dτ W (2.7) simplifies (2.5) to with the additional final data Due to analytic difficulties, using some features of our inverse problem we will further simplify (2.5) to the approximate linearized inverse problem by replacing α with the best approximating y independent function α 1 : In particular for computational purposes we replace R by a finite interval Ω containingω: with the additional final data As in [5], Theorem 4, one can show that uniqueness and stability for the approximate inverse problem imply uniqueness and stability for the linearized inverse problem when the interval ω is small.

The approximate linearized inverse problem
In this section we prove the following unique solvability result Then there is a unique solution to the approximate linearized inverse problem (2.10), (2.11).
Moreover there is a constant C depending only on τ 1 , τ 2 , σ 0 such that a solution (f 0 , f 1 ) to the approximate linearized inverse problem satisfies the bound We use the Fourier transform of a function f (y, τ ) with respect to y and the standard norm of a function f in the Sobolev space H (l )(R).

Proof:
Applying the Fourier transform to the initial value problem (2.10) we yield where we let Solving the initial value problem for the ordinary differential equation we obtain Hence the inverse problem (2.10),(2.11) is equivalent to the system of two integral equations Solving this linear algebraic system we obtain Integrating by parts, By elementary substitution θ = τ j ρ, Using the Taylor expansion for the exponential function and the definition of Ξ = 1 where |O(ξ)| ≤ C|ξ|, when |ξ| < 1. So calculating the integral of ρ Splitting the integration interval and integrating by parts again in the first integral, we will have when 1 ≤ |ξ|, if we repat the integration by parts. Hence Using that Ξ = 1 √ 2 σ 0 ξ, from (3.6) (3.7) we derive that Similarly to the derivation of (3.7) we obtain if |ξ| > 1, and therefore, Using the definition of Sobolev norms via the Fourier transform we obtain the bound (3.1). This bound implies uniqueness and stability. The formulae (3.4) give a solution explicitly for smooth compactly supported data (W 1 , W 2 ), combined with the density of such data in H (2) (R) × H (2) (R) and the stability estimate (3.1) it implies existence of solution in the needed Sobolev space.
The proof is complete. Moreover there is a constant C depending only on τ 1 , τ 2 , B, σ 0 such that a solution (f 0 , f 1 ) to the approximate linearized inverse problem satisfies the bound The proof is similar to the above proof of Theorem 3.1 if instead of the Fourier transform we use orthonormal (eigenfunction) series It suffices to replace ξ by n. By elementary means ( using sharp bounds in Taylor formula ) one can explicitly evaluate constants C. We did not do it to simplify the exposition.

Uniqueness and stability for the linearized inverse problem
The purpose of this section is to prove the following then a solution (f 0 , f 1 ) ∈ L ∞ (ω) × L ∞ (ω) to the linearized inverse option pricing problem (2.5), (2.6) is unique. Moreover, there is a constant C depending only on σ 0 , τ 1 , τ 2 , b, such that We remind that f ∞ (ω) is essential supremum of |f | over ω. In particular, when f is continuous on ω is it just max|f (x)| over x ∈ ω.
From numerical experiments in [19] and in [6] we can see that values of σ outside ω are not essential, due to a very fast decay of the Gaussian kernel α in y.
Proof: When f 1 = 0 this result was obtained in [6]. So by using the linearity with respect to (f 0 , f 1 ) we can assume that f 0 = 0.
The well-known representation [13] of the solution to the Cauchy problem (2.8) for the heat equation yields As in [6] we will simplify K(x, y; τ ) by using the Laplace transform Φ(p) = L(φ)(p) of φ(τ ) with respect to τ . Since the Laplace transform of the convolution is the product of Laplace transforms of convoluted functions, we have where we used the formula for the Laplace transform of τ − 1 2 e − β τ [20], p.526, formula 16, and that L(τ φ(τ )) = − d dp L(φ(τ )), [20], p.499. Applying the formula for the inverse Laplace transform of the function 1 p e −γ √ p [20], p. 528, formula 42, and using the relation between Laplace transforms of a function and its integral [20], p. 499, and [20], p. 526, formula 16, we conclude that where we switched order of integration with respect to ρ and θ and integrated by parts. To find the inverse Laplace transform for (4.4) we will again use the known above relation between Laplace transform of a functions and of its integral, to find where we utilized the substitution θ = γ 2 4ρ 2 and integrated by parts. By using (4.5), (4.6) with γ = √ 2 |x−y|+|y| σ 0 we yield Now the formula of Lemma 4.2 for K 1 follows from (4.4).

Lemma 4.3
The linearized inverse problem (2.8), (2.9) implies the following Fredholm system of the linear integral equations where and

Proof:
The formula for A j1 f 0 was derived in [6]. Now we will consider A j2 f 1 .
Differentiating once more and using that ∂ x sign(x − y) = 2δ(x − y) we yield ∂ 2 x ω The proof is complete. In our opinion, several cancellations in the the proof of Lemma 4.3 suggest that our approach is quite suitable for this problem Proof of Theorem 4.1: Due to Lemma 4.3 to prove Theorem 4.1 it suffices to show uniqueness and stability of solution (f 0 , f 1 ) of (4.7). By linear algebra (4.7) is equivalent to the system Using (4.8), where sup is over x ∈ ω = (−b, b) (or , equivalently, over x ∈ (0, b) ). Similarly, To explain the last inequality we observe that where we used the substitution θ = √ X 2 + σ 2 and the well known probability integral. Using that |x| ≤ |x − y| + |y| and hence From (4.9), (4.10), (4.11) by the triangle inequality we yield Adding the first inequality and the second inequality multiplied by √ τ 1 τ 2 we obtain Using the condition (4.1) we conclude that the factors of f 0 , f 1 on the left side are positive, so (by dividing by smallest of them both sides) we arrive at (4.2).
At this point we do not have an explicit bound on C in (4.2). We can not claim the existence of the solution of the linearized inverse option pricing problem, but we can show that under the conditions of Theorem 4.1 the range of the operator defined by the left side of (4.12) (considered from C(ω) × C(ω) into C 2 (ω) × C 2 (ω)) has the codimension 4. Indeed, for any (W 1 , W 2 ) ∈ C 2 (ω) × C 2 (ω) the Fredholm system (4.7) (which follows from (4.12)) has a unique solution (f 0 , f 1 ) ∈ C(ω) × C(ω), due to uniqueness guaranteed by theorem 4.1. Obviously, using this (f 0 , f 1 ) in the left side of (4.12) produces the right side C 1 + C 2 x with some constant vectors C 1 , C 2 . Observe that we can not conclude that C 1 = 0, C 2 = 0, so while the system (4.7) follows from the system (4.12), it is not equivalent to this equation. In any event, the range of the operator from left side of (4.12) has the codimension not greater than 4 We will show that it is exactly 4.
To do it we will observe that the system (4.12) has no solution when W 1 , W 2 are linear functions of x. Indeed, let (f 0 , f 1 ) be this solution. In the both cases W ′′ j (x) = 0, so the solution (f 0 , f 1 ) to the equation (4.12) is zero due to uniqueness guaranteed by Theorem 4.1. But then W j = 0 and we have a contradiction.
At present we do not know an exact description of the range of the (matrix) linear integral operator defined by the left side of (4.12).
For analytical and especially numerical purposes it can be useful to handle f j ∈ C(R). The results of section 4 can be extended to these functions which are constants on the intervals (−∞, −b], [b, +∞). The analysis will be a little more technical and the condition (4.1) will be replaced by more complicated one.

Conclusion
Uniqueness result of this paper can be certainly extended to the many-dimensional case and probably to general parabolic equations with variable coefficients. The inverse option pricing problem is a particular case of the more general inverse diffusion problem which has a probabilistic interpretation. We are not aware of any uniqueness results about recovery of space dependent diffusion rate from probability of distribution at a fixed moment of time, and moreover of a (special) time dependent rate from probability at few moments of time. The method of this paper can be applied at least to a linearized version of this inverse probabilistic problem. Observe that while in case of one unknown space dependent coefficient (scalar case) one can obtain complete uniqueness results [16] under certain monotonicity conditions on boundary data (which are not satisfied for the inverse option pricing problem) by using maximum or positivity principles, the (vectorial) case of several unknown spacial functions corresponds to systems of equations when maximum principles are only in very special situations. There are almost no uniqueness results in the vectorial case.
The proposed reconstruction algorithm is expected to perform very well when volatility is not changing fast with respect to stock price s and is changing very slow with respect to time. Sudden and dramatic changes of market situations most likely can not be properly described by our model and more generally by the Black-Scholes equation. Probably, a minor modification of the proposed model (replacing R in (2.2) by a finite interval) can eliminate difficulties with existence theorem and generate even better numerical algorithms. Observe, that to find continuous f j the data U j must be at least twice differentiable on ω, so the real market data are in need of a proper interpolation, minimizing the size of second derivatives of U j . A choice of an appropriate smoothing interpolation and an intensive numerical testing will be a subject of future work. Our immediate goal was to develop some mathematical theory of the inverse option pricing problem and to suggest a numerical algorithm. The next step is an additional testing on simulated and real data. We already made preliminary numerical tests which show very fast convergence of iterative methods. The proposed method provides with simplest representation of the inverse problem and hence promises to be fastest and most stable of existing numerical methods.
For simplicity we considered only European options. We hope to adjust the linearization technique to American and more complicated options, which are in particular described by free boundary problems and multidimensional parabolic equations. So far there are actually no analytic results in this very important practical case.
Aknowledgement: This research was in part supported by the NSF grant DMS 10-08902 and by Emylou Keith and Betty Dutcher Distinguished Professorship at the Wichita State University.