\documentclass[reqno]{amsart} \usepackage{hyperref} \usepackage{graphicx} \AtBeginDocument{{\noindent\small Variational and Topological Methods: Theory, Applications, Numerical Simulations, and Open Problems (2012). {\em Electronic Journal of Differential Equations}, Conference 21 (2014), pp. 183--195. ISSN: 1072-6691. http://ejde.math.txstate.edu, http://ejde.math.unt.edu \newline ftp ejde.math.txstate.edu} \thanks{\copyright 2014 Texas State University - San Marcos.} \vspace{9mm}} \begin{document} \setcounter{page}{183} \title[\hfilneg EJDE-2014/Conf/21 \hfil Inverse volatility for European options] {The inverse volatility problem for European options} \author[I. Knowles, L. Feng, A. Mahato \hfil EJDE-2014/Conf/21\hfilneg] {Ian Knowles, Li Feng, Ajay Mahato} % in alphabetical order \address{Ian Knowles \newline Department of Mathematics, University of Alabama at Birmingham, Birmingham AL 35294, USA} \email{iknowles@uab.edu} \address{Li Feng \newline Department of Mathematics, University of Alabama at Birmingham, Birmingham AL 35294, USA} \email{lifeng@uab.edu} \address{Ajay Mahato \newline Department of Mathematics, University of Alabama at Birmingham, Birmingham AL 35294, USA} \email{amahato7@gmail.com} \thanks{Published February 10, 2014.} \subjclass[2000]{34B24, 65L09, 45J40} \keywords{Inverse volatility; European option; Dupire equation; \hfill\break\indent convex functional} \begin{abstract} The problem of determining equity volatility from a knowledge of European call option prices for a range of exercise (strike) prices and expirations is solved by minimization of a convex functional. \end{abstract} \maketitle \numberwithin{equation}{section} \newtheorem{theorem}{Theorem}[section] \allowdisplaybreaks \section{Introduction} The inner workings of financial markets, from a modeling perspective, are still not well understood, despite more than a century of effort dating back to the pioneering work of Bachelier \cite{bachelier}. Modern physics and its associated PDE modeling is supported by the laws of physics, which have withstood the test of time over centuries. Not so the ``laws of finance'', which appear quite flimsy in comparison. What we do know is that a market is a large collection of people acting individually and collectively, each with their own goals and economic reasons for participating. We also know that in transactions associated with future-oriented instruments, such as stock options and other financial derivatives, a huge amount of data is available buried inside of which is the market's best guess as to what the future holds. We are concerned here with the possibility of extracting information from this type of data with the aid of certain computational inverse algorithms. It is common to model a financial asset (such as a stock or a commodity) via a stochastic differential equation \begin{equation}\label{sde} \frac{dS_t}{S_t} = m (S_t,t) dt + \sigma (S_t,t) d B_t, \end{equation} where, for each time $t$, $S_t(\omega)$ is a random variable representing the price of the financial asset for the trial $\omega$, $m$ is the drift, which relates to the ``trend'' of the asset, $\sigma$ is the volatility (``wobble''), and $B_t(\omega)$ is the Brownian motion stochastic process used to model the randomness. Financial derivatives are contracts that derive their value from such an underlying asset. In particular, a {\it European call option} on a stock is the right to buy one share of the stock at a specified price $K$ (the strike, or exercise, price) at a specified future time $T$ (expiration date). In their Nobel Prize winning paper \cite{black-scholes} Black and Scholes showed that, under certain rather severe restrictions, the arbitrage-free price of a European call contract, $v(S,t)$, satisfies a deterministic PDE of diffusion type in time $t$ and the value $S$ of the underlying asset: \begin{equation}\label{bs} \frac{\partial v}{\partial t} + \frac{1}{2} \sigma^2 S^2 \frac{\partial^2v}{\partial S^2} + \mu s\frac{\partial v}{\partial S} -rv =0, \end{equation} where $\sigma$ is the (assumed constant) volatility, $\mu$ is the risk-neutral drift, and $r$ is the short-term interest rate. For practical purposes the latter may be taken to be the interest rate on a 13-week US government treasury bill. Many authors over the years \cite{fama1965,kon1984,madan-seneta1990,mandelbrot1963a,officer1972,praetz1972, press1967} have noted that several of the assumptions laid down by Black and Scholes are basically incompatible with market data. Notable among these are objections to the constancy of $\sigma$. As recently as the early nineteen eighties, assertions such as ``the Black-Scholes volatility is constant'' seemed to hold true, at least while the market believed it so; but then, after the crash of 1987, volatility was anything but constant, and in fact it has recently become fashionable to speak of (and invest in) the volatility of the volatility! So it is common to regard the volatility as a function of $S$ and $t$, $\sigma=\sigma(S,t)$, and we assume this in the sequel. Now, financial data specifying the market price $v$ of an option is readily available in quantity at various strike values $K$ around the current price of the underlying asset (the ``spot'' price), and for values of the expiration $T$ up to around six months into the future. Given that the computer projections currently used by most stock analysts are only valid for a week or so into the future, one is led quite naturally to the so-called {\it inverse volatility problem}: determine a market-inspired estimate of the future volatility function $\sigma(S,t)$ from a knowledge of current market prices $v$ of options with different strikes and future expirations. The solution of this problem generally goes as follows. The value $v$ of an option contract also depends on the exercise (strike) price, $K$, and the expiration date, $T$, of the contract. In 1994 Bruno Dupire \cite{dupire1994} noticed that the function $v(S,t;K,T)$ satisfies the ``dual'' Black-Scholes equation \begin{equation}\label{dual} \frac{\partial v}{\partial T} - \frac{1}{2} K^2 \sigma(K,T) \frac{\partial^2v}{\partial K^2} + \mu K\frac{\partial v}{\partial K} -(\mu-r)v =0, \end{equation} known also as the Dupire equation. If $v$ is known for all strikes $K$ and expirations $T$ then, as was noted first in \cite{dupire1994}, the volatility is uniquely determined in principle from the equation \eqref{dual}. But such a formula for $\sigma$ is of little use in practice, as the market data for $v$ is not only noisy (which would make the estimation of these derivatives highly ill-posed), but even worse, the data is both discrete in $T$ and somewhat sparse in $K$. A number of alternate approaches have been proposed subsequent to the appearance of \cite{dupire1994}, none of which has offered a definitive solution. Minimization methods using regularized least-squares fitting have been proposed in \cite{avellaneda1997,bodurtha1999,lagnado-osher1997}; the possible presence of spurious local minima is always an issue here. An integral equation approach is presented in \cite{bouchouev-isakov1997,bouchouev-isakov1999}, where convergence problems are possible given the underlying ill-posedness, and in \cite{bouchouev-isakov-valdivia2002,isakov2004} linearization of the inherently non-linear inverse problem is discussed. In this article we present a new variational algorithm for computing, via the Dupire equation \eqref{dual}, the volatility $\sigma(K,T)$ from a knowledge of European option prices at various strikes and expirations. The method used is an adaption of the variational approach involving the minimization of convex functionals (with the associated distinct advantage of having unique global minima and stationary points) used in \cite{kw1} for numerical differentiation (formulated as an inverse problem), and in \cite{kly1,krty1} for solving the inverse groundwater modeling problem. \section{Reconstruction of volatility} For simplicity we assume that there is no dividend for the underlying asset. Thus the risk-neutral drift $\mu$ in \eqref{bs} and \eqref{dual} is equal to the interest rate $r$, and the Dupire equation \eqref{dual} can be written as \begin{equation}\label{dupire} \frac{\partial v}{\partial T} - \frac{1}{2} K^2 \sigma^2 (K,T) \frac{\partial^2v}{\partial K^2} + r K\frac{\partial v}{\partial K} =0. \end{equation} Let $T_0 0$ is a parameter. For each such fixed $i$, $1\le i\le n$, we now Laplace transform the Dupire equation over $[T_{i-1},T_i]$ to obtain \begin{equation*} \int_{T_{i-1}}^{T_i} e^{-\lambda T}v_T\, dT -\frac{1}{2} K^2 \sigma^2_i \underbrace{\int_{T_{i-1}}^{T_i} e^{-\lambda T} v_{KK}\, dT}_{w''_\lambda} +r K \underbrace{\int_{T_{i-1}}^{T_i} e^{-\lambda T}v_K\, dT}_{w'_\lambda} =0, \end{equation*} where the primes indicate differentiation with respect to $K$. On integrating the first term by parts we get \begin{equation*} [e^{-\lambda T} v]_{T_{i-1}}^{T_i} +\lambda \underbrace{\int_{T_{i-1}}^{T_i} e^{-\lambda T}v\, dT}_{w_\lambda} -\frac{1}{2} K^2 \sigma_i^2 w''_\lambda+r K w'_\lambda =0, \end{equation*} and rearranging terms gives, \begin{equation*} -\frac{1}{2} K^2 \sigma_i^2 w''_\lambda + r K w'_\lambda + \lambda w_\lambda = -v(K,T_i)e^{-\lambda T_i}+v(K,T_{i-1})e^{-\lambda T_{i-1}}. \end{equation*} Next, dividing by $\frac{1}{2}K^2\sigma_i^2$ throughout, we obtain \begin{equation*} - (w''_\lambda - \frac{2r}{K \sigma_i^2} w'_\lambda) + \frac{\lambda}{\frac{1}{2}K^2\sigma_i^2} w_\lambda = \frac{-v(K,T_i)e^{-\lambda T_i}+v(K,T_{i-1}) e^{-\lambda T_{i-1}}}{\frac{1}{2}K^2\sigma_i^2}. \end{equation*} Finally, on multiplying by the integrating factor \begin{equation}\label{P} P(K)=e^{-2r\int^K \frac{dk}{k\sigma^2_i(k)}}, \end{equation} we now have an equation in Sturm-Liouville form: \begin{equation}\label{sturm-liouville} -(P(K)w_\lambda')'+\lambda Q(K) w_\lambda = \beta(K,\lambda)Q(K), \end{equation} where \begin{gather} Q(K) = (\frac{2}{K^2\sigma^2_i(K)}) P(K), \label{Q}\\ \beta(K,\lambda)=-v(K,T_i)e^{-\lambda T_i}+v(K,T_{i-1})e^{-\lambda T_{i-1}}. \end{gather} If we can recover the functions $P(K)$ and $Q(K)$ for each $i$, $1\le i\le n$, we can find the volatility {$\sigma_i(K)$} from the formula \begin{equation}\label{sigma} \sigma_i(K)= \sqrt{\frac{2 P(K)}{K^2 Q(K)}}. \end{equation} We now focus attention on a variational approach to the recovery of one such pair of positive coefficient functions $P,Q$ defined on an interval $a\le K\le b$. It is assumed that we are given the functions $w_\lambda(K)$ for $K$ in $[a,b]$ and all $\lambda>0$. For positive functions $p$ and $q$ also defined on $[a,b]$, let $c=(p,q)$. Define {$w_{\lambda,c}(K)$} to be the solution to the boundary value problem \begin{gather} L_{p,\lambda q}w_{\lambda,c}={-( p(K)w_{\lambda,c}')'} +\lambda q(K) w_{\lambda,c} = \beta(K,\lambda)q(K) \label{elliptic_eq},\\ {w_{\lambda,c}(a)=w_\lambda(a),\quad w_{\lambda,c}(b)=w_\lambda(b)} \label{elliptic_eq_bc}. \end{gather} Let $\mathcal{D}$ be the set of all positive function pairs $c=(p,q)$ such that boundary value problem \eqref{elliptic_eq}, \eqref{elliptic_eq_bc} is {\it disconjugate} on [a,b], i.e. every non-trivial solution has at most one zero on [a,b]. It is known \cite[Theorem~6.1, p.~351]{hartman} that \eqref{elliptic_eq} is disconjugate if and only if the boundary value problem \eqref{elliptic_eq}, \eqref{elliptic_eq_bc} can always be solved uniquely. It is also known (c.f. \cite[Proposition~2.1]{kw1}) that this set is open and convex in $\mathcal{L}[a,b]\times\mathcal{L}[a,b]$ and $\mathcal{L}^2[a,b]\times\mathcal{L}^2[a,b]$. For each $\lambda>0$ define the functional $G_\lambda$ on the convex set $\mathcal{D}$ by \begin{equation}\label{functional1} G_\lambda(c)=\int_a^b p(K)(w_\lambda'^2-w_{\lambda,c}'^2) +\lambda q(K) (w_\lambda^2-w_{\lambda,c}^2) -2\beta q(K) (w_\lambda-w_{\lambda,c})\, dK. \end{equation} \section{Properties of the Functional $G_\lambda$} The main properties of the functional $G_\lambda$ are summarized in the following \begin{theorem}\label{th31} \begin{enumerate} \item[(a)] For any $c=(p,q)$ in $\mathcal{D}$, \begin{equation}\label{functional2} G_\lambda(c)=\int_a^b p(w'_\lambda-w'_{\lambda,c})^2 + \lambda q(w_\lambda-w_{\lambda,c})^2. \end{equation} \item[(b)] $G_\lambda(c)\geq0$ for all $c=(p,q)$ in $\mathcal{D}$, and $G_\lambda$(c)=0 if and only if $w_\lambda=w_{\lambda,c}$. \item[(c)] The first G\^{a}teaux derivative of $G_\lambda$ is given by \begin{equation}\label{G'} G_\lambda'(p,q)[h_1,h_2] =\int_a^b \underbrace{(w_\lambda'^2-w_{\lambda,c}'^2)}_{\rm {{L^2\:gradient\:in\:p}}}h_1 +\underbrace{[\lambda(w_\lambda^2-w_{\lambda,c}^2)-2\beta (w_\lambda-w_{\lambda,c})]} _{\rm {{L^2\:gradient\:in\:q}}}h_2. \end{equation} \item[(d)] The second G\^{a}teaux derivative of $G_\lambda$ is given by \begin{equation} G''_\lambda(c)[h,k] = 2(L_{p,\lambda q}^{-1}(e(h)),e(k)), \end{equation} where $h=(h_1,h_2)$ , $k=(k_1,k_2)$, \begin{equation*} e(h)= -(h_1 w'_{\lambda,c})'+\lambda h_2 w_{\lambda,c} - \beta h_2 , \end{equation*} and $(\cdotp,\cdotp)$ denotes the usual inner product in $L^2[a,b]$. \end{enumerate} \end{theorem} \begin{proof} (a) If $v\in W^{1,2}[a,b]$ and $\phi\in W_0^{1,2}[a,b]$ then by integration by parts we have \begin{equation}\label{IBP} \int_a^b p(x) v' \phi'\, dx =\underbrace{p(x) v' \phi |_a^b}_{=0} -\int_a^b \phi (p(x) v')'\,dx = -\int_a^b \phi (p(x) v')'\,dx. \end{equation} Consequently, from \eqref{IBP} using $\phi=w_\lambda-w_{\lambda,c}\in W^{1,2}_0[a,b]$, \begin{align*} G_\lambda(c)&= \int_a^b p(w'^2_\lambda-w'^2_{\lambda,c}) + \lambda q((w_\lambda^2-w_{\lambda,c}^2)-2\beta q(w_\lambda-w_{\lambda,c})\\ &= \int_a^b p(w'_\lambda-w'_{\lambda,c})^2 + 2pw'_{\lambda,c}(w'_\lambda-w'_{\lambda,c}) \\ &\quad + \lambda q((w^2_\lambda-w^2_{\lambda,c}) -2\beta q(w_\lambda-w_{\lambda,c}) \\ &= \int_a^b p(w'_\lambda-w'_{\lambda,c})^2 - 2(w_\lambda-w_{\lambda,c})(pw'_{\lambda,c})' \\ &\quad + \lambda q((w^2_\lambda-w^2_{\lambda,c}) -2\beta q(w_\lambda-w_{\lambda,c}), \\ &\quad \text{using \eqref{elliptic_eq} for $(pw'_{\lambda,c})'$ }, \\ &= \int_a^b p(w'_\lambda-w'_{\lambda,c})^2 - 2(w_\lambda-w_{\lambda,c})(\lambda q w_{\lambda,c} - \beta q)\\ &\quad + \lambda q((w^2_\lambda-w^2_{\lambda,c}) -2\beta q(w_\lambda-w_{\lambda,c}) \\ &= \int_a^b p(w'_\lambda-w'_{\lambda,c})^2 + \lambda q(w_\lambda-w_{\lambda,c})^2, \end{align*} after some rearrangement. (b) As $p$ and $q$ are chosen to be positive and $\lambda>0$, from (a) we get (b). (c) The first G\^{a}teaux derivative of the functional $G_\lambda$ is given by \begin{align*} &G_\lambda'(p,q)[h_1,h_2]\\ &=\lim_{\varepsilon\rightarrow0} \frac{G_\lambda(c+\varepsilon h)- G_\lambda(c)}{\varepsilon} \\ &=\lim_{\varepsilon\rightarrow0}\frac{1}{\varepsilon} \int_a^b (p+\varepsilon h_1) (w_\lambda'^2-w_{\lambda,c+\varepsilon h}'^2)+ \lambda (q+\varepsilon h_2)(w_\lambda^2-w_{\lambda,c +\varepsilon h}^2) \\ &\quad\-2\beta (q+\varepsilon h_2) (w_\lambda-w_{\lambda,c+\varepsilon h}) -p(w_\lambda'^2-w_{\lambda,c}'^2) \\ &\quad-\lambda q(w_\lambda^2-w_{\lambda,c}^2)+2\beta q(w_\lambda-w_{\lambda,c}) \\ &=\lim_{\varepsilon\rightarrow0}\frac{1}{\varepsilon} \int_a^b \varepsilon (w_\lambda'^2-w_{\lambda,c+\varepsilon h}'^2)\, h_1 + \varepsilon \lambda (w_\lambda^2-w_{\lambda,c +\varepsilon h}^2)\, h_2 \\ &\quad- 2 \varepsilon \beta (w_\lambda-w_{\lambda,c})\, h_2 + p (w_\lambda'^2-w_{\lambda,c+\varepsilon h}'^2) \\ &\quad+ \lambda q (w_\lambda^2-w_{\lambda,c+\varepsilon h}^2) - 2 q \beta (w_\lambda-w_{\lambda,c+\varepsilon h}) \\ &=\lim_{\varepsilon\rightarrow0}\int_a^b (w_\lambda'^2-w_{\lambda, c+\varepsilon h}'^2)\, h_1 +[ \lambda (w_\lambda^2-w_{\lambda, c+\varepsilon h}^2)-2\beta (w_\lambda-w_{\lambda,c})]\, h_2 \\ &\quad+\lim_{\varepsilon\rightarrow0} \int_a^b \frac{1}{\varepsilon} p (w_{\lambda,c}'^2-w_{\lambda, c+\varepsilon h}'^2) +\frac{1}{\varepsilon} \lambda q (w_{\lambda,c}^2 -w_{\lambda,c+\varepsilon h}^2)\\ &\quad -\frac{1}{\varepsilon} 2 q \beta (w_{\lambda,c} -w_{\lambda,c+\varepsilon h} ) \end{align*} If we can show the second term is zero we get \eqref{G'}. Let the integral in the second term be denoted by $I$. Now, \begin{gather} {-(pw_{\lambda,c}')'} +\lambda q w_{\lambda,c} = \beta q \label{Lw_c},\\ {-((p+\varepsilon h_1)w_{\lambda,c+\varepsilon h}')'}+\lambda (q+\varepsilon h_2) w_{\lambda,c+\varepsilon h} = \beta (q +\varepsilon h_2) \label{Lw_c+eh}. \end{gather} The first term in the integral $I$ can be expanded as \begin{align*} &\varepsilon^{-1} \int_a^b p (w_{\lambda,c}'^2-w_{\lambda, c+\varepsilon h}'^2)\\ &= \varepsilon^{-1} \int_a^b p(w_{\lambda,c}'+w'_{\lambda, c+\varepsilon h})(w_{\lambda,c}'-w'_{\lambda,c+\varepsilon h}), \\ &\quad \text{from \eqref{IBP} using $\phi=w_{\lambda,c}-w_{\lambda, c+\varepsilon h}$ } ,\\ &= \varepsilon^{-1} \int_a^b (w_{\lambda,c+\varepsilon h}-w_{\lambda, c})(p(w_{\lambda,c}'+w'_{\lambda,c+\varepsilon h}))'\\ &= \varepsilon^{-1} \int_a^b (w_{\lambda,c+\varepsilon h}-w_{\lambda, c})[(pw_{\lambda,c}')'+(pw'_{\lambda,c+\varepsilon h})'], \\ &\quad \text{using \eqref{Lw_c} and \eqref{Lw_c+eh}}, \\ &= \varepsilon^{-1} \int_a^b (w_{\lambda,c+\varepsilon h}-w_{\lambda,c}) [\lambda q w_{\lambda,c} - \beta q + \lambda (q+\varepsilon h_2) w_{\lambda,c+\varepsilon h} \\ &\quad -\beta (q +\varepsilon h_2) -\varepsilon(h_1 w'_{\lambda, c+\varepsilon h})'] \\ &= \int_a^b (w_{\lambda,c+\varepsilon h}-w_{\lambda,c}) [ \lambda h_2 w_{\lambda,c+\varepsilon h} -\beta h_2 -(h_1 w'_{\lambda,c+\varepsilon h})'] \\ &\quad+\varepsilon^{-1} \int_a^b (w_{\lambda,c+\varepsilon h} -w_{\lambda,c})[\lambda q (w_{\lambda,c} + w_{\lambda,c+\varepsilon h}) -2\beta q ] \\ &= \int_a^b (w_{\lambda,c+\varepsilon h}-w_{\lambda,c}) [ \lambda h_2 w_{\lambda,c+\varepsilon h} -\beta h_2 -(h_1 w'_{\lambda,c+\varepsilon h})'] \\ &\quad +\varepsilon^{-1} \int_a^b \lambda q(w_{\lambda, c+\varepsilon h}^2-w_{\lambda,c}^2)(-2\beta q (w_{\lambda,c+\varepsilon h} -w_{\lambda,c}) ). \end{align*} Substituting the above for $\varepsilon^{-1} \int_a^b p (w_{\lambda,c}'^2 -w_{\lambda,c+\varepsilon h}'^2)$ in $I$ we obtain \[ I=\int_a^b (w_{\lambda,c+\varepsilon h}-w_{\lambda,c}) [ \lambda h_2 w_{\lambda,c+\varepsilon h} -\beta h_2 -(h_1 w'_{\lambda,c+\varepsilon h})']. \] It follows that $I \rightarrow 0$ as $\varepsilon \rightarrow 0$. (d) To find the second G\^{a}teaux derivative of the functional $G_\lambda$ we will need the following result: \begin{equation} \begin{aligned} L_{p,\lambda q}(w_{\lambda,c+\varepsilon h} -w_{\lambda,c}) &= -(p(w_{\lambda,c+\varepsilon h} -w_{\lambda,c})')' + \lambda q(w_{\lambda,c+\varepsilon h} -w_{\lambda,c}) \\ &=-(pw'_{\lambda,c+\varepsilon h})'+ \lambda q w_{\lambda,c+\varepsilon h} -[-(p w_{\lambda,c}')' + \lambda q w_{\lambda,c}], \\ &\text{\quad using \eqref{Lw_c} and \eqref{Lw_c+eh}}, \\ &=\varepsilon [(h_1 w'_{\lambda,c+\varepsilon h})'-\lambda h_2 w_{\lambda,c+\varepsilon h} + \beta h_2] \end{aligned}\label{L_pq} \end{equation} The second G\^{a}teaux derivative of the functional $G_\lambda$ is given by \begin{align*} & G''_\lambda(c)[h,k]\\ &= \lim_{\varepsilon\rightarrow0}\frac {G'(c+\varepsilon h)[k]-G'(c)[k]}{\varepsilon}\\ &= \lim_{\varepsilon\rightarrow0}\frac{1}{\varepsilon} \int_a^b(w_\lambda'^2-w_{\lambda,c+\varepsilon h}'^2)k_1 +[\lambda(w_\lambda^2-w_{\lambda,c+\varepsilon h}^2)-2\beta (w_\lambda-w_{\lambda,c+\varepsilon h})]k_2\\ &\quad -(w_\lambda'^2-w_{\lambda,c}'^2)k_1 -[\lambda(w_\lambda^2-w_{\lambda,c}^2)-2\beta (w_\lambda-w_{\lambda,c})]k_2 \\ &= \lim_{\varepsilon\rightarrow0} \frac{1}{\varepsilon} \int_a^b (w_{\lambda,c}'^2-w_{\lambda,c+\varepsilon h}'^2) k_1 +[\lambda (w_{\lambda,c}^2-w_{\lambda,c+\varepsilon h}^2) -2\beta(w_{\lambda,c}-w_{\lambda,c+\varepsilon h})] k_2\\ &= \lim_{\varepsilon\rightarrow0} \frac{1}{\varepsilon} \int_a^b k_1(w_{\lambda,c}'+w_{\lambda,c+\varepsilon h}') (w_{\lambda,c}'-w_{\lambda,c+\varepsilon h}') \\ &\quad +[\lambda (w_{\lambda,c}^2-w_{\lambda,c+\varepsilon h}^2) -2\beta(w_{\lambda,c}-w_{\lambda,c+\varepsilon h})] k_2, \\ &\quad \text{from \eqref{IBP} using $\phi=w_{\lambda,c}-w_{\lambda,c+\varepsilon h}$ }, \\ &= \lim_{\varepsilon\rightarrow0} \frac{1}{\varepsilon} \int_a^b (w_{\lambda,c}-w_{\lambda,c+\varepsilon h}) (-k_1(w_{\lambda,c}'+w_{\lambda,c+\varepsilon h}') )'\\ &\quad +[\lambda (w_{\lambda,c}^2-w_{\lambda,c+\varepsilon h}^2) -2\beta(w_{\lambda,c}-w_{\lambda,c+\varepsilon h})] k_2, \\ &\quad \text{factoring $(w_{\lambda,c}-w_{\lambda,c+\varepsilon h})$ }, \\ &= \lim_{\varepsilon\rightarrow0} \frac{1}{\varepsilon} \int_a^b (w_{\lambda,c}-w_{\lambda,c+\varepsilon h})\, [(-k_1(w_{\lambda,c}'+w_{\lambda,c+\varepsilon h}') )'\\ &\quad+(\lambda (w_{\lambda,c}+w_{\lambda,c+\varepsilon h}) -2\beta)k_2], \\ &\quad \text{using \eqref{L_pq}}, \\ &= \lim_{\varepsilon\rightarrow0} \int_a^b L_{p,\lambda q}^{-1} [-(h_1 w'_{\lambda,c+\varepsilon h})'+\lambda h_2 w_{\lambda,c+\varepsilon h}- \beta h_2] \\ &\quad\times [(-k_1(w_{\lambda,c}'+w_{\lambda,c+\varepsilon h}') )' +(\lambda (w_{\lambda,c}+w_{\lambda,c+\varepsilon h}) -2\beta)k_2]\\ &= \lim_{\varepsilon\rightarrow0} \int_a^b L_{p,\lambda q}^{-1} [-(h_1(w'_{\lambda,c+\varepsilon h}-w'_{\lambda,c}))'+\lambda h_2(w_{\lambda,c+\varepsilon h}-w_{\lambda,c})] \\ &\quad\times [(-k_1(w_{\lambda,c}'+w_{\lambda,c+\varepsilon h}') )' +(\lambda (w_{\lambda,c}+w_{\lambda,c+\varepsilon h}) -2\beta)k_2]\\ &\quad+ \lim_{\varepsilon\rightarrow0} \int_a^b L_{p,\lambda q}^{-1} [-(h_1 w'_{\lambda,c})'+\lambda h_2 w_{\lambda,c}) - \beta h_2] \\ &\quad\times [(-k_1(w_{\lambda,c}'+w_{\lambda,c+\varepsilon h}') )' +(\lambda (w_{\lambda,c}+w_{\lambda,c+\varepsilon h}) -2\beta)k_2], \\ &\quad \text{expanding the second integral},\\ &= \lim_{\varepsilon\rightarrow0} \int_a^b L_{p,\lambda q}^{-1} [-(h_1(w'_{\lambda,c+\varepsilon h}-w'_{\lambda,c}))'+\lambda h_2(w_{\lambda,c+\varepsilon h}-w_{\lambda,c})] \\ &\quad\times[(-k_1(w_{\lambda,c}'+w_{\lambda,c+\varepsilon h}') )' +(\lambda (w_{\lambda,c}+w_{\lambda,c+\varepsilon h}) -2\beta)k_2]\\ &\quad + \lim_{\varepsilon\rightarrow0} \int_a^b L_{p,\lambda q}^{-1} [-(h_1 w'_{\lambda,c})'+\lambda h_2 w_{\lambda,c}) - \beta h_2] \\ &\quad\times [(-k_1(w_{\lambda,c+\varepsilon h}'-w_{\lambda,c}') )' +\lambda (w_{\lambda,c+\varepsilon h}-w_{\lambda,c})k_2]\\ &\quad + 2\int_a^b L_{p,\lambda q}^{-1} [-(h_1 w'_{\lambda,c})' +\lambda h_2 w_{\lambda,c}- \beta h_2] \\ &\quad\times [-(k_1 w'_{\lambda,c})'+\lambda k_2 w_{\lambda,c} - \beta k_2]. \end{align*} The first and second terms equal zero. Thus we obtain \begin{equation} G''_\lambda(c)[h,k] = 2(L_{p,\lambda q}^{-1}(e(h)),e(k)), \end{equation} where \begin{gather*} e(h)= -(h_1w'_{\lambda,c})'+\lambda h_2 w_{\lambda,c}- \beta h_2, \\ e(k)= -(k_1w'_{\lambda,c})'+\lambda k_2 w_{\lambda,c}- \beta k_2. \end{gather*} This completes the proof of the theorem. \end{proof} With some additional work one can show that the first and second G\^{a}teaux derivatives of $G_\lambda$ are also Fr\'{e}chet derivatives. As $L_{p,\lambda q}$ is a positive operator on $W^1_0[a,b]$, we have from Theorem~\ref{th31}(d) that $G''_\lambda(c)\geq0$ for all $c$ in the convex set $\mathcal{D}$. By \cite[Corollary~42.8]{Zeidler3} the functional $G_\lambda$ is therefore convex on $\mathcal{D}$. We know from Theorem~\ref{th31}(b) that $G_\lambda$ has a global minimum (zero) at $c=(p,q)$ if and only if $w_\lambda=w_{\lambda,c}$. Choose $N\ge 3$ positive distinct real numbers $\lambda_j$, $1\le j\le N$, so that $$ 0<\lambda_j T<2,\quad T\in[T_{i-1},T_i]. $$ Define a convex functional $G$ on the domain $\mathcal{D}$ (defined above) by \begin{equation} G(c)= \sum_{j=1}^{N} G_{\lambda_j}(c). \end{equation} From the uniqueness theorem \cite[Theorem~3.5]{k3} we know that, under certain (computer-verifiable) conditions on the nature of the flows of certain associated vector fields (which amount here to an admissibility restriction on the data $v(K,T)$), the condition $w_\lambda=w_{\lambda,c}$ for at least three distinct values of $\lambda$ implies that $c=(p,q)=(P,Q)$. By \cite[Proposition~42.6(1)]{Zeidler3} we know that if the convex functional $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 process getting stuck in spurious local minima. That this cannot happen here is one of the significant advantages of our approach. \section{The Algorithm} $G(c)$ is a nonnegative convex functional since it is the sum of nonnegative convex functionals, and it also has a unique stationary point at $c=(P,Q)$. The idea here is that by using $G$ rather than just one of the $G_\lambda$, 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. We note in passing from \cite{klar1} that this inverse recovery is conditionally well-posed in the weak-$L^2$ sense, so from a theoretical standpoint, the recoveries are expected to be quite stable, which indeed is the case. We minimize this functional for $N=20$ using the steepest descent method to recover the coefficients $P(K)$ and $Q(K)$. The $L^2$-direction of steepest descent for $G$ at $c_0=(p_0,q_0)$ with respect to $p$ is $$ -\nabla_{L^2,p} G(c_0) = \sum_{j=1}^N\,(w_{\lambda_j}'^2-w_{\lambda_j,c_0}'^2), $$ and the $L^2$-direction of steepest descent for $G$ at $(p_0,q_0)$ with respect to the variable $q$ is given by $$ -\nabla_{L^2,q} G(c_0) = \sum_{j=1}^N\,[\lambda_j(w_{\lambda_j}^2 -w_{\lambda_j,c_0}^2)-2\beta (w_{\lambda_j}-w_{\lambda_j,c_0})]. $$ Instead of using these $L^2$-gradients we use the corresponding Neuberger-gradients (see \cite{neuberger}) as the $L^2$-gradient has numerical problems that are extensively discussed in \cite{kw1}. In particular, the $L^2$-gradient with respect to $q$ is zero on the boundary of [a,b] given that $w_\lambda$ and $w_{\lambda,c}$ are equal there, and thus the algorithm is unable to properly recover $Q$. The Neuberger-gradient smooths the $L^2$-gradient and preserves boundary data during the descent, an important property not shared by other descent techniques. Our Neuberger-gradient $g=\nabla_{H^1}G$ can be found from an $L^2$-gradient $\nabla_{L^2}G$ by solving the boundary value problem \begin{equation} \begin{gathered} -g''+g=\nabla_{L^2} G,\\ g(a)=g(b)=0. \end{gathered} \end{equation} Below is the steepest descent algorithm used to get one descent step in $p$: \begin{enumerate} \item Initialize ${p(K)}$ and ${q(K)}$ with $c_0=(p_0,q_0)$. \item Find $w_{\lambda,c_0}$ and $w'_{\lambda,c_0}$ by solving \eqref{elliptic_eq},\eqref{elliptic_eq_bc}. \item Find the $L^2$ gradient of $G$ in $p$, $\nabla_{L^2,p} G(c_0)$. \item Find the Neuberger gradient in $p$, {$\nabla_{H^1,p} G( c_0)$}. \item Evaluate {$ p_{new}(K) = p_0(K) - \alpha \nabla_{H^1,p} G( c_0) $}. \item Find ${G(p,q_0)}$ using {$ p_{new}(K) $} for {$p(K)$}. \item Find ${\alpha}$ that gives the lowest value of ${G(p,q_0)}$. \item Set {$ p(K)= p_{new}(K)$}. \end{enumerate} The descent in $q$ is similar to that of descent in $p$. Here we find corresponding gradients in $q$. The $q_{new}(K)$ is given by $$ {q_{new}(K) = q(K) - \alpha \nabla_{H^1,q} G(c_0)} $$ The specific order of descent is somewhat problem dependent, and different combinations of descents in $p$ and $q$ were tried to get the best minimization. Typically one needs more $p$-descent steps relative to $q$-descent steps as the descent progresses. \section{Results} One of the most popular European options traded on US exchanges is the option on the Standard \& Poors 500 (SPX) index. Call option prices on the SPX index were taken from the official website of the Chicago Board Options Exchange (CBOE), for two consecutive maturities on the 22nd of February, 2012. The data includes only the near-the-money options as they are the most heavily traded. To recover the coefficient functions $P(K)$ and $Q(K)$ in \eqref{sturm-liouville} a computer code code was written in the programming language C. The volatility recovered was compared to the ``implied volatility'' obtained directly from the standard formula of Black and Scholes by substituting the known option price and solving for the implied volatility $\sigma$ as an unknown. \begin{figure}[ht] \begin{tabular}{ | c | r | r |} \hline \multicolumn{2}{|c|}{Spot Price($S_0$)} & \$ 1357.66 \\ \hline \multicolumn{2}{|c|}{Maturity Time ($T_1$)} & 2 days \\ \hline \multicolumn{2}{|c|}{Maturity Time ($T_2$)} & 22 days \\ \hline Strike Price($K$) & $v(K,T_1)$ & $v(K,T_2)$ \\ \hline 1300 & 59.40 & 63.00 \\ \hline 1305 & 54.40 & 58.60 \\ \hline 1310 & 49.20 & 54.20 \\ \hline 1315 & 44.60 & 50.00 \\ \hline 1320 & 39.40 & 45.80 \\ \hline 1325 & 34.80 & 41.00 \\ \hline 1330 & 29.80 & 37.90 \\ \hline 1335 & 25.10 & 34.10 \\ \hline 1340 & 20.60 & 30.50 \\ \hline 1345 & 16.40 & 27.00 \\ \hline 1350 & 12.5 & 23.7 \\ \hline 1355 & 7.9 & 20.6 \\ \hline 1360 & 5.1 & 17.7 \\ \hline 1365 & 2.85 & 14.5 \\ \hline 1370 & 1.6 & 12.7 \\ \hline 1375 & 0.9 & 10 \\ \hline 1380 & 0.55 & 8.7 \\ \hline 1385 & 0.3 & 7 \\ \hline 1390 & 0.25 & 5.6 \\ \hline 1395 & 0.2 & 4.5 \\ \hline 1400 & 0.2 & 3.6 \\ \hline \end{tabular} \end{figure} We have option prices for discrete sets of strikes and expirations. We generated the function $v(K,T)$ by linearly interpolating the option price in both strike and expiration. The function $v(K,T)$ was mollified (c.f. \cite[\S6]{kly1}) so that it could be differentiated, and the derivative $v_K(K,T)$ was found using central differences. For 20 fixed values of $\lambda$ the functions $v(K,T)$ and $v_K'(K,T)$ were Laplace transformed using \eqref{Laplace_tranform} to $w_\lambda(K)$ and $w'_\lambda (K)$ respectively. The functions $p(K)$ and $q(K)$ were initialized using \eqref{P} and \eqref{Q} with the initial $\sigma_i$ chosen to be the implied volatility. We performed a series of descents in $p$ using the aforementioned Neuberger steepest descent algorithm such that the functional could not be minimized any further. Then a series of descents in $q$ were performed to the point where functional likewise could not be lowered any further. We repeated this sequence of descents in $p$ and $q$. The minimization of $G(c)$ in $\alpha$ was done using the well known Brent minimization technique, by adapting the one-variable code in the Numerical Recipes in C function {\tt brent()}. To avoid possible catastrophic cancellation in the Simpson rule formula used in the calculation of the integrals in the formula \eqref{functional1} for the functional $G_\lambda$, we used the alternate formula \eqref{functional2} instead. After running the code we recovered the functions $P(K)$ and $Q(K)$ graphed below. \begin{figure}[ht] \begin{center} \includegraphics[width=0.7\textwidth]{fig1} % P.pdf \end{center} \end{figure} \begin{figure}[ht] \begin{center} \includegraphics[width=0.7\textwidth]{fig2} % Q.pdf \end{center} \end{figure} From \eqref{sigma} we calculated the volatility and compared it to the implied volatility of the option at first and second expirations, as shown in the graph below. On taking subsequent maturity intervals a volatility surface can in principle be plotted. \begin{figure}[ht] \begin{center} \includegraphics[width=0.7\textwidth]{fig3} % recovered_n_implied_vol.pdf \end{center} \end{figure} Finally, from the recovered volatility we calculated the option price in MATLAB using the Binomial method and compared it to the actual price, as shown in the figure below. \begin{figure}[ht] \begin{center} \includegraphics[width=0.7\textwidth]{fig4} % compare_option_price_at_first_expiry.pdf \end{center} \end{figure} \subsection*{Conlusion} We have shown that volatility can be recovered from published option prices using a steepest descent minimization technique. This provides a ``market view'' of future volatility which in principle can used to trade options more efficiently. The results obtained look promising. The analogous work on recovering volatility for the much more ubiquitous American options is in progress. It would be interesting to consider interest rate $r$ (also known in this context as the risk-neutral drift) as function of time and asset price, instead of treating it as a constant, and recover it in similar fashion from the option price. \begin{thebibliography}{10} \bibitem{avellaneda1997} {M.} Avellaneda, {C.} Friedman, {L.} Holmes, and {L.} Sampieri; \newblock Calibrating volatility surfaces via relative entropy minimization. \newblock {\em Appl. Math. Finance}, 4:37--64, 1997. \bibitem{bachelier} Louis Bachelier; \newblock {\em Th\'{e}orie de la Sp\'{e}culation}. \newblock PhD thesis, L'\'{E}cole Normale Sup\'{e}riere, 1900. \newblock Translation: Cootner, 1964. \bibitem{black-scholes} Fischer Black and Myron Scholes; \newblock The pricing of options and corporate liabilities. \newblock {\em Journal of Political Economy}, 81:637--654, 1973. \bibitem{bodurtha1999} {J. N.} Bodurtha and {M.} Jermakyan; \newblock Non-parameric estimation of an implied volatility surface. \newblock {\em J. Computational Finance}, 2:29--61, 1999. \bibitem{bouchouev-isakov1997} Ilia Bouchouev and Victor Isakov; \newblock The inverse problem of option pricing. \newblock {\em Inverse Problems}, 13:L11--L17, 1999. \bibitem{bouchouev-isakov1999} Ilia Bouchouev and Victor Isakov; \newblock Uniqueness, stability and numerical methods for the inverse problem that arises in financial markets. \newblock {\em Inverse Problems}, 15:R95--R116, 1999. \bibitem{bouchouev-isakov-valdivia2002} Ilia Bouchouev and Victor Isakov; \newblock Recovery of volatility coefficient by linearization. \newblock {\em Quant. Finance}, 2:257--263, 2002. \bibitem{dupire1994} {B.} Dupire; \newblock Pricing with a smile. \newblock {\em RISK}, 7:18--20, 1994. \bibitem{fama1965} {Eugene F.} Fama; \newblock The behavior of stock-market prices. \newblock {\em Journal of Business}, 38:34--105, 1965. \bibitem{hartman} Philip Hartman; \newblock {\em Ordinary differential equations}. \newblock S. M. Hartman, Baltimore, Md., 1973. \newblock Corrected reprint. \bibitem{isakov2004} Victor Isakov; \newblock The inverse problem of option pricing. \newblock Preprint, 2004. \bibitem{k3} Ian Knowles; \newblock Uniqueness for an elliptic inverse problem. \newblock {\em SIAM J. Appl. Math.}, 59(4):1356--1370, 1999. \bibitem{klar1} Ian Knowles, Mary A. LaRussa; \newblock Conditioinal well-posedness for an elliptic inverse problem. \newblock {\em SIAM J. Appl. Math.}, 71:952--971, 2011. \newblock Available online at http://www.math.uab.edu /knowles/pubs.html. \bibitem{kly1} Ian Knowles, {Tuan A.} Le, Aimin Yan; \newblock On the recovery of multiple flow parameters from transient head data. \newblock {\em J. Comp. Appl. Math.}, 169:1--15, 2004. \bibitem{krty1} Ian Knowles, Michael Teubner, Aimin Yan, Paul Rasser, {Jong Wook} Lee; \newblock Inverse groundwater modelling in the {W}illunga {B}asin, {S}outh {A}ustralia. \newblock {\em Hydrogeology Journal}, 15:1107--1118, 2007. \bibitem{kw1} Ian Knowles, Robert Wallace; \newblock A variational method for numerical differentiation. \newblock {\em Numerische Mathematik}, 70:91--110, 1995. \bibitem{kon1984} {S. J.} Kon; \newblock Models of stock returns - a comparison. \newblock {\em Journal of Finance}, 39(1):147--165, 1984. \bibitem{lagnado-osher1997} {R.} Lagnado and {S.} Osher; \newblock A technique for calibrating derivation of the security pricing models: numerical solution of the inverse problem. \newblock {\em J. Computational Finance}, 1:13--25, 1997. \bibitem{madan-seneta1990} {D. B.} Madan and {E.} Seneta; \newblock The variance gamma model for share market returns. \newblock {\em Journal of Business}, 63(4):511--524, 1990. \bibitem{mandelbrot1963a} {Benoit B.} Mandelbrot. \newblock The variation of certain speculative prices. \newblock {\em Journal of Business}, 36:394--41, 1963. \bibitem{neuberger} J.~W. Neuberger; \newblock {\em Sobolev gradients and differential equations}, volume 1670 of {\em Lecture Notes in Mathematics}. \newblock Springer-Verlag, Berlin, second edition, 2010. \bibitem{officer1972} {R. R.} Officer; \newblock The distribution of stock returns. \newblock {\em Journal of the American Statistical Association}, 67(340):807--812, 1972. \bibitem{praetz1972} {P. D.} Praetz. \newblock The distribution of share price changes. \newblock {\em Journal of Business}, 45(1):49--55, 1972. \bibitem{press1967} {S. J.} Press; \newblock A compound events model for security prices. \newblock {\em Journal of Business}, 40(July):317--335, 1967. \bibitem{Zeidler3} Eberhard Zeidler; \newblock {\em Nonlinear functional analysis and its applications. {III}}. \newblock Springer-Verlag, New York, 1985. \newblock Variational methods and optimization, Translated from the German by Leo F. Boron. \end{thebibliography} \end{document}