An inverse problem of reconstructing option drift rate from market observation data

Drift rate is a very important parameter in the evolution of stock price, which has significant impact on the corresponding option pricing. This paper deals with an inverse problem of recovering the drift function by current market prices of options. Different from the usual inverse volatility problem, our mathematical model does not tend to zero at infinity, which may bring great trouble to theoretical analysis and numerical calculation. To overcome this difficulty, we use an artificial boundary and homogenization technique to transform the original problem into a homogeneous initial boundary value problem on a bounded domain. Then, based on the optimal control framework, we construct the corresponding optimization problem and strictly prove the well-posedness of the minimizer. Finally, we design an iterative algorithm to obtain the numerical solution. We give some typical examples to verify the validity of our method.

In general, the underlying asset S t at time t is modeled by the stochastic differential equation where the process W t is the standard Brownian motion. The parameters μ and σ are called the real drift and local volatility of the underlying asset, respectively. The drift term μ indicates the expected return of stock price changes, whereas the volatility is used to measure the variability of variables. The stock price movement is determined by the drift rate, volatility, and Brownian motion. The randomness of volatility, drift rate, and Brownian motion determines the process of stock price, which is full of randomness and uncertainty. In practice, the drift rate is quite difficult to measure, but it has an important impact on the stock price trend. Therefore it is of great financial significance to use some indirect approaches to estimate the drift rate, which represents people's prediction of the expected return rate of stock price.
Black and Scholes [1] first discovered how to construct a dynamic portfolio t of a derivative security and the underlying asset. By Itô's lemma the stochastic behavior of the derivative security u(t, S) is driven by the stochastic differential equation In the absence of arbitrage opportunities, the instantaneous return of this portfolio must be equal to the interest rate r > 0, that is, the return of a riskless asset such as a bank deposit. Therefore this equality takes the form of the following partial differential equation (the Black-Scholes equation): where r and the dividend rate δ are known constants. However, the theoretical prices of options with different strike prices calculated by the Black-Scholes model differ from real market prices. Under the noarbitrage property of the financial market, the real drift μ does not enter the above equation. In [23], taking this into account, the following new binary option model is derived: This model is a form of arbitrage model. Further, considering the property of binary option, the final condition at the maturity is specified by We would like to determine the drift function μ using the current market prices of options with different strikes K and fixed maturity T.
Using the Dupire technique, we deduce that the price u(T, K) of the binary option with maturity T and strike price K satisfies the adjoint equation: Making the changes in variablesW (τ , y) = u(T, K), y = ln(K/S * ), τ = Tt, we identify the drift a(y) := μ(K) satisfying the following:

Problem P1
The Cauchy problem of second-order parabolic equation where the local volatility σ 0 is a constant, a(y) is an unknown coefficient in (1.2). The extra Determine the functionsW and a(y) satisfying (1.2)-(1.3).
In this problem, y ∈ R, so the problem is an unbounded domain problem which is not conducive to numerical calculations. With this in mind, we translate the problem into an approximation problem on a larger bounded domain y ∈ [-L, L], where L is a large positive constant. The problem can be rewritten in the following form.

Problem P2
The initial-boundary value problem of the second-order parabolic equation The boundary conditions in P2 is nonhomogeneous, which is not conductive to integration by parts. We try to convert the nonhomogeneous equation to a homogeneous one. Let Then Problem P2 is transformed into the following form.
Problem P The initial-boundary value problem of second-order parabolic equation (1.6) The additional condition is given by Inverse coefficient problems for parabolic equations are well studied in the literature, and abundant theoretical and numerical results are obtained. The inverse problem of identifying coefficient q(y) in the parabolic equation from final overdetermination data u(y, T) has been investigated by several authors, for example, in [5,6,17,18,25,27,29,31]. The purely time-dependent case, that is, determining the unknown radiative coefficient in heat conduction equations independent of the spatial variable, has been extensively studied by several authors (see, e.g., [3, 4, 8-11, 24, 28]). For the space-and time-dependent case, we refer the readers to [12,13].
In financial mathematics, there is often a kind of inverse problems of using market observation data to reconstruct the implied volatility. Using the optimal control theory, the recovery of volatility in the Black-Scholes equation is studied from the current options market in [21]. In [2] the authors reduce the identification of volatility to an inverse parabolic problem with terminal observation and establish uniqueness and stability results by using the Carleman estimates. In [7] the local volatility surface is recovered by nonlinear Landweber iteration using simulated data. A new continuous-time model is proposed in [19] to recover the volatility, and the corresponding numerical results are obtained by solving a couple of fully nonlinear parabolic equations. In [16] the volatility is parameterized by five special numbers, and (nonlinear) minimization of the mismatch functional is implemented. Compared with the inverse volatility problem, there are few documents on the inverse problem of drift rate. In [23] the authors consider an inverse problem of recovering the real drift of binary call options from market prices. By using the linearization method the inverse problem is transformed to an integral equation, and numerical results are also obtained.
In this paper, we use the optimal control method to discuss Problem P. Compared with other papers concerning volatility identification problems, our work has the following unusual features. First, the unknown function to be identified is the first-order derivative coefficient rather than the principal one. Second, our mathematical model does not tend to zero at infinity, which may bring great trouble to theoretical analysis and numerical calculation.
The rest of the paper is organized as follows. In Sects. 2 and 3, we transform Problem P into an optimal control Problem P3 and prove the existence of the minimum for the control functional. The necessary condition, which must be satisfied by the minimum, is deduced in Sect. 4. In Sect. 5, we prove that the minimum is locally unique under some assumptions. In Sect. 6, we design an iterative algorithm to obtain the numerical solution and give some typical examples. Section 7 ends this paper with a summary.

