Pricing Perpetual American put options with asset-dependent discounting

The main objective of this paper is to present an algorithm of pricing perpetual American put options with asset-dependent discounting. The value function of such an instrument can be described as \begin{equation*} V^{\omega}_{\text{A}^{\text{Put}}}(s) = \sup_{\tau\in\mathcal{T}} \mathbb{E}_{s}[e^{-\int_0^\tau \omega(S_w) dw} (K-S_\tau)^{+}], \end{equation*} where $\mathcal{T}$ is a family of stopping times, $\omega$ is a discount function and $\mathbb{E}$ is an expectation taken with respect to a martingale measure. Moreover, we assume that the asset price process $S_t$ is a geometric L\'evy process with negative exponential jumps, i.e. $S_t = s e^{\zeta t + \sigma B_t - \sum_{i=1}^{N_t} Y_i}$. The asset-dependent discounting is reflected in the $\omega$ function, so this approach is a generalisation of the classic case when $\omega$ is constant. It turns out that under certain conditions on the $\omega$ function, the value function $V^{\omega}_{\text{A}^{\text{Put}}}(s)$ is convex and can be represented in a closed form; see Al-Hadad and Palmowski (2021). We provide an option pricing algorithm in this scenario and we present exact calculations for the particular choices of $\omega$ such that $V^{\omega}_{\text{A}^{\text{Put}}}(s)$ takes a simplified form.


Introduction
In this paper we consider a perpetual American put option with asset-dependent discounting. We consider a standard stochastic background for this problem, i.e. we define a complete filtered risk-neutral probability space (Ω, F, {F t } t≥0 , P), on which we define the asset price process S t . Then F t is a natural filtration of S t satisfying the usual conditions and P is a risk-neutral measure under which the discounted (with respect to the risk-free interest rate r > 0) asset price process e −rt S t is a local martingale. A family of F t -stopping times is denoted by T while E s denotes the expectation with respect to P when S 0 = s = e x . The value function of the perpetual American put option with asset-dependent discounting can be represented by The asset-dependent discounting is reflected in the ω function what is our key concept considered in this article. We underline here that the discount function ω for various economical reasons can be different from the risk-free interest rate r > 0; see [1] for further explanations. The way we choose discounting is to model strong dependence of discount factor with the asset price. The goal is to understand various economical phenomena that might appear in this extreme case. Our approach differs from typical studies considered in the literature, where the interest rate is independent from the asset price or there is a weak dependence between these two factors. Therefore, this research is noteworthy not only in the context of option pricing, but also in other areas where optimisation problems appear. Moreover, we assume that the asset price process S t is a geometric Lévy process with negative exponential jumps, i.e.
(2) S t := se Xt with (3) where ζ and σ > 0 are constant, N t is the Poisson process with intensity λ ≥ 0 independent of Brownian motion B t and {Y i } i∈N are i.i.d. random variables independent of B t and N t having exponential distribution with mean 1/ϕ > 0. Under the martingale measure P the drift parameter is of the form ζ := r − σ 2 2 + λ ϕ + 1 .
Note that when λ = 0 then we end up with the classical Black-Scholes model. However, empirical studies show that stock prices have heavier left tail than normal distribution. Therefore, nowadays many books and articles concern, as we do in this work, pricing of derivative securities in market models based on Lévy processes; see [3] for more details.
The main objective of this paper is to present an algorithm of pricing perpetual American put options with asset-dependent discounting with the value function defined in (1) and the asset price process S t given in (2). Furthermore, we take into account some specific scenarios (e.g. when σ = 0 or λ = 0) and for these cases we are able to derive analytical forms of the value function, while for more complex examples we show how to handle them numerically. Detailed theoretical results of the analysed problem was already developed in [1], where the authors presented the approach of deriving a closed form of value function (1) for even a more general setting than it is considered here. Therefore, in this paper we focus more on numerical side of this problem and analyse in detail few particular cases where more explicit results can be derived.
Still, before we present the option pricing method in our set-up, we recall the most important theoretical issues on which our article is based on. A key step in deriving a closed form of (1) is identifying the form of the optimal stopping rule τ * for which the supremum in (1) is attained. It turns out that under certain conditions on the discount function ω, which are presented in the next section, the value function is convex. By combining this fact with the classical optimal stopping theory presented e.g. in [10], it allows us to conclude that the optimal stopping region is an interval [l * , u * ] and hence for some optimal thresholds l * ≤ u * . Observe that for the nonnegative discount function ω we have l * = 0 (since waiting is not beneficial). Therefore, in this case a single continuation region appears. In general, for the negative ω we can observe a double continuation region; for more details see [4]. The optimal boundary levels l * and u * can be found by application of standard methods of maximising the over l and u > l. To find v ω A Put (s, l, u) we use exit identities for spectrally negative Lévy processes containing so-called omega scale functions introduced in [7]. As shown in [1,Theorem 9], another way of finding the optimal thresholds l * < u * is to apply the classical smooth and continuous fit conditions. Typically, a price of the option is a solution to a certain Hamiltonian-Jacobi-Bellman (HJB) system and the optimal thresholds are identified using the smooth fit conditions. We want to underline that our approach is different, although still finding the omega scale functions is done via solving certain ordinary differential equations. The paper is organised as follows. In Section 2 we introduce basic theory and notation. Section 2.4 provides main theoretical results of this paper. In Section 3 we present some specific examples where the option price can be expressed in the explicit way. Section 4 focuses on the purely numerical analysis. We also show there that these two approaches are consistent. The last section includes our conclusions.

