Pricing American Options with a Non-Constant Penalty Parameter

Abstract

As the American early exercise results in a free boundary problem, in this article we add a penalty term to obtain a partial differential equation, and we also focus on an improved definition of the penalty term for American options. We replace the constant penalty parameter with a time-dependent function. The novelty and advantage of our approach consists in introducing a bounded, time-dependent penalty function, enabling us to construct an efficient, stable, and adaptive numerical approximation scheme, while in contrast, the existing standard approach to the penalisation of the American put option-free boundary problem involves a constant penalty parameter. To gain insight into the accuracy of our proposed extension, we compare the solution of the extension to standard reference solutions from the literature. This illustrates the improvement of using a penalty function instead of a penalising constant.

1. Introduction

As American options give the holder the right to exercise the option before and at maturity, this leads to free boundary value problems which have to be solved numerically. Several schemes for solving American option problems have been proposed, such as the projected SOR scheme (Cryer 1971; Wilmott et al. 1993), the binomial method, front-fixing schemes (Nielsen et al. 2002), the power penalty method (Wang et al. 2006), and Monte Carlo simulation techniques. These schemes compute the free boundary value implicitly. Other researchers focused on an explicit representation of the free boundary value (Lauko and Ševčovič 2010; Stamicar et al. 1999; Zhu 2006).

Another approach to solving American option problems is to add a penalty term to the problem (Günther and Jüngel 2010; Nielsen et al. 2002; Wang et al. 2006). A penalty term forces the problem to fulfill the free boundary constraint asymptotically. If the free boundary constraint is fulfilled, the penalty term is zero, otherwise it penalizes the problem with a factor. Until now, the penalty term included a penalisation constant being roughly estimated by an educated guess or an optimization. The recent approach of the power penalty function (Wang et al. 2006) improves the convergence and accuracy by introducing a non-constant penalization function of a given power form. Here, we further improve the performance of the penalization by replacing the penalty constant with a general penalty function, allowing for an adaptive penalization to be determined in a consistent way while solving the partial differential equation (PDE). An explicit formulation for the derivation of the parameters of the penalty function is part of our further research.

Once we add a penalty term to the American option problem, the problem reduces to a PDE on a fixed domain and we can apply standard numerical methods, such as finite difference methods or finite element methods. In front-fixing schemes, a change of variables is used to obtain a PDE on a fixed domain. Through the change of variables, the free boundary value is tracked implicitly (Han and Wu 2003; Schwartz 1977). In our approach, the domain is fixed by the initial guess of the free boundary value obtained by the intersection point of the payoff and the exact solution of the Black–Scholes equation. The penalisation of a time-step depends on the solution itself in comparison with the payoff (Günther and Jüngel 2010; Nielsen et al. 2002). Our presented penalty term is bounded by an initial guess of the free boundary value, and is thus independent from the solution of the last time-step. The usual approach to penalisation of the American put option-free boundary problem involves a small parameter, making the numerical analysis harder. The novelty and advantage of our approach consists in introducing a bounded penalty function, enabling us to construct an efficient and stable numerical approximation scheme.

The outline of this article is as follows. Section 2 reviews the mathematical modeling for an American option with and without a penalty term. As a benchmark, we chose the classical Black–Scholes equation and use a variable transformation to simplify our computations and add the penalty function. In Section 3 the model is discretized, and the numerical results of the different test cases are presented in Section 4. In Section 5 we conclude this work with a brief outlook.

2. Mathematical Modeling

American options are more expensive than European options, as American options give the holder the right to exercise the option before the maturity T. In the following, we will focus on American put options for clearness of the idea, but as it can be seen in the numerical results, all assumptions also hold analogously for the American call options. For pricing an American put option P with the Black–Scholes model, we seek a pair of functions $\left(P(S,t),S_f(t)\right)$ such that

$$\begin{array}{cc}\hfill {\mathcal{L}}_{BSM}[P]\equiv \frac{\partial P}{\partial t}& +\frac{{\sigma }^2}{2}{S}^2\frac{{\partial }^{2}P}{\partial S^2}+(r-q)S\frac{\partial P}{\partial S}-rP\le 0,0\le t\le T,\hfill \\ \hfill {(K-S)}^{+}& =P(S,t)\mathrm{for}S\le S_f(t),\hfill \\ \hfill {(K-S)}^{+}& <P(S,t)\mathrm{for}S>S_f(t),\hfill \end{array}$$

