An iterative algorithm for evaluating approximations to the optimal exercise boundary for a nonlinear Black–Scholes equation

In the past years, the Black–Scholes equation for pricing derivatives has attracted a lot of attention from both theoretical as well as practical point of view. Recall that a European Call (Put) option is the right but not obligation to purchase (sell) an underlying asset at the expiration price $E$ at the expiration time $T$. In an idealized financial market the price of an option can be computed from the well-known Black–Scholes equation derived by Black and Scholes in [4], and, independently by Merton (see also [21], Dewynne et al. [7], Hull [16]). Assuming that the underlying asset follows a geometric Brownian motion one can derive a governing partial differential equation for the price of an option. We remind ourselves that the equation for option’s price $V(S,t)$ is the following parabolic PDE:

$$\partial_{t}V+(r-q)S\partial_{S}V+{\scriptstyle\frac{1}{2}}\sigma^{2}S^{2}\partial_{S}^{2}V-rV=0$$

(1)

where $\sigma$ is the volatility of the underlying asset price process, $r>0$ is the interest rate of a zero-coupon bond, $q\geq 0$ is the dividend yield rate. A solution $V=V(S,t)$ represents the price of an option at time $t\in[0,T]$ if the price of an underlying asset is $S>0$. In this paper we shall focus our attention to the case when the diffusion coefficient $\sigma^{2}$ may depend on the time $T-t$ to expiry, the asset price $S$ and the second derivative $\partial^{2}_{S}V$ of the option price (referred to as $\Gamma$), i.e.

$$\sigma=\sigma(S^{2}\partial^{2}_{S}V,S,T-t)\,.$$

(2)

A motivation for studying the nonlinear Black–Scholes equation (1) with a volatility $\sigma$ having a general form (2) arises from option pricing models taking into account nontrivial transaction costs, market feedbacks and/or risk from a volatile (unprotected) portfolio. Recall that the linear Black–Scholes equation with $\sigma$ constant has been derived under several restrictive assumptions like e.g. frictionless, liquid, complete markets, etc. We also recall that the linear Black–Scholes equation provides a perfectly replicated hedging portfolio. In recent years, some of these assumptions have been relaxed in order to model, for instance, the presence of transaction costs (see e.g. Leland [22], Hoggard et al. [17], Avellaneda and Paras [2]), feedback and illiquid market effects due to large traders choosing given stock-trading strategies (Frey and Patie [12], Frey and Stremme [13], During et al.[8], Schönbucher and Wilmott [29]), imperfect replication and investor’s preferences (Barles and Soner [5]), risk from unprotected portfolio (Kratka [20], Jandačka and Ševčovič [18]). One of the first nonlinear models is the so-called Leland model (c.f. [22]) for pricing Call and Put options under the presence of transaction costs. It has been generalized for more complex options by Hoggard, Whaley and Wilmott in [17]. In this model the volatility $\sigma$ is given by $\sigma^{2}(S^{2}\partial^{2}_{S}V,S,\tau)=\hat{\sigma}^{2}(1+\hbox{Le}\,\hbox{sgn}(\partial^{2}_{S}V))$ where $\hat{\sigma}>0$ is a constant historical volatility of the underlying asset price process and $\hbox{Le}>0$ is the so-called Leland constant given by $\hbox{Le}=\sqrt{2/\pi}C/(\hat{\sigma}\sqrt{\Delta t})$ where $C>0$ is a constant round trip transaction cost per unit dollar of transaction in the assets market and $\Delta t>0$ is the time-lag between portfolio adjustments.

If transaction costs are taken into account perfect replication of the contingent claim is no longer possible and further restrictions are needed in the model. By assuming that investor’s preferences are characterized by an exponential utility function Barles and Soner (c.f. [5]) derived a nonlinear Black–Scholes equation with the volatility $\sigma$ given by

$$\sigma^{2}(S^{2}\partial^{2}_{S}V,S,\tau)=\hat{\sigma}^{2}\left(1+\Psi(a^{2}e^{r\tau}S^{2}\partial^{2}_{S}V)\right)$$

(3)