Optimal control problem
Consider the following optimal control problem.
is the solution to problem (1.6) for given a ∈ A, and N is the regularization parameter.

Lemma 2.1 If U(y, τ ) is a solution to the initial-boundary value problem
Integrating by parts, we obtain From the Cauchy-Schwarz inequality we have This completes the proof of Lemma 2.1.

Existence Theorem 3.1 There exists a minimizerā ∈ A of J(a), that is,
Proof Let (U n , a n ) be a minimizing sequence. Since J(a n ) ≤ C, thanks to the particular structure of J, we deduce ∇a n L 2 (-L,L) ≤ C (C is independent of n).
By the Sobolev imbedding theorem we obtain a n C 1/2 (-L,L) ≤ C. Thus Therefore we can select subsequences of a n and U n , again denoted by a n and U n , such that a n (y) →ā( We easily check that (ā(y),Ū(y, τ )) satisfies (1.6). By the Lebesgue control convergence theorem and the weak semicontinuity of the L 2 norm we obtain This completes the proof of Theorem 3.1.

Necessary condition
Proof For any h ∈ A and 0 ≤ δ ≤ 1, we have Let U δ be the solution to the initial-boundary value problem (4.1) with given a = a δ . Since a is an optimal solution, by (2.2) we have Denoting U δ ≡ ∂U δ ∂δ and using (4.1), by direct calculation we get the following equation: Suppose V is the generalized solution of the following adjoint problem: (4.9) Combining (4.7) and (4.9), we easily obtain that This completes the proof of Theorem 4.1. Let a 1 (y) and a 2 (y) be two minimizers of the control Problem P3, and let {U i , V i } (i = 1, 2) be solutions of system (4.1)-(4.2) in whichā = a i , respectively. Set

Uniqueness
Then U and V satisfy Using the maximum principle, we obtain V i ∞ ≤ U i (y, τ * ) -U * i (y) ∞ (i = 1, 2). This completes the proof of Lemma 5.2.
Using the boundedness of a 1 (y), we get the following inequality: that is, This completes the proof of Lemma 5.3.

By Lemma 5.3 this yields
This completes the proof of Lemma 5.4.

Theorem 5.5
Let a 1 (y) and a 2 (y) be two minimizers of the optimal control Problem P3. If there exists a point y 0 such that a 1 (y 0 ) = a 2 (y 0 ), then for τ * 1, we have a 1 (y) ≡ a 2 (y) for any y ∈ [-L, L].

Numerical examples
In this section, we give some numerical examples to test the validity of the proposed methods for reconstruction of the drift. In this paper, we consider the gradient iteration algorithm to obtain the numerical solutions. The key ingredient of this iteration algorithm is the Gâteaux derivative of J (a), which is given as follows.
The proof is similar to that of Theorem 4.1.
Remark 6.1 The main difference between Problems P2 and P is the boundary condition. Problem P is homogeneous, whereas Problem P2 is nonhomogeneous. Problem P is convenient for theoretical analysis, such as integration by parts, but there is no difference between the two problems in calculation. In this section, we use Eq. where M and P are two integers.
Based on the above analysis, the detailed procedure of iteration algorithm can be summarized as follows: Step 1. Choose an initial value of iteration a = a 0 (y).
Step 4. Compute the Gâteaux derivative J (a)ψ j = c j for j = 0, 1, 2, . . . , 2M, where the functions ψ j are taken as the base function under current grid:  Then the iteration direction from the jth step to the (j + 1)th step is given by c j ψ j (y).
Step 5. Compute the norm of C 0 (y) at the jth step: where h is the spatial step size.
Step 6. Choose an arbitrary small positive constant ε as the stopping parameter. Go on or stop the iteration is determined by the following steps: Step 6.1. Let k = 1.
Step 6.3. Compare it with 0; if err ≤ 0, then go to Step 6.4; Otherwise, let k . = μk and go to Step 6.2, where μ is an adjusting parameter.
Step 6.4. Take a 1 (y) = a 0 (y) + kC 0 (y); if kC 0 (y) ≤ ε, then exit and stop the iteration scheme. Otherwise, set j = j + 1 and go to Step 2. Let a 1 (y) be the new initial value of iteration and go on computing by the induction rules until the iterations meet the termination conditions.
We have performed three numerical experiments to check the stability of our iteration algorithm. Since real-world data are not available, we would like to use artificial data to test the stability of the proposed numerical algorithm, that is, the extra condition is obtained by solving the direct problem. In all experiments, we used the basic parameters r = 0.5, σ 0 = 1, ε = 10 -4 , N = 10 -5 .
Equation (1.4) is solved by the classical finite difference method, the Crank-Nicolson difference scheme: The discrete equations of the initial boundary value are The difference scheme is absolutely stable, and its truncation error is O(τ 2 + h 2 ).
Example 1 In the first numerical experiment, we take L = 10 The exact solution and the reconstruction results for different time (denoted by τ * ) are shown in Fig. 1. The initial guess is taken to be 0.5. We can see that the drift coefficient a(y) is well recovered after 260 iterations. The iterative procedure converges quickly, and the effect is satisfactory. However, since a(y) is a segmented function, the values on the right boundary are quite difficult to reconstruct well. The algorithm converges rather slowly near the right boundary. Also, since the prices near the strike are the most interesting for practitioners, we investigate the recovered function around y = 0 (S = K ), where we use three different observation times τ * = 0.5, 0.8, and 1. From the figure we observe that the reconstruction is numerically good when τ * = 1. For all times τ * , the reconstructed a(y) is near-perfect on interval y ∈ (0, 2), as shown in Fig. 1.  The initial guess is taken to be 0. We only needed 380 iterations to achieve a satisfactory result. From the figure we observe that the reconstruction is numerically near-perfect when τ * = 0.5. The iterative procedure converges quickly and the reconstruction solution seems to be very satisfactory. Since the function a(y) changes drastically around zero, it is quite difficult to guarantee the convergence at this point. In fact, this form of a function has a larger fluctuation, which also shows that the drift function changes more violently in the real market. However, our algorithm still performs well, and the shape reconstruction of cusp is very good. In the experiment, we find that 400 iterations are not as good as 380 iterations. Since the observation data contains error, which comes from the calculation of forward problem, to obtain stable numerical results, we will cease the iteration at some suitable time. where L, h, τ are the same as in the first experiment. The numerical results for exact input data can be seen in Fig. 4.
To investigate the stability of the numerical solution, we employ the following noisy data: W δ (y, 1) = W (y, 1) 1 + δ × random(y) with δ = 0.001 and δ = 0.01. The reconstructed results are shown in Figs. 5 and 6. From these two figures we can see that the reconstruction of a(y) with the noisy data is also satisfactory. Like in Fig. 1, for all times τ * , the reconstructed a(y) is near-perfect on the interval y ∈ (0, 2).
From Fig. 5 we observe that the reconstruction of the 0.1% relative random noise data is almost identical to the reconstruction using noiseless data shown in Fig. 4. From Fig. 6 we observe that the reconstruction of the 1% relative random noise data has a noticeable gap when compared to the reconstruction of the noiseless data. We can see that the reconstruction of the 1% has an upward trend. Nonetheless, even in this case, the form of the reconstruction of the noiseless data is maintained, and we think that changes depend on the size of the error (i.e., noise).
From the above results we conclude that our numerical algorithm for reconstructing trend coefficients is indeed stable for data containing 0.1% and 1% relative random noise data.
Remark 6.2 The parameters of numerical examples are taken as r = 0.5 and σ 0 = 1. However, in the real market, the range of volatility is usually [0.2, 0.8], and the risk-free interest rate r hardly reaches 0.5. So, to better adapt to the real market, we consider the more suitable parameters r = 0.09 and σ 0 = 0.5. Moreover, to observe the effect of τ * on the reconstruction process, we also consider different values of τ * . The numerical results are shown in Fig. 7. We can see that our algorithm still performs well for different parameters r, σ , and τ * , and the simulated drift rate is in good agreement with the real one.

Concluding remarks
In this paper, we discuss an inverse problem of reconstructing the drift rate coefficient of stock index options using market observation data. Considering the boundlessness and nonhomogeneity of the original model, we use the artificial boundary method and homogenization technique to transform the original problem into a terminal control problem of homogeneous initial-boundary value equation on a bounded domain. We strictly prove the well-posedness of the minimizer of the control problem and give an iterative calculation scheme. Numerical results show that our algorithm is fast and robust.
This paper focuses on the reconstruction of the drift coefficient. To simplify the problem, we assume that the volatility is constant. This assumption usually does not hold in practice. Therefore an interesting question is what conditions need to be imposed to reconstruct the volatility and drift rate simultaneously, which is also the future work for our research.