2.1.
Assumptions. It the beginning, we note that the analysed American put option will not be realised when its payoff is equal to 0. Hence, we can transform the form of the value function given in (1) into the following one We work under the same assumptions as those formulated in [1, Section 2.2]. However, this time we consider slightly more specific assumptions on the ω function, namely that Assumption 1. A discount function ω is concave, nondecreasing and bounded from below.
From [1, Remark 3] we can then conclude that under Assumption 1 the value function V ω A Put (s) is convex.

2.2.
Optimal stopping time. From [1, Section 2.4] it follows that then the optimal exercise time is the first entrance of the process S t into some interval, that is, it has the following form Hence, we can represent value function (4) as Moreover, we denote the optimal stopping time by where l * and u * realise the supremum above. As shown in [1,Theorem 9], another way of identifying the critical points l * and u * can be done via application of the smooth fit property. In that case, l * and u * satisfy 2.3. Scale functions. By applying the fluctuation theory of Lévy processes, we can find a closed form of (5) and hence of (4) in terms of the so-called omega scale functions.
To introduce them formally, firstly let us define the Laplace exponent of X t via which is well-defined for θ ≥ 0 since our X t is a spectrally negative Lévy process. In the case of X t given in (3) the Laplace exponent takes the form By Φ(q) we denote the right inverse of ψ(θ), i.e.
where q ≥ 0. The first scale function W (q) (x) is defined as a continuous and increasing function such that W (q) (x) = 0 for all x < 0, while for x ≥ 0 it is defined via the following Laplace transform We define also the related scale function Z (q) (x) by [2] we know that for X t given in (3) we have where {γ 1 , γ 2 , Φ(q)} is the set the real solutions to ψ(θ) = q. In turn, Z (q) (x) is as follows .
where ξ is an arbitrary measurable function. They are defined as the unique solutions to the following equations where W (x) = W (0) (x) is a classical zero scale function. To simplify notation, we introduce also the following S t counterparts of the scale functions (11) and (12) where ξ • exp(x) := ξ(e x ). For α for which the Laplace exponent is well-defined we can define a new probability measure By [9] and [5, Cor. 3.10], under P (α) , the process X t is again spectrally negative Lévy process with the new Laplace exponent For the new probability measure P (α) we can define the ξ-scale functions which are denoted by the adding subscript α to the regular counterparts, i.e. W Theorem 1. Let Assumption 1 holds and assume that ω is nonnegative. Then the optimal stopping region is of the form (0, u * ] and we have (1) For σ = 0 and λ > 0 (2) For λ = 0 and σ > 0 (3) For σ > 0 and λ > 0 where .
The optimal boundary u * in (18) and (19) can be determined by the smooth fit condition (17) can be determined by the continuous fit condition Remark 2. Let us note that using (13) and (14)   To find the option price V ω A Put (s) we have to identify where ω α u (s) is given in (16) and the case of α = 0 corresponds to W (ωu) (s) and Z (ωu) (s). Observe that we need to find the above ξ-scale functions for ξ = ω α u • exp under measure P (α) , i.e. we have to identify W . This is equivalent to taking our asset price process S t of the form of (2) but with the new parameters given in (15). The second key result for our numerical analysis follows straightforward from [1, Theorem 16] and allows to identify the above omega scale functions using ordinary differential equations.

Option pricing -analytical approach
In this section we present some examples of discount functions for which we are able to determine the analytical form of the value function V ω A Put (s).