where $\Psi$ is a solution to the ODE: $\Psi^{\prime}(x)=(\Psi(x)+1)/(2\sqrt{x\Psi(x)}-x),\Psi(0)=0,$ and $a>0$ is a given constant representing risk aversion. Notice that $\Psi(x)=O(x^{\frac{1}{3}})$ for $x\to 0$ and $\Psi(x)=O(x)$ for $x\to\infty$.

Another popular model has been derived for the case when the asset dynamics takes into account the presence of feedback and illiquid market effects. Frey and Stremme (c.f. [13, 12]) introduced directly the asset price dynamics in the case when a large trader chooses a given stock-trading strategy (see also [29]). The diffusion coefficient $\sigma$ is again nonconstant and it can be expressed as:

$$\sigma^{2}(S^{2}\partial^{2}_{S}V,S,\tau)=\hat{\sigma}^{2}\left(1-\varrho\lambda(S)S\partial^{2}_{S}V\right)^{-2}$$

(4)

where $\hat{\sigma}^{2},\varrho>0$ are constants and $\lambda(S)$ is a strictly convex function, $\lambda(S)\geq 1$.

The last example of the Black–Scholes equation with a nonlinearly depending volatility is the so-called Risk Adjusted Pricing Methodology model proposed by Kratka in [20] and revisited by Jandačka and Ševčovič in [18]. In order to maintain (imperfect) replication of a portfolio by the delta hedge one has to make frequent portfolio adjustments leading to a substantial increase in transaction costs. On the other hand, rare portfolio adjustments may lead to an increase of the risk arising from a volatile (unprotected) portfolio. In the RAPM model the aim is to optimize the time-lag $\Delta t$ between consecutive portfolio adjustments. By choosing $\Delta t>0$ in such way that the sum of the rate of transaction costs and the rate of a risk from unprotected portfolio is minimal one can find the optimal time lag $\Delta t>0$ (see [18] for details). In the RAPM model, it turns out that the volatility is again nonconstant and it has the following form:

$$\sigma^{2}(S^{2}\partial^{2}_{S}V,S,\tau)=\hat{\sigma}^{2}\left(1+\mu(S\partial^{2}_{S}V)^{\frac{1}{3}}\right)\,.$$

(5)

Here $\hat{\sigma}^{2}>0$ is a constant and $\mu=3(C^{2}R/2\pi)^{\frac{1}{3}}$ where $C,R\geq 0$ are nonnegative constant representing the transaction cost measure and the risk premium measure, resp. (see [18] for details).

Notice that all the above mentioned nonlinear models are consistent with the original Black–Scholes equation in the case the additional model parameters (e.g. Le, $a$, $\varrho$, $\mu$) are vanishing. If plain Call or Put vanilla options are concerned then the function $V(S,t)$ is convex in $S$ variable and therefore each of the above mentioned models has a diffusion coefficient strictly larger than $\hat{\sigma}^{2}$ leading to a larger values of computed option prices. They can be therefore identified with higher Ask option prices, i.e. offers to sell an option. Furthermore, these models have been considered and analyzed mostly for European options, i.e. options can be exercised only at the maturity $t=T$. On the other hand, American options are more common in financial markets as they allow for exercising of an option anytime before the expiry $T$. In the case of an American Call option a solution to equation (1) is defined on a time dependent domain $0<t<T,\ 0<S<S_{f}(t)$. It is subject to the boundary conditions

$$V(0,t)=0\,,\ \ V(S_{f}(t),t)=S_{f}(t)-E\,,\ \ \partial_{S}V(S_{f}(t),t)=1\,,$$

(6)

and terminal pay-off condition at expiry $t=T$

$$V(S,T)=\max(S-E,0)$$

(7)

where $E>0$ is a strike price (c.f. [7, 21]). One of important problems in this field is the analysis of the early exercise boundary $S_{f}(t)$ and the optimal stopping time (an inverse function to $S_{f}(t)$) for American Call options on stocks paying a continuous dividend $q>0$. However, an exact analytical expression for the free boundary profile is not even known for the case when the volatility $\sigma$ is constant. Many authors have investigated various approximation models leading to approximate expressions for valuing American Call and Put options: analytic approximations (Barone–Adesi and Whaley [3], Kuske and Keller [19], Dewynne et al. [7], Geske et al. [14, 15], MacMillan [23], Mynemi [26]); methods of reduction to a nonlinear integral equation (Alobaidi [1], Kwok [21], Mallier et al. [24, 25], Ševčovič [28], Stamicar et al. [30]). Concerning numerical methods for solving free boundary problem we refer to the book by Kwok [21] and recent papers by Ehrhardt and Mickens [9] and Zhao et al. [31].