(1)

where K denotes the predefined strike price, r is the risk-free interest rate, q is the dividend rate, $\sigma$ is the volatility, S is the price of an asset, and $S_f(t)$ is the free boundary value at time t with $0\le t\le T$. In (1), we used the standard notation ${(f)}^{+}:=max(f,0)$. The differential operator appearing in the Black–Scholes PDE (1) is abbreviated by ${\mathcal{L}}_{BSM}$.

The terminal condition at maturity $t=T$ reads

$$P(S,T)={(K-S)}^{+},$$

(2)

and the "spatial" boundary conditions at $S=S_f(t)$, $S\to \infty$, are given by

$$P\left(S_f(t),t\right)={\left(K-S_f(t)\right)}^{+},\frac{\partial P}{\partial S}\left(S_f(t),t\right)=-1,\underset{S\to \infty }{lim}P(S,t)=0,0\le t\le T.$$

The American Put option problem can be split into two regions: the exercise region $0\le S\le S_f(t)$, and the holding region $S>S_f(t)$. In the holding region, the American Put option problem follows the Black–Scholes equation for the European Put option

$${\mathcal{L}}_{BSM}[P](S,t)=0,S>S_f(t).$$

For the exercise region, $P(S,t)={(K-S)}^{+}$ has to hold. By inserting this condition into ${\mathcal{L}}_{BSM}[P](S,t)$, we obtain

$${\mathcal{L}}_{BSM}[{(K-S)}^{+}]=-rS-r(K-S)=-rK<0,$$

(3)

fulfilling the condition (1). We obtain a case system for the American put option problem