3.1.
Constant discount function. The case when ω function is constant, i.e. ω(s) = q is the standard example which appears in the literature quite extensively. However, this case is quite special, as it turns out that the second term of the sum in (19) simplifies and we do not need to deal with the measure P α (and thus to calculate the limit for α → ∞) to find V ω A Put (s). This fact is stated in the below theorem.
Theorem 5. Assume that ω(s) = q. Then Proof. Note that α s u corresponds to the continuous transition of the process S t to the interval (0, u] or, in other words, continuous exit from half-line (u, ∞). We define for details see [1, Proof of Theorem 7].
Ultimately, value function (19) for the constant discount function ω(s) = q can be written as Remark 6. For the case of λ = 0, using (22), one can show that value function (18) simplifies to the well known formula for the value function in the Black-Scholes model, i.e.
where we substituted q = r and ζ = r − σ 2 2 . Therefore, we are not forced to apply the smooth fit condition in order to find the optimal value of u. We can do this analytically by finding the maximum of V ω A Put (s) with respect to u.  Based on Figure 1 we can simply note that a higher value of the discount function ω results in a smaller value of V ω A Put (s) which is in line with the financial intuition. In turn, Figure 2 shows a comparison of (17), (18) and (19) for the same value of q = 0.5. The resulting relation between these functions is again consistent with the economical expectations.

3.2.
Linear discount function. In this subsection we consider a linear discount function of the form ω(s) = Cs for some positive constant C.
3.2.1. σ = 0. Let us consider the case of σ = 0. Then the asset price process S t jumps into the interval (0, u], which means from Theorem 1 that where ω u s u = ω(s) = Cs. Equivalently, (24) can be rewritten as To find a closed form of value function (25) we need to identify the scale functions W (ηu) (x − log u) and Z (ηu) (x − log u). From Theorem 4 it follows that both W (η) (x) and Z (η) (x) solve the following ordinary differential equation where K 1 and K 2 are the constants that can be found based on the initial conditions, is the Kummer confluent hypergeometric function. We denote by K W 1 , K W 2 and K Z 1 , K Z 2 the constants corresponding to W (η) (x) and Z (η) (x), respectively. Using initial conditions (27) and (28) we can simply calculate these constants for both W (η) (x) and Z (η) (x). By shifting these functions by log u we simply produce W (ηu) (x − log u) and Z (ηu) (x − log u). The asymptotic behaviour of 1 F 1 (a, b; t) for t → ∞ is as follows Based on (31) we calculate the constant c Z (η) /W (η) (or equivalenly c Z (ω) /W (ω) ) occuring in (25). It has the following form Combining all the obtained results and substituting them into (25), we can present value function (25) graphically for some sample parameter values. This is done in Section 4.
3.2.2. λ = 0. Let us consider the case of λ = 0. In this case the asset price process S t enters the interval (0, u] in a continuous way only. Therefore, from Theorem 1 .
It suffices to find now W The initial conditions have the following form (35) we obtain the Bessel differential equation of the form The general solution to (38) is equal to and therefore Based on the form of (39) and the fact that D α does not depend on α we can simply note that value function (34) is also independent of α. Therefore, we can take an arbitrary value of α in (35). Thanks to this key observation, its solution (39) can take a simplified form. Indeed, for α = 0 equation (35) is equal to Hence, the general solution to (40) takes the following form If we take the following sample parameters r = 0.05 and σ = 0.2, we obtain n = 2 and therefore Applying initial conditions (36) and (37) we can simply obtain K W 1 , K W 2 and K Z 1 , K Z 2 . Using equality (42) which hold for both W Taking into account all the obtained results, we can obtain value function (34) for sample data. This is done in Section 4 as well.
3.3. Power discount function. This time, we take into account a power function of the form ω(s) = Cs n for n ∈ (0, 1] and C being some positive contant. This case is a generalisation of a linear discount function. 3.3.1. σ = 0. Similarly to the case of a linear discount function, the scale functions W (η) (x) and Z (η) (x) solve with A = C ζ , B = λ−ϕζ ζ and D = C n+ϕ ζ , while the initial conditions are the same as those provided in (27) and (28). Applying a substitution t = A n e nx and F (t) = f (x) we transform (44) into (45) where b = 1 − B n and a = D An . The general solution to (45) has the same form as was provided in (30). Therefore, for both the linear and the power discount function ω, the form of the value function V η A Put (x) is identical.
3.3.2. λ = 0. As in the above case, the idea of finding a closed form of the value function can be borrowed from the linear case. This time, the scale functions W Therefore, we have where v = √ Bα 2 +4Eα n = 2ζ nσ 2 . We can show, as in the prevoius section, that the value function which arises in this scenario does not depend on α. Thus, for α = 0, (46) takes the form Again, having exact formulas for the scale functions, we can easily represent the form of the value function.
All results are presented in Section 4.