The main goal of this paper is to propose a new iterative numerical algorithm for solving the free boundary problem for an American Call option in the case the volatility $\sigma$ may depend on the option and asset values as well as on the time $T-t$ to expiry. The key idea of the proposed algorithm consists in transformation of the free boundary problem into a semilinear parabolic equation defined on a fixed spatial domain coupled with a nonlocal algebraic constraint equation for the free boundary position. Since the resulting parabolic equation contains a strong convective term we make use of the operator splitting method in order to overcome numerical difficulties. Full space-time discretization of the problem leads to a system of semi-linear algebraic equations that can be solved by an iterative procedure at each time level.

The paper has the following plan: in the next section we transform the free boundary problem into a system consisting of a nonlinear parabolic equation defined on a fixed domain with time depending coefficients and an algebraic constraint equation for the free boundary position. In section 3 we present a numerical discretization scheme based on the idea of operator splitting. In the last section 4 we present several numerical results for nonlinear Black–Scholes equations with volatility functions $\sigma$ defined in (3) and 5). We make a comparison to well-known methods in the case the volatility $\sigma$ is constant. Finally, we discuss dependence of the free boundary position with respect to various parameters entering expressions (3) and (5).

The main goal of this section is to perform a fixed domain transformation (referred to as Landau’s transformation) of the free boundary problem for the nonlinear Black–Scholes equation (1) into a parabolic equation defined on a fixed spatial domain. For the sake of simplicity we will present a detailed derivation of an equation only for the case of an American Call option. Derivation of the corresponding equation for the American Put option is similar.

Let us consider the following change of variables:

$$\tau=T-t,\quad x=\ln\left(\varrho(\tau)/S\right)\ \ \hbox{where}\ \ \varrho(\tau)=S_{f}(T-\tau).$$

Then $\tau\in(0,T)$ and $x\in(0,\infty)$ iff $S\in(0,S_{f}(t))$. The value $x=0$ corresponds to the free boundary position $S=S_{f}(t)$ whereas $x\approx+\infty$ corresponds to the default value $S=0$ of the underlying asset. Following Stamicar et al. [30] and Ševčovič [28] we construct the so-called synthetic portfolio function $\Pi=\Pi(x,\tau)$ defined as follows:

$$\Pi(x,\tau)=V(S,t)-S\partial_{S}V(S,t)\,.$$

(8)

It corresponds to a synthetic portfolio consisting of one long positioned option and $\Delta=\partial_{S}V$ underlying short stocks. Clearly, we have

$$\partial_{x}\Pi=S^{2}\partial_{S}^{2}V,\ \ \partial_{\tau}\Pi+\frac{\dot{\varrho}}{\varrho}\partial_{x}\Pi=-\partial_{t}\left(V-S\partial_{S}V\right)$$

where we have denoted $\dot{\varrho}=d\varrho/d\tau$. Assuming sufficient smoothness of a solution $V=V(S,t)$ to (1) we can deduce from (1) a parabolic equation for the synthetic portfolio function $\Pi=\Pi(x,\tau)$

$$\partial_{\tau}\Pi+(b(\tau)-{\scriptstyle\frac{1}{2}}\sigma^{2})\partial_{x}\Pi-{\scriptstyle\frac{1}{2}}\partial_{x}\left(\sigma^{2}\partial_{x}\Pi\right)+r\Pi=0$$

defined on a fixed domain $x\in R,\,t\in(0,T),$ with a time-dependent coefficient

$$b(\tau)={\dot{\varrho}(\tau)\over\varrho(\tau)}+r-q$$

(9)

