Recovery of local volatility for financial assets with mean-reverting price processes

This article is concerned with the model calibration for financial assets with mean-reverting price processes, which is an important topic in mathematical finance. The discussion focuses on the recovery of local volatility from market data for Schwartz(1997) model. It is formulated as an inverse parabolic problem, and the necessary condition for determining the local volatility is derived under the optimal control framework. An iterative algorithm is provided to solve the optimality system and a synthetic numerical example is provided to illustrate the effectiveness.


1.
Introduction. Stochastic models for the evolution of the financial assets' prices are at the core of the mathematical finance. Pricing of a derivative product begins with a reasonable modelling of its underlying asset price process.
In the celebrated Black-Scholes model (Black & Scholes (1973) [3]), the price of a stock S, is assumed to follow the Geometric Brownian motion dS t = µS t dt + σS t dW t , (1.1) where µ is the drift, σ is the volatility and W t is a standard Brownian motion. For other financial asset such as a physical commodity, its spot price exhibits mean-reverting nature, because of the supply and demand fluctuations and the fixed costs of adjusting the long-term supply of a commodity or the fixed costs of shifting consumption from one commodity to another. In this case, a Geometric Brownian motion as in the Black-Scholes model cannot properly capture its price behavior, a mean-reverting stochastic process is expected.
A simple mean-reverting model for an asset price is represented by the following Ornstein-Uhlenbeck process (Uhlenbeck & Ornstein (1930) [27]): where α is the mean-reverting rate, µ is the long-term mean-reverting level, W t is a standard Brownian motion and σ is the volatility. In this model, the price mean reverts to the long-term level µ at a speed given by the mean-reverting rate α, which is taken to be strictly positive. If the spot price is above the long-term level µ, then the drift of the spot price will be negative and the price will tend to revert back towards the long-term level. Similarly, if the spot price is below the long-term level then the drift will be positive and the price will tend to move back towards µ. Note that at any point in time the spot price will not necessarily move back towards the long-term level as the random change in the spot price may be of the opposite sign and greater in magnitude than the drift component.
The stochastic models typically incorporate some parameters (e.g. volatility, mean-reverting rate and long-term mean-reverting level, etc, in the above (1.1) and (1.2)) that may be constants or deterministic functions. The successful application of the model in pricing, hedging and other trading activities will critically depend on how the parameters are specified.
Usually, these parameters are not directly observable in the market. One needs to estimate the parameters by calibrating the model-produced prices to the observed market prices of actively traded derivatives, such as vanilla European options. This "mark-to-market" model calibration serves the purpose to make the pricing and hedging of the over-the-counter exotic derivative assets be within a consistent framework with the prices of more liquid exchange-traded derivative assets.
The recovery of local volatility function from market data is one very important example of model calibration. In this article, we shall consider this problem for the Schwartz(1997) mean-reverting underlying process ( [23]).
The remaining sections are organized as follows: In section 2, we introduce the socalled local volatility model. In section 3, we first set up the recovery problem, and then reformulate it into a more standard inverse parabolic problem with terminal observation. In section 4, we derive the necessary condition that the solution of the inverse problem should satisfy under an optimal control framework. Computational issues on how to solve the optimality system are described in section 5, a synthetic numerical example is provided to illustrate the effectiveness and accuracy of the proposed algorithm as well.
2. Local volatility model. A quantity of fundamental importance to the pricing of a financial derivative asset is the stochastic component in the evolution process of its underlying price, the so-called volatility, which is a measure of the amount of fluctuation in the asset prices, i.e., a measure of the randomness. Obtaining estimates for the volatility is a major challenge in market finance.
As is well known, much evidence suggests that the constant volatility model is not adequate (cf. Rubinstein (1994) [22]; Dupire (1994) [12]; Derman & Kani (1994) [11] and the references therein). Indeed, numerically inverting the Black-Scholes formula on real market data supports the notion of asymmetry with stock price (volatility skew), as well as dependence on the time to expiration (volatility term structure). The challenge is to accurately (and efficiently) model these effects.
A few different approaches have been proposed for modelling the volatility effects seen in the market. Merton (1976) [21] considered a jump-diffusion process for the underlying asset. Dupire (1994) [12] and Derman & Kani (1994) [11] made a direct extension to the Black-Scholes model -instead of assuming volatility to be a constant σ, they assumed that the volatility is a deterministic function of asset price and time, σ(S, t). Their approach is called "local volatility (LV) model". Also, there is a class of "stochastic volatility (SV) model", starting with Stein & Stein (1991) [25] and Heston (1993) [15], where volatility itself follows a stochastic process. Recently, some "local stochastic volatility (LSV) models" were introduced in the literature, which have become the industry market standard for FX and equity markets (cf. Lipton, Gal, & Lasis (2014) [20]).
Among these approaches, local volatility model is the easiest one to implement. Its advantages over the jump or stochastic model include that, no non-traded source of risk such as the jump or stochastic volatility is introduced (Dupire (1994) [12]), so that the "completeness" of the model, i.e., the ability to hedge the derivative using its underlying asset, is maintained. Completeness is ultimately important since it allows for arbitrage pricing and hedging (Dupire (1994) [12]).
In recent years, the research on local volatility model has attracted a good deal of attention in both academic and practical area. Various approaches have been proposed on how to recover the local volatility function from the most liquid market data, usually the data used is the market European option price on the asset (cf. Avellaneda, Friedman, Holmes & Samperi (1997) [1]; Lagnado& Osher (1997) [18]; Bouchouev & Isakov (1997 [14] and the references therein).
However, all the available research has been done based on a log-normal underlying process, which is suitable for equity prices, but not for the mean reverting asset prices considered here.
In this article, we shall suggest a heuristic approach to recover the local volatility function from the market data for asset prices which follow a mean-reverting process.
3. Formulation as an inverse parabolic problem. Assume that under a riskneutral measure the underlying asset follows a mean-reverting process (Schwartz (1997) [23]) dS = α(κ − ln S)Sdt + σSdW, (3.1) Denote by C(S, t; K, T ) the European call option price on the underlying asset S at time t, with strike price K and expiration T .
From the standard risk-neutral pricing results, under the dynamics (3.1), the discounted value e −rt C(S t , t; K, T ) is a martingale, where r is risk-free interest rate. By using the Itô-Doeblin formula (cf. Shreve (2004) [24]), the differential of this martingale is Setting the dt term in the above to be zero (due to the property of the martingale), we conclude that the option value C(S, t; K, T ) satisfies the following partial differential equation (cf. Clewlow & Strickland (2000) [7]): The problem of recovering the local volatility amounts to, in the continuous time setting, determining the coefficient σ, such that the solution of (3.2) fits the current market prices of European options at (s * , t * ) for different strikes K and/or maturities T . From the mathematical point of view, this is an inverse problem of partial differential equation (PDE), but it is not a standard one, since it requires determining the coefficient of the pricing equation from a series of observed values of the solution corresponding to various parameters K and/or T , at a fixed point (s * , t * ). In a standard inverse problem, the coefficient is determined from the observed values of the solution corresponding to various values of the variables (S, t), for fixed parameters K and T .
In the following, we make use of the well-known property of the fundamental solution and convert the original problem into a standard inverse parabolic problem with terminal observation, for the case where the local volatility σ is a function of S only.
Let G(S, t; K, T ) = C KK (S, t; K, T ), directly differentiate on both sides of (3.2), we have LG where δ(S − K) is the Dirac delta function concentrated at K. G(S, t; K, T ), as a function of (S, t), is called the fundamental solution of operator L. We can show that, G(S, t; K, T ), as a function of (K, T ), is the fundamental solution of its dual operator, i.e., G satisfies where H is the Heviside function. Now, recovering σ(S) in (3.2) from the market European option prices C is equivalent to recovering σ(K) in (3.5) such that U (S * , t * ; K, τ * ) = C K (S * , t * ; K, τ * ), where τ * = T − t * and C K (S * , t * ; K, τ * ) is the derivative of C with respect to K obtained from the market price of C. In theory, we need the price of C for all strikes to compute this derivative, which is not realistic. We shall use some interpolation on C across the available strikes and then get the derivative.
Perform the following change of variables:
Let us add the regularization term to the cost functional and consider the following optimization problem: where A = a ∈ C(R)|0 < a 0 ≤ a(y), a y ∈ L 2 (R) (here the condition 0 < a 0 ≤ a(y) is specified to ensure unique solvability of the Cauchy problem (3.6).) In J(a), rather than a 2 , a y 2 is chosen as the penalty function to get the "smoothest" minimizer, and β, as a weight coefficient, represents the trade-off between the accuracy and the smoothness of the minimizer, this follows the lines of Tikhonov's method.
The recovery of the local volatility from the market data is reduced to seeking for a minimizer of (4.1).
In the following, we shall derive the necessary optimal condition that the minimizer should satisfy.
Letā ∈ A be a minimizer. Since A is a convex set, for any h ∈ A and θ ∈ [0, 1], we haveā θ = (1 − θ)ā + θh ∈ A. Then, for any given h ∈ A, the function is well defined and reaches its minimum at θ = 0, so we must have where V θ (y, τ ) is the solution of (3.6) corresponding to a =ā θ . By direct differentiation with respect to θ on both sides of the equation for V θ , we get where ξ satisfies whereV (y, τ ) is the solution of (3.6) corresponding to a =ā. Letφ be the solution of the adjoint equation Moreover, by (4.4), the left hand side in the above is equal to
Remark. Actually, we can prove that any minimizerā is a weak solution to the following complementarity problem: Proof. Supposeā ∈ A solves (4.5) andā yy ∈ L 2 (R). Then,ā − a 0 ≥ 0 and, for any It follows from (4.5) that Fix any D ⊂⊂ R, let {χ n } be a sequence of functions from C ∞ c (R) such that 0 ≤ χ n ≤ 1 and χ n → χ D a.e. in R, where χ D is the characteristic function of D, i.e., This, combined with (4.9) and the arbitrariness of D gives, Thus,ā is a solution of complementarity problem (4.6).

5.
1. An iterative algorithm. From the last section, we know that the minimizer a(y), together with the corresponding stateV (y, τ ) and costateφ(y, τ ), satisfies the optimality condition consisting of a state equation, an adjoint equation and a complementarity problem: We shall use an iterative method to solve the above system and get the numerical solution ofā(y). The local volatility,σ(S), is then obtained byσ(S) = 2ā(ln S/S * ). The current underlying price is S * and the current time is t * . We are given the market option price with respect to different strikes K i 's for the maturity T . Using a smooth interpolation along the K i 's, we get the market price curve C * (K) and thus its derivative curve with respect to K, C * K (K). Then V * (y) is determined by its definition V * (y) = e r(T −t * ) C * K (S * e y ). The iterative procedure to get σ(K) is as follows: Step 1: Make an initial guess forā(y), say, a old (y). Substituteā = a old into equation (5.1) and solve forV (y, τ ) using a finite-difference method; Step 2: The difference betweenV (y, τ * ) and V * (y) gives the terminal condition of (5.2). With this condition and the guess forā = a old in step 1, we solve (5.2) forφ(y, τ ) using a finite-difference method; Step 3: WithV (y, τ ) andφ(y, τ ), obtain f (y;V ,φ) by numerical integration.
Substitute f (y;V ,φ) into the complementarity problem (5.3) and solve for a newā(y), call it a new (y); Step 4: Justify the convergence of the iteration. if |a new (y) − a old (y)| is smaller than some prescribed tolerance , the iteration is terminated. Otherwise, we use a new (y) as a new guess and return to Step 1 above for another round of iteration. In Step 3, the complementarity problem, after a finite-difference discretization, is solved by the progressive over-relaxation (PSOR) method. A good discussion on the progressive over-relaxation (PSOR) method can be found in (Cottle, Pang & Stone (1993) [9]). From the reference we see that our discretized problem satisfies some nice properties which ensures the suitability of applying the particular progressive over-relaxation (PSOR) method.

5.2.
Removing the initial discontinuity. The initial condition in (5.1) is not continuous at y = 0, which may cause oscillation in the numerical solution. To remove the discontinuity, we employ a solution of the heat equation with the same initial condition as that in (5.1), which is of the form where N (·) is the cumulative density function of the standard normal distribution, i.e., which has a smooth initial condition. Instead of directly solving forV , we solve for W from (5.4) first and then add back V 0 to getV . 5.3. Numerical example. In order to demonstrate the effectiveness and accuracy of the proposed iterative algorithm. we consider a synthetic example, in which the "true" local volatility is assumed to be σ(S) = −(ln S) 3 /10 + 0.2, and the "market data" are taken to be the solution of (5.1) at τ * corresponding to the "true" volatility.
We use a flat line σ(S) ≡ 0.2 as our initial guess. Figure 1 shows the result obtained through our algorithm after 5 iterations. The blue solid line is the "true" local volatilities and the black dashed line is the recovered ones from our algorithm. We see that satisfactory result is achieved, especially at the near-the-money region. 6. Conclusion. This article suggests a heuristic approach to recover the local volatility from the market data for asset prices which follows a mean-reverting process. To our knowledge, this is the first time that a local volatility recovery problem is considered for a mean-reverting underlying process.
We formulate the problem under an optimal control framework and derive the necessary condition that a minimizer should satisfy. An iterative solving algorithm is given, where the discontinuity in the initial condition is removed, and a numerical example is provided.
Our focus is on suggesting an approach for the reconstruction of local volatility function; as a future direction of research, we hope to establish a rigorous theoretical ground for this approach, e.g., to investigate the existence and uniqueness of the minimizer, etc. Moreover, we shall also seek the extension of our discussion to allow σ be a function of both S and t.