Option pricing -numerical approach
In this section we show how to numerically identify the value function V ω A Put (s) for arbitrary discount function ω. We present some figures corresponding to the various discount functions as well.
4.1. Different discount functions. For some discount functions ω we are unable to find the analytical forms of which are the solutions to the ordinary differential equations occurring in Theorem 4. That is, formally we cannot identify explicitly the value function either. In such a situation, we can proceed a numerical analysis of these equations. In general, solving a high-order ordinary differential equation consists in transforming it into first-order vector form and then applying an appropriate algorithm that returns us the numerical solution of the n + 1−dimensional system of first-order ordinary differential equations. For practical purposes, howeversuch as in financial engineering -numeric approximations to the solutions of ordinary differential equations are often sufficient. In this paper we focus on the Higher-Order Taylor Method. This method employs the Taylor polynomial of the solution to the equation. It approximates the 0−th order term by using the previous step's value (which is the initial condition for the first step), and the subsequent terms of the Taylor expansion by using the differential equation.
4.1.1. σ = 0. For the case of σ = 0 we show in the previous section how to derive the value function for the linear discount function ω. Figure 3 presents comparison of the value function given in (24) when the scale functions were calculated analytically (as shown in subsection 3.2.1) and numerically by solving differential ordinary equation (26). In this example, a linear discount function was chosen, i.e. ω(s) = Cs. The difference between these functions is so small that this is negligible. Moreover, Figure 4 illustrates the constant c Z (η) /W (η) obtained in (32) together with the quotient of the functions Z (η) (x) and W (η) (x). As we mentioned at the beginning of this section, we are not always able to get an analytical solution to an ordinary differential equation. This is the case for example when the discount function is of the form ω(s) = C arctan(s) for some positive C. Then we can only obtain the scale functions numerically. Figure  5 shows two value functions, for ω(s) = Cs and ω(s) = C arctan(s), respectively. Since for all positive s we have s > arctan(s), then we expect that the value function corresponding to ω(s) = C arctan(s) takes greater values rather than for ω(s) = Cs. We can also note that the difference between these functions becomes greater for higher values of s, which is in line with economical intuition since the difference between ω(s) = s and ω(s) = arctan(s) also extends as s increases.   section we avoid this problem by selecting discount functions for which the value function is independent of the α parameter. Figure 6 presents comparison of the value function given in (33) for ω(s) = Cs and for the scale functions obtained analytically and numerically. As shown in subsection 3.2.2, in this case we can take an arbitrary value of α and obtain a simplified form of the value function and ordinary differential equation that the scale functions solve. Similarly to the previous example, we again can observe a negligible difference between these two value functions. In turn, Figure 7 shows the constant c Z (η) /W (η) given in (43) with the ratio of Z (η) (x) and W (η) (x). Lastly, in Figure 8 we can observe the value functions for both ω(s) = C √ s and ω(s) = C √ s + Z for some positive Z, i.e. we compare two discount functions that differ in shift. This time, we can see that the value functions obtained in this way also differ only in shift, which confirms financial intuition. 4.1.3. σ > 0 and λ > 0. The most general case is when σ > 0 and λ > 0. Then the considered value function V ω A Put (s) is given by (19). It can be also represented as the function of x variable: For the linear discount function ω(s) = Cs, the scale functions W (η) (x) and Z (η) (x) occurring in (47) are the solutions to the following ordinary differential equation  The initial conditions for W (η) (x) and Z (η) (x) are as follows with A α = γ α2 +γ α3 , B α = C [Υ α2 (γ α2 − γ α3 ) − Υ α1 γ α3 ], D α = −Υ α2 (γ α2 −γ α3 )ψ(α)−γ 2 γ 3 +Υ α1 γ α3 ψ(α), E α = C [Υ α2 (γ α2 − γ α3 ) + Υ α1 γ α2 γ α3 − Υ α1 γ α3 ], F α = −Υ α1 γ α2 γ α3 ψ(α). The initial conditions for W (η α ) α (x) and Z (η α ) α (x) are as follows In this case, we are not able to identify explicit solutions to third order ordinary differential equations (48) and (49). So we are forced to use only the numerical method to find the scale functions and hence the value function. Figure 9 shows several graphs of the value function (19) for different values of α parameter together with the first and second component occurring in (19).

Conclusions
In this paper, we have presented the novel approach to pricing the perpetual American put options with asset-dependent discounting. For the asset price process S t being the geometric spectrally negative Lévy process we have shown that value function (1) can be represented in a closed form based the omega-type scale functions that solve some ordinary differential equations given in Theorem 4. We have used these theoretical results to perform extended numerical analysis for some key financial examples. In particular, for some cases we have managed to produce some explicit formulas for the value function. In the cases where it was impossible to do so, we have used the numerical analysis of the above mentioned ordinary differential equations based on Higher-Order Taylor Method. We have presented many figures of the value functions that arise in various scenarios. One can think of further generalisations. For example when the discount factor is randomised. It can be done in different ways, e.g. by introducing additional Markov economical environment. This type of research is left for future investigations.