and a diffusion coefficient given by: $\sigma^{2}=\sigma^{2}(\partial_{x}\Pi(x,\tau),\varrho(\tau)e^{-x},\tau)$ depending on $\tau,x$ and the gradient $\partial_{x}\Pi$ of a solution $\Pi$. Now the boundary conditions $V(0,t)=0,V(S_{f}(t),t)=S_{f}(t)-E$ and $\partial_{S}V(S_{f}(t),t)=1$ imply

$$\Pi(0,\tau)=-E,\quad\Pi(+\infty,\tau)=0\,,\quad 0<\tau<T\,,$$

(10)

and, from the terminal pay-off diagram for $V(S,T)$, we deduce

$$\Pi(x,0)=\left\{\begin{array}[]{lll}-E&\hbox{for}\ x<\ln\left({\varrho(0)\over E}\right)\hfil\\ \ \ 0&\hbox{otherwise.}\hfil\end{array}\right.$$

(11)

In order to close up the system of equations that determines the value of a synthetic portfolio $\Pi$ we have to construct an equation for the free boundary position $\varrho(\tau)$. Indeed, both the coefficient $b$ as well as the initial condition $\Pi(x,0)$ depend on the function $\varrho(\tau)$. Similarly as in the case of a constant volatility $\sigma$ (see [28, 30]) we proceed as follows: since $S_{f}(t)-E=V(S_{f}(t),t)$ and $\partial_{S}V(S_{f}(t),t)=1$ we have ${\scriptstyle\frac{d}{dt}}S_{f}(t)=\partial_{S}V(S_{f}(t),t){\scriptstyle\frac{d}{dt}}S_{f}(t)+\partial_{t}V(S_{f}(t),t)$ and so $\partial_{t}V(S,t)=0$ along the free boundary $S=S_{f}(t)$. Moreover, assuming $\partial_{x}\Pi$ is continuous up to the boundary $x=0$ we obtain $S^{2}\partial^{2}_{S}V(S,t)\to\partial_{x}\Pi(0,\tau)$ and $S\partial_{S}V(S,t)\to\varrho(\tau)$ as $S\to S_{f}(t)^{-}$. Now, by taking the limit $S\to S_{f}(t)^{-}$ in the Black–Scholes equation (1) we obtain $(r-q)\varrho(\tau)+{\scriptstyle\frac{1}{2}}\sigma^{2}\partial\Pi(0,\tau)-r(\varrho(\tau)-E)=0$. Therefore

$$\varrho(\tau)={rE\over q}+\frac{1}{2q}\sigma^{2}(\partial_{x}\Pi(0,\tau),\varrho(\tau),\tau)\partial_{x}\Pi(0,\tau)$$

for $0<\tau\leq T$. Recall that in the linear case when $\sigma>0$ is constant the initial position of the interface $\varrho(0)$ is given by $\varrho(0)=rE/q$ if $r\geq q>0$ and $\varrho(0)=E$ if $0<r<q$ (see Dewynne et al. [7] or Ševčovič [28]). We also recall that the value of $\varrho(0)$ in the case $r\geq q>0$ can derived easily from the smoothness assumption made on $\partial_{x}\Pi$ at the origin $x=0,\tau=0$. Indeed, continuity of $\partial_{x}\Pi$ at the origin $(0,0)$ implies $\lim_{\tau\to 0^{+}}\partial_{x}\Pi(0,\tau)=\partial_{x}\Pi(0,0)=\lim_{x\to 0^{+}}\partial_{x}\Pi(x,0)=0$ because $\Pi(x,0)=-E$ for $x$ close to $0^{+}$. From the above equation for $\varrho(\tau)$ we deduce $\varrho(0)={rE\over q}$ by taking the limit $\tau\to 0^{+}$. Henceforth, we restrict our attention to the case when the interest rate is greater then the dividend yield rate, i.e.

$$0<q\leq r$$

(12)

leading to the initial position of the free boundary $\varrho(0)=rE/q$. Putting all the above equations together we end up with a closed system of equations for $\Pi=\Pi(x,\tau)$ and $\varrho=\varrho(\tau)$

	$\displaystyle\partial_{\tau}\Pi+(b(\tau)-{\scriptstyle\frac{1}{2}}\sigma^{2})\partial_{x}\Pi-{\scriptstyle\frac{1}{2}}\partial_{x}\left(\sigma^{2}\partial_{x}\Pi\right)+r\Pi=0\,,$
	$\displaystyle\Pi(0,\tau)=-E\,,\ \ \Pi(+\infty,\tau)=0\,,\ x>0\,,\tau\in(0,T)\,,$
	$\displaystyle\Pi(x,0)=\left\{\begin{array}[]{lll}-E&\quad\hbox{for}\ x<\ln(r/q)\hfil\\ \ \ \ 0&\quad\hbox{otherwise\,,}\hfil\end{array}\right.$		(15)

where $\sigma=\sigma(\partial_{x}\Pi(x,\tau),\varrho(\tau)e^{-x},\tau)\,,\ b(\tau)={\dot{\varrho}(\tau)\over\varrho(\tau)}+r-q$ and the free boundary position $\varrho(\tau)=S_{f}(T-\tau)$ satisfies an implicit algebraic equation

$$\varrho(\tau)={rE\over q}+{\sigma^{2}(\partial_{x}\Pi(0,\tau),\varrho(\tau),\tau)\over 2q}\partial_{x}\Pi(0,\tau)\,,\quad\hbox{with }\ \varrho(0)={rE\over q}$$

(16)

where $\tau\in(0,T)$. Notice that, in order to guarantee parabolicity of equation (15) we have to assume that the function $p\mapsto\sigma^{2}(p,\varrho(\tau)e^{-x},\tau)p$ is strictly increasing. More precisely, we shall assume that there exists a positive constant $\gamma>0$ such that

$$\sigma^{2}(p,\xi,\tau)+p\partial_{p}\sigma^{2}(p,\xi,\tau)\geq\gamma>0$$

(17)

for any $\xi>0,\tau\in(0,T)$ and $p\in R$.

In this section we derive a full space-time discretization scheme for a numerical approximation of the problem (15), (16). Recall that in the case of a constant volatility there are, in principle, two ways how to solve numerically the free boundary problem for the value of an American Call resp. Put option and the position of the early exercise boundary. The first class of algorithms is based on reformulation of the problem in terms of a variational inequality (see Kwok [21] and references therein). The variational inequality can be then solved numerically by the so-called Projected Super Over Relaxation method (PSOR for short). An advantage of this method is that it gives us immediately the value of a solution. A disadvantage is that one has to solve large systems of linear equations iteratively taking into account the obstacle for a solution, and, secondly, the free boundary position should be deduced from the solution a posteriori. Moreover, the PSOR method is not directly applicable for solving the problem (1)-(6) when the diffusion coefficient $\sigma$ may depend on the second derivative of a solution itself. The second class of methods is based on derivation of a nonlinear integral equation for the position of the free boundary without the need of knowing the option price itself (see e.g. Kuske-Keller [19], Mallier and Alobaidi [1, 24, 25], Ševčovič et al. [28, 30]. In this approach an advantage is that only a single equation for the free boundary has to be solved; a disadvantage is that the method is based on integral transformation techniques and therefore the assumption $\sigma$ is constant is crucial.

In our approach we make an attempt to take advantages of both above mentioned methods. As it was mentioned in the previous section we are going to solve the system of nonlinear equations (15) with constraint (16). Because the volatility $\sigma$ may be nonconstant we cannot use integral transformation techniques in order to derive a single integral equation for $\varrho(\tau)$. However, the form of the system (15), (16) allows for an efficient and fast numerical algorithm for computing of the early exercise boundary position $\varrho(\tau)=S_{f}(T-\tau)$.

The idea of the iterative numerical algorithm for solving the problem (15), (16) is rather simple: we use the backward Euler method of finite differences in order to discretize the parabolic equation (15) in time. In each time level we find a new approximation of a solution pair $(\Pi,\varrho)$. First we determine a new position of $\varrho$ from the algebraic equation (16). We remind ourselves that (even in the case $\sigma$ is constant) the free boundary function $\varrho(\tau)$ behaves like $rE/q+O(\tau^{1/2})$ for $\tau\to 0^{+}$ (see e.g. [7] or [28]) and so $b(\tau)=O(\tau^{-1/2})$. Hence the convective term $b(\tau)\partial_{x}\Pi$ becomes a dominant part of equation (15) for small values of $\tau$. In order to overcome this difficulty we use the operator splitting method for successive solving of the convective and diffusion parts of equation (15). Since the diffusion coefficient depends on the solution $\Pi$ itself we make several micro-iterates to find a solution of a system of nonlinear algebraic equations.

Now we present our algorithm in more details. We restrict the spatial domain $x\in(0,\infty)$ to a finite interval of values $x\in(0,L)$ where $L>0$ is sufficiently large. For practical purposes one can take $L\approx 3$ as it corresponds to the interval $S\in(S_{f}(t)e^{-L},S_{f}(t))$ in the original asset price variable $S$. The value $S_{f}(t)e^{-L}$ is then a good approximation for the default value $S=0$ if $L\approx 3$. Let us denote by $k>0$ the time step, $k=T/m$, and, by $h>0$ the spatial step, $h=L/n$ where $m,n\in N$ stand for the number of time and space discretization steps, resp. We denote by $\Pi_{i}^{j}$ an approximation of $\Pi(x_{i},\tau_{j})$, $\varrho^{j}\approx\varrho(\tau_{j})$, $b^{j}\approx b(\tau_{j})$ where $x_{i}=ih,\tau_{j}=jk$. We approximate the value of the volatility $\sigma$ at the node $(x_{i},\tau_{j})$ by finite difference as follows:

$$\sigma_{i}^{j}=\sigma_{i}^{j}(\varrho^{j},\Pi^{j})=\sigma((\Pi_{i+1}^{j}-\Pi_{i}^{j})/h,\varrho^{j}e^{-x_{i}},\tau_{j})\,.$$

Then for the backward in time Euler finite difference approximation of equation (15) we have

$$\frac{\Pi^{j}-\Pi^{j-1}}{k}+\left(b^{j}-\frac{1}{2}(\sigma^{j})^{2}\right)\partial_{x}\Pi^{j}-\frac{1}{2}\partial_{x}\left((\sigma^{j})^{2}\partial_{x}\Pi^{j}\right)+r\Pi^{j}=0$$

(20)

and the solution $\Pi^{j}=\Pi^{j}(x)$ is subject to Dirichlet boundary conditions at $x=0$ and $x=L$. We set $\Pi^{0}(x)=\Pi(x,0)$. Now we decompose the above problem into two parts - a convection part and a diffusive part by introducing an auxiliary intermediate step $\Pi^{j-{\scriptstyle\frac{1}{2}}}$:

(Convective part)

$$\frac{\Pi^{j-{\scriptstyle\frac{1}{2}}}-\Pi^{j-1}}{k}+b^{j}\partial_{x}\Pi^{j-{\scriptstyle\frac{1}{2}}}=0\,,$$

(21)

(Diffusive part)

$$\frac{\Pi^{j}-\Pi^{j-{\scriptstyle\frac{1}{2}}}}{k}-(\sigma^{j})^{2}\partial_{x}\Pi^{j}-\frac{1}{2}\partial_{x}\left((\sigma^{j})^{2}\partial_{x}\Pi^{j}\right)+r\Pi^{j}=0\,.$$

(22)

The idea of the operator splitting technique now consists in comparison the sum of solutions to convective and diffusion part to a solution of (20). Indeed, if $\partial_{x}\Pi^{j}\approx\partial_{x}\Pi^{j-{\scriptstyle\frac{1}{2}}}$ then it is reasonable to assume that $\Pi^{j}$ computed from the system (21)–(22) is a good approximation of the system (20).

The convective part can be approximated by an explicit solution to the transport equation:

$$\partial_{\tau}\tilde{\Pi}+b(\tau)\partial_{x}\tilde{\Pi}=0\qquad\hbox{for }x>0,\tau\in(\tau_{j-1},\tau_{j}]$$

(23)

subject to the boundary condition $\tilde{\Pi}(0,\tau)=-E$ and initial condition $\tilde{\Pi}(x,\tau_{j-1})=\Pi^{j-1}(x)$. Since the free boundary $\varrho(\tau)=S_{f}(T-\tau)$ must be an increasing function in $\tau$ and we have assumed $0<q\leq r$ we have $b(\tau)=\dot{\varrho}(\tau)/\varrho(\tau)+r-q>0$ and so the in-flowing boundary condition $\tilde{\Pi}(0,\tau)=-E$ is consistent with the transport equation. Let us denote by $B(\tau)$ the primitive function to $b(\tau)$, i.e. $B(\tau)=\ln\varrho(\tau)+(r-q)\tau$. Equation (23) can be integrated to obtain its explicit solution:

$$\tilde{\Pi}(x,\tau)=\left\{\matrix{\Pi^{j-1}(x-B(\tau)+B(\tau_{j-1}))\hfill&\quad\hbox{if }x-B(\tau)+B(\tau_{j-1})>0\,,\hfill\cr-E\hfill&\quad\hbox{otherwise.}\hfill}\right.$$

(24)

Thus the spatial approximation $\Pi^{j-{\scriptstyle\frac{1}{2}}}_{i}$ can be constructed from the formula

$$\Pi^{j-{\scriptstyle\frac{1}{2}}}_{i}=\left\{\matrix{\Pi^{j-1}(\xi_{i})\hfill&\quad\hbox{if }\xi_{i}=x_{i}-\ln\varrho_{j}+\ln\varrho_{j-1}-(r-q)k>0\,,\hfill\cr-E\hfill&\quad\hbox{otherwise,}\hfill}\right.$$

(25)

where a linear approximation between discrete values $\Pi^{j-1}_{i},i=0,1,...,n,$ is being used to compute the value $\Pi^{j-1}(x_{i}-\ln\varrho_{j}+\ln\varrho_{j-1}-(r-q)k)$.

The diffusive part can be solved numerically by means of finite differences. Using central finite difference for approximation of the derivative $\partial_{x}\Pi^{j}$ we obtain

$\displaystyle\frac{\Pi_{i}^{j}-\Pi_{i}^{j-{\scriptstyle\frac{1}{2}}}}{k}+r\Pi_{i}^{j}-(\sigma^{j}_{i})^{2}\frac{\Pi_{i+1}^{j}-\Pi_{i-1}^{j}}{2h}-\frac{1}{2h}\left((\sigma_{i}^{j})^{2}\frac{\Pi_{i+1}^{j}-\Pi_{i}^{j}}{h}-(\sigma_{i-1}^{j})^{2}\frac{\Pi_{i}^{j}-\Pi_{i-1}^{j}}{h}\right)=0\,.$

Hence, the vector of discrete values $\Pi^{j}=\{\Pi_{i}^{j},i=1,2,...,n\}$ at the time level $j\in\{1,2,...,m\}$ satisfies the tridiagonal system of equations

$$\alpha_{i}^{j}\Pi_{i-1}^{j}+\beta_{i}^{j}\Pi_{i}^{j}+\gamma_{i}^{j}\Pi_{i+1}^{j}=\Pi_{i}^{j-{\scriptstyle\frac{1}{2}}}$$

(26)

for $i=1,2,...,n,$ where

	$\displaystyle\alpha_{i}^{j}$	$\displaystyle\equiv$	$\displaystyle\alpha_{i}^{j}(\varrho^{j},\Pi^{j})=-\frac{k}{2h^{2}}(\sigma_{i-1}^{j})^{2}+\frac{k}{2h}\frac{(\sigma_{i}^{j})^{2}}{2}$
	$\displaystyle\gamma_{i}^{j}$	$\displaystyle\equiv$	$\displaystyle\gamma_{i}^{j}(\varrho^{j},\Pi^{j})=-\frac{k}{2h^{2}}(\sigma_{i}^{j})^{2}-\frac{k}{2h}\frac{(\sigma_{i}^{j})^{2}}{2}$		(27)
	$\displaystyle\beta_{i}^{j}$	$\displaystyle\equiv$	$\displaystyle\beta_{i}^{j}(\varrho^{j},\Pi^{j})=1+rk-(\alpha_{i}^{j}+\gamma_{i}^{j})\,.$

The initial and boundary conditions at $\tau=0$ and $x=0,L,$ resp., can be approximated as follows:

$$\Pi_{i}^{0}=\left\{\begin{array}[]{lll}-E&\ \ \hbox{for}\ x_{i}<\ln\left({r/q}\right)\hfil\\ \ \ \ 0&\ \ \hbox{for}\ x_{i}\geq\ln\left({r/q}\right)\hfil\end{array}\right.$$

for $i=0,1,...,n,$ and $\Pi_{0}^{j}=-E,\quad\Pi_{n}^{j}=0$.

Next we proceed by approximation of equation (16) which introduces a nonlinear constraint condition between the early exercise boundary function $\varrho(\tau)$ and the trace of the solution $\Pi$ at the boundary $x=0$ ($S=S_{f}(t)$ in the original variable). Taking a finite difference approximation of $\partial_{x}\Pi$ at the origin $x=0$ we obtain

$$\varrho^{j}=\frac{rE}{q}+\frac{1}{2q}\sigma^{2}\left((\Pi_{1}^{j}-\Pi_{0}^{j})/h,\varrho^{j},\tau\right)\frac{\Pi_{1}^{j}-\Pi_{0}^{j}}{h}\,.$$

(28)

Now, equations (25), (26) and (28) can be written in an abstract form as a system of nonlinear equations:

	$\displaystyle\varrho^{j}$		$\displaystyle={\mathcal{F}}(\Pi^{j},\varrho^{j})$
	$\displaystyle\Pi^{j-{\scriptstyle\frac{1}{2}}}$		$\displaystyle={\mathcal{T}}(\Pi^{j},\varrho^{j})$		(29)
	$\displaystyle{\mathcal{A}}(\Pi^{j},\varrho^{j})\Pi^{j}$		$\displaystyle=\Pi^{j-{\scriptstyle\frac{1}{2}}}$

where ${\mathcal{F}}(\Pi^{j},\varrho^{j})$ is the right-hand side of the algebraic equation (28), ${\mathcal{T}}(\Pi^{j},\varrho^{j})$ is the transport equation solver given by the right-hand side of (25) and ${\mathcal{A}}={\mathcal{A}}(\Pi^{j},\varrho^{j})$ is a tridiagonal matrix with coefficients given by (3). The system (29) can be approximately solved by means of successive iterates procedure. We define, for $j\geq 1,$ $\Pi^{j,0}=\Pi^{j-1},\varrho^{j,0}=\varrho^{j-1}$. Then the $(p+1)$-th approximation of $\Pi^{j}$ and $\varrho^{j}$ is obtained as a solution to the system:

	$\displaystyle\varrho^{j,p+1}$		$\displaystyle={\mathcal{F}}(\Pi^{j,p},\varrho^{j,p})$
	$\displaystyle\Pi^{j-{\scriptstyle\frac{1}{2}},p+1}$		$\displaystyle={\mathcal{T}}(\Pi^{j,p},\varrho^{j,p+1})$		(30)
	$\displaystyle{\mathcal{A}}(\Pi^{j,p},\varrho^{j,p+1})\Pi^{j,p+1}$		$\displaystyle=\Pi^{j-{\scriptstyle\frac{1}{2}},p+1}\,.$

Notice that the last equation is a linear tridiagonal equation for the vector $\Pi^{j,p+1}$ whereas $\varrho^{j,p+1}$ and $\Pi^{j-{\scriptstyle\frac{1}{2}},p+1}$ can be directly computed from (28) and (25), resp. Now, if the sequence of approximate solutions $\{(\Pi^{j,p},\varrho^{j,p})\}_{p=1}^{\infty}$ converges to some limiting value $(\Pi^{j,\infty},\varrho^{j,\infty})$ then this limit is a solution to a nonlinear system of equations (29) at the time level $j$ and we can proceed by computing the approximate solution the next time level $j+1$.

In this section we focus on numerical experiments based on the iterative scheme described in the previous section. The main purpose is to compute the free boundary profile $S_{f}(t)=\varrho(T-t)$ for different (non)linear Black–Scholes models and for various model parameters. We used $n=750$ spatial points and $m=225000$ time discretization steps. In average we needed $p\leq 6$ micro-iterates (30) in order to solve the nonlinear system (29) with the precision $10^{-7}$. A solution $(\Pi,\varrho)$ has been computed by our iterative algorithm for the following basic model parameters: $E=10,T=1,r=0.1,q=0.05,$ and $\hat{\sigma}=0.2$.