In a complete financial market without transaction costs, the Black–Scholes (B–S) no-arbitrage argument provides a rational option pricing formula and a hedging portfolio that replicates the contingent claim. Under the transaction costs, the continuous trading required by the hedging portfolio is prohibitively expensive, [1]. Several alternatives lead to option prices that are equal to Black–Scholes price but with an adjusted volatility. In 1992, Boyle and Vorst [2], derived from a binomial model an option price taking into account transaction costs and that is equal to a B–S price but with a modified volatility of the form Turn MathJax on
Here, μμ is the proportional transaction cost, ΔtΔt the transaction period, and σ0σ0 is the original volatility constant. Leland [3] computed View the MathML sourcec=(2π)1/2. Kusuoka [4] then showed that the “optimal” c depends on the risk structure of the market. Paras and Avellaneda [5] derived the modified volatility
View the MathML sourceσ=σ0(1+Asign(VSS))1/2,
Turn MathJax on
from a binomial model using the algorithm of Bensaid et al. [6]. Whalley and Wilmott [7] using an asymptotic analysis based on [8] propose the same adjusted volatility. A comparison of the exact hedging strategy of [8] and the asymptotic hedging strategy of [7] has been studied in [9]. Here, VV is the option price, S the price of the underlying asset, and VSSVSS denotes the second derivative of VV with respect to SS (the “Gamma”). In particular, the option price does not need to be convex.
Kratka in [10] and Jandačka and Ševčovič in [11] propose a correction of volatility of the form
View the MathML sourceσ2=σ02(1+μ(SVSS)13),
Turn MathJax on
where View the MathML sourceμ=3(C2R/2π)13 and C,RC,R are nonnegative constants representing the transaction cost measure and the risk premium measure, respectively.
A more complex model has been proposed by Barles and Soner [1], assuming that investor’s preferences are characterized by an exponential utility function. In their model the nonlinear volatility reads
equation(1.1)
View the MathML sourceσ2=σ02(1+Ψ[exp(r(T−t)a2S2VSS)]),
Turn MathJax on
where rr is the risk-free interest rate, TT the maturity, and View the MathML sourcea=μγN, with risk aversion factor γγ and the number NN of options to be sold. The function ΨΨ is the solution of the nonlinear initial value problem
equation(1.2)
View the MathML sourceΨ′(A)=Ψ(A)+12AΨ(A)−A,A≠0,Ψ(0)=0.
Turn MathJax on
In the mathematical literature, only a few results can be found on the numerical discretization of B–S equation, mainly for linear B–S equations. The numerical approaches vary from finite element discretizations [12] and [13], finite-difference approximations [14], [15] and [16]. The numerical discretization of the B–S equations with the nonlinear volatility (1.2) has been performed using explicit finite-difference schemes [1]. However, explicit schemes have the disadvantage that restrictive conditions on the discretization parameters (for instance, the ratio of the time and the space step) are needed in order to obtain stable, convergent schemes [17]. Moreover, the order of convergence is only one in time and two in space. [18] combines high-order compact difference schemes derived by [19] and techniques to construct numerical solutions with frozen values of the nonlinear coefficient of the nonlinear B–S equation to make the formulation linear.
In this paper we use a semidiscretization technique by using fourth-order difference approximations of the partial derivatives VSVS and VSSVSS arising in the nonlinear B–S equation
equation(1.3)
View the MathML sourceVt+12σ(VSS)2S2VSS+rSVS−rV=0.
Turn MathJax on
Then we achieve an ordinary system of nonlinear ordinary differential equations with respect to the time, that is solved numerically. Apart form (1.3), in the Barles–Soner model one has the terminal condition
equation(1.4)
View the MathML sourceV(S,T)=max(0,S−E),S>0,
Turn MathJax on
and the boundary conditions
equation(1.5)
View the MathML sourceV(0,t)=0,lims→∞V(S,t)S−Ee−r(T−t)=1.
Turn MathJax on
In order to compute the numerical solution, it is necessary to work in a bounded domain. Once this numerical domain has been chosen, the boundary conditions can be translated from the asymptotic condition (1.5), as it is done for instance in [1] or [18], or the boundary values must be found together with the solution and they are linked with the rest of the numerical solution in the interior of the numerical domain by using extrapolation techniques. This last approach is used in this paper in accordance with the used scheme.
Using the change of variable View the MathML sourceτ=T−t,U(S,τ)=V(S,t) Eq. (1.3) together with the initial condition (1.4) is transformed into
equation(1.6)
View the MathML sourceUτ−S22σ2USS−rSUS+rU=0,0<S<∞,0<τ≤T,
Turn MathJax on
equation(1.7)
U(S,0)=max(0,S−E).U(S,0)=max(0,S−E).
Turn MathJax on
This paper is organized as follows. Section 2 is addressed to the study of the properties of the volatility correction function ΨΨ after obtaining the implicit solution of (1.2). In Section 3, by using semidiscretization with respect to SS one gets a nonlinear system of ordinary differential equations with respect to the time, and then it is discretized using a forward explicit scheme. This approach allows us to study the stability and consistency of the nonlinear scheme in Sections 4 and 5 without using linearization strategies as it is done in [18]. Section 6 includes illustrative examples of European call option pricing where the computed numerical solution and their properties are checked.
If AA is a matrix in Rp×pRp×p and AtAt denotes its transposed matrix, we denote by ‖A‖‖A‖ the spectral norm of AA defined as, [20],
View the MathML source‖A‖=max{λ;λ eigenvalue of AtA}.
Turn MathJax on
If qq is an integer with |q|≤p−1|q|≤p−1, and AqAq is a band matrix in Rp×pRp×p such that Aq=(aij)Aq=(aij) with aij=0aij=0 everywhere outside of the diagonal j=i+qj=i+q, then it is easy to show that
equation(1.8)
View the MathML source
Dealing with reliable numerical computations of FD schemes, the consistency of the difference-scheme with the equation
is a necessary requirement because this means that the exact theoretical solution of the partial differential equation
approximates well to the exact solution of the difference equation as the stepsizes tend to zero, [23]. The strategy developed
by the authors in [1, p. 383], of using a very small time step near the maturity is an advisable decision but by no means a
guarantee that numerical results are reliable.
Let us represent Eq. (1.6) by L.U/ D 0, and let F .unj
/ D 0 represent the approximating difference equation defined by
(3.19) with exact solution u. In accordance with [23, p.100] where Un
j denotes the theoretical solution of (1.6) evaluated at the point .Sj; nk/, i.e., Un
j D U.Sj; nk/. If in (5.1) one has that
T n
j .U/ D O.hp/ C O.kq/, then we say that the FD scheme is consistent of order .p; q/.
Taking into account (3.19), (3.2) and (3.4) for the internal points we have