$${\mathcal{L}}_{BSM}[P](S,t)=g(S,t)\left\{\begin{array}{cc}-rK,\hfill & 0<S\le S_f(t),\hfill \\ 0,\hfill & S>S_f(t).\hfill \end{array}\right.$$

(4)

To simplify the computational effort, we transform the problem (4) into the result that we obtain on the left-hand side of the heat equation. This transformation is well-known for the European options, since the heat equation has an analytical solution and is analysed very often. Since the right-hand side is in comparison to the European Put options unequal to zero, we sketch the transformation in more detail.

2.1. Transformation

For the rewriting of the problem (4), we use the standard transformations

$$x=ln\left(\frac{S}{K}\right),\tau =\frac{{\sigma }^2}{2}(T-t),v(x,\tau )=\frac{P(S,t)}{K}$$

(5)

and use the chain rule

$$-\frac{{\sigma }^2}{2}K{v}_{\tau }+\frac{{\sigma }^2}{2}{S}^2\frac{K}{S}(v_{xx}-v_x)+rS\frac{K}{S}{v}_x-rKv=\left\{\begin{array}{cc}-rK,\hfill & x\le x_f(\tau ),\hfill \\ 0,\hfill & x>x_f(\tau ),\hfill \end{array}\right.$$

where $x_f(\tau )=ln(S_f(\tau )/K)$ denotes the free boundary value for x. We simplify this equation and divide by $-\sigma /2$

$$v_{\tau }-v_{xx}+(1-k)v_x+kv=\left\{\begin{array}{cc}k,\hfill & x\le x_f(\tau ),\hfill \\ 0,\hfill & x>x_f(\tau ),\hfill \end{array}\right.$$

(6)

with $k=\frac{2r}{{\sigma }^2}$. In the last step, we use the transformation $u(x,\tau )=exp(\alpha x+\beta \tau )v(x,\tau )=e(x,\tau )v(x,\tau )$, where

$$\alpha =-\frac{1}{2}(k-1)\mathrm{and}\beta =-\frac{1}{4}{(k+1)}^2.$$

(7)

We obtain the transformed problem

$$u_{\tau }-u_{xx}=\tilde{g}(x,\tau )=exp(-\alpha x-\beta \tau )\left\{\begin{array}{cc}k\hfill & x\le x_f(\tau ),\hfill \\ 0,\hfill & x>x_f(\tau ).\hfill \end{array}\right.$$

(8)

with the transformed terminal condition (2)

$$f(x,\tau )=exp\left(-\alpha x-\beta \tau \right){\left(1-exp(x)\right)}^{+},$$

(9)

supplied with the initial and boundary conditions

$$u(x,0)=f(x,0),x\in \mathbb{R},\underset{x\to \pm \infty }{lim}\left(u(x,\tau )-f(x,\tau )\right)=0,0\le \tau \le T.$$

(10)

2.2. The Penalty Term

For the penalty term, we chose an affine function,

$$\delta (t)=at+b,$$

(11)

in a way that the penalisation function is given by $p(S,t)=\delta (t)g(S,t)$. This choice preserves the novelty of this approach, as the known penalty terms are neither bounded nor independent from the solution itself (Günther and Jüngel 2010; Nielsen et al. 2002). The penalised case system for (4) reads

$${\mathcal{L}}_{BSM}[P](S,t)=\delta (t)·g(S,t).$$

(12)

At $t=T$, the terminal condition (2) and (3) has to hold, and we set $p(t)=a(T-t)+1$. Similarly, we chose an affine function $\tilde{\delta }(\tau )=\tilde{a}\tau +b$ as a penalisation term for the penalisation function $\tilde{p}=\tilde{\delta }(\tau )\tilde{g}(x,\tau )$ for the transformed case system (8). Since the relation between $u(x,\tau )$ and $P(S,t)$ is given by

$$u(x,\tau )=\frac{1}{K}exp\left(-\alpha x-\beta \tau \right)P(S,t),$$

the same relation has to hold for the penalised right-hand side. We focus only on the exercise region since the right-hand side of the holding region is always zero. As with t, we step forward in time, and with $\tau$ backward, $\frac{1}{K}exp(-\alpha x-\beta \tau )p(S,t)=-\tilde{p}(x,\tau )$ can be assumed. We thus obtain

$$\begin{array}{ccc}\hfill \frac{-rK}{K}exp(-\alpha x-\beta \tau )(a(T-t)+1)& =& -(\tilde{a}\tau +\tilde{b})exp(-\alpha x-\beta \tau )\frac{2r}{{\sigma }^2}\hfill \\ \hfill a(T-t)+1& =& (\tilde{a}\frac{{\sigma }^2}{2}(T-t)+b)\frac{2}{{\sigma }^2}\hfill \\ \hfill \frac{{\sigma }^2}{2}& =& b.\hfill \end{array}$$

The transformed case system with the bounded penalisation function is given by

$$u_{\tau }-u_{xx}=\tilde{p}(x,\tau )=(\tilde{a}+\frac{{\sigma }^2}{2})\left\{\begin{array}{cc}kexp(-\alpha x-\beta \tau ),\hfill & x\le x_f(\tau )\hfill \\ 0,\hfill & x>x_f(\tau )\hfill \end{array}\right.,$$

(13)

supplied with the initial and boundary conditions

$$u(x,0)=f(x,0),x\in \mathbb{R},\underset{x\to \pm \infty }{lim}\left(u(x,\tau )-f(x,\tau )\right)=0,0\le \tau \le T.$$

(14)

Our approach requires an initial guess for the free boundary value. Since the price of an American option is larger than or equal to the price of the corresponding European option, we can use the intersection point ${\overline{x}}_f(\tau )$ between the payoff and the solution of the European Put option as the initial guess for the free boundary value $x_f(\tau )$, with the result that the bounded penalty term forces the PDE (12) to fulfill the conditions (1) asymptotically. As the approximation of the free boundary and the analysis for the approximation is one of the main research areas in the field of American options, well-known approximation formulas can be found in the literature, such as the approximations from Zhu (2006), Stamicar et al. (1999), and Evans et al. (2002). Each of those formulations can be considered to gain the initial values for the free boundary. Since the penalisation is based on the initial choice of the free boundary, the choice of the formula for the computation of the initial free boundary has a large effect on the accuracy of the method. A detailed analysis of the choice of the initial free boundary computed by different approximation formulas from the literature is part of our future research. After computing the penalty term with the initial free boundary value, we solve the penalised heat Equation (13). The obtained solution is the solution of the American option problem.

3. Discretization

Let us introduce a temporal discretization ${\tau }_j=T-j\Delta \tau$, $\Delta \tau =T/M$, $j=0,\dots ,M$, and a spatial grid between the points $x_{min}$ and $x_{max}$

$$x_i=x_{min}+i\Delta x,\Delta x=\frac{x_{max}-x_{min}}{N},i=0,\dots ,N.$$

We use the finite difference $\theta$-scheme for discretization and simplify the notation by

$${\alpha }_1=\theta \frac{\Delta \tau }{{(\Delta x)}^2}\mathrm{and}{\alpha }_2=(1-\theta )\frac{\Delta \tau }{{(\Delta x)}^2}\mathrm{with}0\le \theta \le 1.$$

We obtain $w^j={(w_1^j,\dots ,w_{N-1}^j)}^{\top }$ with $w_i^j$ as the approximation for $u(x_i,{\tau }_j)$, $f^j={(f_1^j,\dots ,f_{N-1}^j)}^{\top }$ with $f_i^j\sim f(x_i,{\tau }_j)$ and the diagonal matrices A and B

$$A=diag(-{\alpha }_1,2{\alpha }_1+1,-{\alpha }_1),B=diag({\alpha }_2,-2{\alpha }_2+1,{\alpha }_2),$$

as well as the vector $d^j$ containing the boundary values

$${(d^j)}^{\top }=\left({\alpha }_1{w}_0^{j+1}+{\alpha }_2{w}_0^j,0,\dots ,0,{\alpha }_1{w}_{N+1}^{j+1}+{\alpha }_2{w}_{N+1}^j\right).$$

The discretized penalty term (13) is given by

$$\begin{array}{c}\hfill p^j={\delta }^j·g^j\mathrm{with}{\delta }^j=\left(a{\tau }_j+\frac{{\sigma }^2}{2}\right)\mathrm{and}\mathrm{with}\\ \hfill g_i^j=\left\{\begin{array}{cc}kexp\left(-\alpha x_i-\beta {\tau }_j\right),\hfill & x_i<{\overline{x}}_f^j,\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}\right.\end{array}$$

where ${\overline{x}}_f^j$ is the unique solution to the European put option problem

$$\overline{u}\left({\overline{x}}_f(\tau ),\tau \right)=f\left({\overline{x}}_f(\tau ),\tau \right),$$

(15)

where $\overline{u}$ is a solution to the Cauchy problem: ${\overline{u}}_{\tau }-{\overline{u}}_{xx}=0$, $\overline{u}(x,0)=f(x,0)$, which is given in a closed explicit form. Including all the components, we obtain the $\theta$-scheme discretized formulation for the penalised heat Equation (12)

$$A{w}^{j+1}-B{w}^j-d^j=\Delta \tau ·p^j,$$

where the multiplication of $\Delta \tau$ results from the discretization.

4. Numerical Results

In this section, we consider the example for pricing American put options from Nielsen et al. (2002). All results were computed on an Intel^® Core™ i7-5557U CPU with 3.10 GHz. We chose $x_{\min }=-4$, $x_{\max }=4$, $M=5000$, and used the parameter sets from

Table 1. Numerical results of the corresponding parameter sets.

Example	T	K	r	$\mathit{\sigma }$	N	$\overline{\mathit{a}}$ (×10−4)	${\mathit{S}}_{\mathit{f}}$	MSE
1	3	100	0.08	0.2	1000	7.5	81.87	$9.6\times {10}^{-3}$
					2500	7.4	82.00	$6.2\times {10}^{-3}$
2	1	1	0.1	0.2	1000	10.2	0.862	$5.2\times {10}^{-5}$
					2500	10.0	0.863	$3.3\times {10}^{-5}$
3	0.05	10	0.1	0.25	1000	8.0	9.158	$3.5\times {10}^{-5}$
					2500	8.5	9.142	$3.5\times {10}^{-5}$
4	0.1	100	0.1	0.3	1000	5.5	86.59	$1.2\times {10}^{-3}$
					2500	5.4	86.87	$7.6\times {10}^{-4}$
5	1	100	0.1	0.4	1000	2.6	66.49	$1.4\times {10}^{-2}$
					2500	2.63	66.60	$9.2\times {10}^{-3}$
6	0.05	50	0.1	0.4	1000	3.0	42.61	$3.8\times {10}^{-4}$
					2500	3.1	42.61	$2.5\times {10}^{-4}$

. To facilitate the optimization, we summarised $\Delta \tau p^j$ as

$$\Delta \tau p^j=\left(\overline{a}{\tau }_j+\Delta \tau ·r\right)·\left\{\begin{array}{cc}exp\left(-\alpha x_i-\beta {\tau }_j\right),\hfill & x_i<{\overline{x}}_f^j,\hfill \\ 0,\hfill & \mathrm{otherwise},\hfill \end{array}\right.$$

(16)

where $\overline{a}=\Delta \tau ·k·a$. Through the initial guess, the only unknown parameter is $\tilde{a}$. Since a deterministic expression for $\tilde{a}$ is a goal of our future research, the penalty parameter $\tilde{a}$ was obtained by an optimization. The optimization was done by minimizing the mean square error (MSE) of the solution corresponding to $\tilde{a}$. The MSE is given by

$$MSE=\frac{1}{N}\sum _{i=1}^N{(P^{\mathrm{PSOR}}(S_i,0)-P^{\mathrm{Pen}}(S_i,0))}^2,S_i=Kexp(x_i),$$

where $P^{\mathrm{PSOR}}$ is the solution obtained by the projected SOR algorithm and $P^{\mathrm{Pen}}$ is the maximum of the solution of the penalised system and the payoff function. The maximum was used to gain comparable results to the PSOR algorithm. Since the parameter set 2 is widely used in research, we compared the free boundary value of the parameter set 2 with the free boundary solution of Nielsen et al. (2002), Fazio et al. (2019), and Company et al. (2014). The value for the free boundary solution obtained by Nielson is 0.8622, Fazio obtained 0.86274, and the free boundary value of Company et al. (2014) is 0.8628, where our approach gives 0.86269 with a finer optimization and $\tilde{a}=10.11\times {10}^{-4}$. This comparison illustrates the high accuracy of this method in comparison with related research.

Since the PSOR method does not consider a penalty term, we also computed the MSE of the solution of the penalised PDE

$$u_{\tau }-u_{xx}=\frac{1}{\widehat{\delta }}(f-u),$$

(17)

using a well-known penalty term (Günther and Jüngel 2010) with $\widehat{\delta }=1\times {10}^{-4}$ and $N=1000$. The results are presented in

Table 2. MSE for the parameter sets from the Table 1 computed with the penalty term (17).

1	2	3	4	5	6
$3.3\times {10}^{-1}$	$1.8\times {10}^{-3}$	$1.1\times {10}^{-3}$	$2.2\times {10}^{-2}$	$1.8\times {10}^{-1}$	$5.5\times {10}^{-4}$

, and the comparison of the numerical results illustrate the accuracy of the method. As we included $\Delta \tau$ in $\overline{a}$, the obtained values for $\overline{a}$ are small, and the corresponding values for $\tilde{a}$ are larger than 1. The best results were obtained by the sample sets with little volatility and short time maturity. The observation of the short time maturity is based on the fact that the number of points were different. The dependence on the volatility was caused by the simplification of the term p, since we cancelled out ${\sigma }^2/2$ and included $2/{\sigma }^2$ in $\tilde{a}$. We observed that the differences were in the range between an estimated free boundary value and the final free boundary value. They were caused by the time-dependent movement of the free boundary position. There are several ways to analyse this approach in detail. Since a sensitive point of the presented method is the choice of the initial free boundary value, an interesting approach for future research is a detailed analysis of the effect of the choice of the approximation formulas to the solution. As approximation formulas for the initial guess, one can choose the formulas in (Evans et al. 2002; Stamicar et al. 1999; Zhu 2006). Another idea is the consideration of an iterative scheme. In this idea, the free boundary value of the obtained solution is used as an initial guess for a second iteration of solving the penalised American Put problem.

5. Conclusions

The numerical results give clear evidence that using a non-constant penalty parameter $\delta$ is both feasible and beneficial. Future work will focus on the inclusion of the free boundary movement, a deterministic penalty function, and on the extension to multi-asset American options.