Pricing methods for $\alpha$-quantile and perpetual early exercise options based on Spitzer identities

We present new numerical schemes for pricing perpetual Bermudan and American options as well as $\alpha$-quantile options. This includes a new direct calculation of the optimal exercise barrier for early-exercise options. Our approach is based on the Spitzer identities for general L\'evy processes and on the Wiener-Hopf method. Our direct calculation of the price of $\alpha$-quantile options combines for the first time the Dassios-Port-Wendel identity and the Spitzer identities for the extrema of processes. Our results show that the new pricing methods provide excellent error convergence with respect to computational time when implemented with a range of L\'evy processes.


Introduction
Much of the recent literature in finance concerns exotic options where the payoff depends on the path of the underlying asset. Fluctuation identities such as the ones by Spitzer have applications in the pricing of many of these contracts. Here, we introduce novel methods for pricing perpetual Bermudan and American options, and α-quantile options, using these identities.
Option pricing is a classic problem in financial management and has been subject to many different approaches such as the one by Tse et al. (2001) for discretely monitored hindsight and barrier options. However, this approach was limited to Brownian motion and, as demonstrated by Kou (2002) and Du and Luo (2017), a process with both jump and diffusion processes is a more appropriate model. The method we outline here is based on probabilistic identities and, as it is valid for general Lévy processes, it is applicable to a wide range of jump-diffusion models.
Furthermore, many decision-making processes can be modelled mathematically as American options (Battauz et al., 2015) and thus the new methods we present here also have a more general application in finance.
While European options can only be exercised at a single maturity, American options can be exercised at any time up to expiry. Bermudan options are halfway between the two in that they can be exercised at a finite set of dates, i.e they are a discretely monitored version of American options. These can have an expiry date, beyond which the contract is worthless, or no expiry date, in which case they are called perpetual options. The valuation of American options is a long-standing problem in mathematical finance (Merton, 1973;Brennan and Schwartz, 1977) as it combines a pricing problem with an optimisation problem (i.e the computation of the early exercise barrier) and a pricing problem. A closed-form solution has not been found for a finite expiry date.
In contrast, perpetual American options permit a closed-form solution when the underlying is driven by geometric Brownian motion (Merton, 1973), as the perpetual nature of the option means that the optimal exercise barrier is constant rather than a function of time to expiry. However, this closed-form solution cannot be extended to Bermudan options. Moreover, as explained by Boyarchenko and Levendorskii (2002a), the smooth pasting method used to price perpetual American options can fail under jump processes such as those in the general Lévy class.
Several approximate methods suggested for fixed-expiry American options, such as finite differences (Brennan and Schwartz, 1977), trees (Cox et al., 1979), Monte Carlo (Rogers, 2002) and recursive Hilbert transforms (Feng and Lin, 2013) are inherently discrete and thus lend themselves to Bermudan options with finite expiry. However, they are not particularly accurate or efficient for perpetual Bermudan options as the computational load increases with the number of monitoring dates which is, of course, infinite for these types of contract. Moreover, in times of very low interest rates such as the last decade, the existence of exercise dates even many years in the future can have a significant effect on the option value and thus the number of dates cannot generally be truncated to a sufficiently low value for these methods to be computationally efficient for perpetual options. Boyarchenko and Levendorskii (2002a) published a method to price perpetual American options for many Lévy processes using analytic approximations of the Wiener-Hopf factors. This was a step forward in showing the applicability of the Wiener-Hopf method to price perpetual options with Lévy processes. However, it is an approximate solution as there is no general closed-form method to calculate the Wiener-Hopf factors; in addition, the proposed method is not applicable to all Lévy processes and specifically excludes the variance gamma (VG) process. In Boyarchenko and Levendorskii (2002b) this method was adapted to perpetual Bermudan options; however, the Wiener-Hopf factorisation again requires approximation in some cases and calculations were presented for simple jump-diffusion and normal inverse Gaussian (NIG) processes only. Mordecki (2002) also devised a pricing approximation for perpetual American put options with Lévy processes which is based on the optimal stopping problem for partial sums by Darling et al. (1972); therefore it intrinsically operates in discrete time and thus is useful for Bermudan options. However, this method has restrictions on the jump measure used in the Lévy-Khinchine representation of the characteristic function and therefore cannot be used for general Lévy processes.
Here we fill the gap in the literature for a direct numerical pricing method for perpetual Bermudan and American options which can be used for general Lévy models. In Section 4 we describe a novel a pricing scheme based on the Spitzer identities which includes a new way to directly calculate the optimal exercise barrier. We also provide the first implementation of the method by Green (2009) based on the Spitzer identities and an expression for prices of first-touch and overshoot options which uses residue calculus.
In this article we also deal with quantile options which belong to the class of hindsight options. These have a fixed expiry where the payoff at expiry is determined by the path of the option up to that date. Two such examples are lookback options, as priced by Fusai et al. (2016), and quantile options. Fixed-strike lookback options have a similar payoff to European options except that, instead of being a function of the underlying asset price at expiry, it uses the maximum or minimum over the monitoring period. Quantile options can be considered an extension of lookback options because, rather than using the maximum or minimum of an asset price, the payoff is based on the value which the asset price spends α% of the time above or below. For this reason they are often described as α-quantile options.
They were first suggested by Miura (1992) as a way to design a hindsight option which is less sensitive to very extreme, but short lived, phenomena in the asset price process compared to simple lookback options. Akahori (1995) and Yor (1995) published analytic methods for pricing these contracts with Brownian motion. Since then, most work on pricing these types of options has been based on the remarkable identity by Dassios (1995) which used work by Wendel (1960) and Port (1963) and relates the probability distribution of the α-quantile to the probability distribution of the maximum and minimum of two independent processes. A note by Dassios (2006) showed that this can also be extended to general Lévy processes. Several solutions for Lévy processes have been developed, such as Monte-Carlo methods for jumpdiffusion processes by Ballotta (2002) and an analytic method for the Kou (2002) doubleexponential process by Cai et al. (2010). For discrete monitoring, Fusai and Tagliani (2001) explored the relationship between continuously and discretely monitored α-quantile options however, whilst they made recommendations for selecting an optimal pricing method depending on the value of the monitoring interval ∆t, they did not produce a correction term like the one by Broadie et al. (1997) for barrier options. Atkinson and Fusai (2007) provided closed-form prices for discretely monitored quantile options by using the z-transform to write the problem in terms of a Wiener-Hopf equation. They solved this analytically for prices driven by geometric Brownian motion and also demonstrated the relationship between their results and the Spitzer identities. This approach was extended to general Lévy processes by Green (2009), who developed direct methods, based on his formulation of the Spitzer identities, for calculating the distribution of the supremum and infimum of processes. These methods were implemented for lookback options in Fusai et al. (2016) which showed exponential convergence for general Lévy processes with a CPU time which is independent of the number of monitoring dates. In this paper we address the gap in the literature for a general method for pricing α-quantile options with general Lévy processes and a date-independent computational time.
The article is organised as follows: Section 2 provides the background to transform techniques, with special reference to the use of the Fourier transform in option pricing. Sections 3, 4 and 5 describe the new techniques for pricing α-quantile, perpetual Bermudan and perpetual American options and present results for these new methods for a range of Lévy processes. All the numerical results were obtained using MATLAB R2016b running under OS X Yosemite on a 2015 Retina MacBook Pro with a 2.7GHz Intel Core i5 processor and 8GB of RAM following the step-by-step procedures described in Appendix A and using the process parameters listed in Appendix B.

Transform methods for option pricing
In this paper we make extensive use of the Fourier transform (see e.g. Polyanin and Manzhirov, 1998;Kreyszig, 2011), an integral transform with many applications. Historically, it has been widely employed in spectroscopy and communications, therefore much of the literature refers to the function in the Fourier domain as its spectrum. According to the usual convention in financial literature, we define the forward and inverse Fourier transforms as Let S(t) be the price of an underlying asset and x(t) = log(S(t)/S 0 ) its (rescaled) log-price. To find the price v(x, t) of an option at time t = 0 when the initial price of the underlying is S(0) = S 0 , and thus its log-price is x(0) = 0, we need to discount the expected value of the undamped option payoff φ(x(T ))e −α d x(T ) at maturity t = T with respect to an appropriate riskneutral probability distribution function (PDF) p(x, T ) whose initial condition is p(x, 0) = δ(x).
Here α d is the damping parameter, the selection of which was described in more details by Feng and Linetsky (2008). As shown by Lewis (2001), and comprehensively used in Fusai et al. (2016) and Phelan et al. (2018b), this can be done using the Plancherel relation, Here, p * (ξ + iα d , T ) is the complex conjugate of the Fourier transform of e −α d x p(x, T ). To price options using this relation, we need the Fourier transforms of both the damped payoff and the PDF. A double-barrier option has the damped payoff where e α d x is the damping factor, θ = 1 for a call, θ = −1 for a put, 1 A (x) is the indicator function of the set A, k = log(K/S 0 ) is the log-strike, u = log(U/S 0 ) is the upper log-barrier, l = log(L/S 0 ) is the lower log-barrier, K is the strike price, U is the upper barrier and L is the lower barrier. The Fourier transform of the damped payoff φ(x) is available analytically, where for a call option a = u and b = max(k, l), while for a put option a = l and b = min(k, u). The Fourier transform of the PDF p(x, t) of a stochastic process X(t) is the characteristic function For a Lévy process the characteristic function can be written as Ψ(ξ, t) = e ψ(ξ)t , where ψ(ξ) is known as the characteristic exponent and uniquely defines the process. The fluctuation identities we use in the pricing methods for discretely monitored options are expressed in the Fourier-z, with the Fourier transform applied to the log-price domain and the z-transform applied to the discrete monitoring times. The z-transform is defined as where f (n) = f (t n ), n ∈ N 0 , is a function of a discrete variable. For some of the pricing methods we also require the inverse z-transform. This is not generally available in closed form, so we use the well established method by Abate and Whitt (1992) which has been used for option pricing by Fusai et al. (2006Fusai et al. ( , 2016; Phelan et al. (2018b).
For continuously monitored options the Laplace transform is applied to the time domain. The forward Laplace transform is an integral transform similar to the Fourier transform and it isf Similarly to the discrete time case, the inverse transform is not generally available in closed form and so we use the numerical inverse Laplace transform by Abate and Whitt (1995), used for option pricing by Phelan et al. (2018a). The relationship between the forward Laplace and z transforms is This relationship can be exploited to derive versions of the fluctuation identities with continuous monitoring (Baxter and Donsker, 1957;Green et al., 2010;Fusai et al., 2016;Phelan et al., 2018a).

Decomposition and factorisation
The fluctuation identities require the decomposition and factorisation of functions in the Fourier domain. Decomposition splits a function in the Fourier domain into two parts which are nonzero in the log-price domain only above or below a fixed value respectively, i.e.
We use the Plemelj-Sokhotsky relations for the decomposition, where H[·] denotes the Hilbert transform, which is implemented numerically using the sinc method described by Stenger (1993) and used by Feng and Linetsky (2008) and Fusai et al. (2016) for option pricing. The latter paper describes the implementation of this method to compute fluctuation identities for option pricing; a detailed error analysis of the technique within this application was subsequently carried out by Phelan et al. (2018a,b). The calculations also require the factorisation of a function, i.e obtaining g ⊕ (ξ) and g (ξ) such that g(ξ) = g ⊕ (ξ) g (ξ). This is done by decomposing the logarithm of g(ξ) and then exponentiating the result.

Quantile options
The α-quantile X α of a random process X(t) is the value which it is below α% of the time; α-quantile options have a payoff which depends on this value. They were designed by Miura (1992) as a variation on lookback options which would be less susceptible to the effects of very large, but short lived, swings in the price of the underlying asset. They are therefore more resistant to market manipulation as it is easier to cause a brief swing in an asset price compared to a longer term price movement. We propose a pricing method for α-quantile options with general exponential Lévy processes.

Pricing discretely monitored quantile options
For α-quantile options the form of the payoff is the same as in Eq. (4), but in this case it is calculated as a function of the quantile, i.e. x = X α rather than the process value at expiry x = X(T ). Therefore, with the characteristic function of X α , we can price an α-quantile option in the Fourier domain using the Plancherel relation in Eq.
(3) and the Fourier transform of the damped payoff φ(ξ) from Eq. (5). The Dassios-Port-Wendel identity (Wendel, 1960;Port, 1963;Dassios, 1995) states that the α-quantile of a Brownian motion over the interval [0, T ] has the same distribution as the sum of the infimum of a Brownian motion over the time [0, 1 − αT ] and the supremum of an independent Brownian motion over the time where X α is the α-quantile of the Brownian motion. We can understand the link between the value of X α and the supremum and infimum intuitively. Firstly, if a process is split into a section for t ∈ [0, αT ] and one for t ∈ (αT, T ] then, by the property of independent increments, they represent two independent processes over t ∈ [0, αT ] and t ∈ (0, (1 − α)T ]. Moreover, the supremum of X 1 (t) is the value that the process has spent αT time below and conversely the infimum of X 2 (t) is the value which this process has spent (1 − α)T time above. Although the basic idea behind the relationship in Eq. (13) is quite clear, the mathematical proof is quite involved and we refer the interested reader to the original paper by Dassios (1995) for the details. A note by Dassios (2006) showed that this identity also extends to general Lévy processes. Green et al. (2010) devised Spitzer-based formulations for the probability distributions of the maximum and minimum of a process which were used by Fusai et al. (2016) to price fixedstrike lookback options with exponential error convergence. These are defined in the Fourier-z domain as where Φ ⊕ (ξ, q) and Φ (ξ, q) are the Wiener-Hopf factors of 1−qΨ(ξ, ∆t) as described in Section 2.1. The inverse z-transform can be applied to obtain in the Fourier-domain, where N is the number of discrete monitoring dates and j is the approximation of αN to the nearest integer, although contracts are often written to avoid this approximation by selecting α and N so that αN is an integer. As X M and X m are the maximum and minimum of mutually independent processes, they are mutually independent random variables. It is a basic result in probability theory that the PDF of the sum of two independent random variables is equal to the convolution of their PDFs, i.e.
by convolution theorem this becomes a plain product in the Fourier domain,  Figure 1: PDF of the maximum of a discretely monitored Brownian motion for t > 0 (left) and t ≥ 0 (right) for 50 monitoring dates with σ=0.4, risk-free interest rate r = 0.05 and 2 16 price grid points.
The option price can then be obtained from the Plancherel relation in Eq.
(3) using the Fourier transform of the damped payoff function φ(ξ) from Eq. (5). The calculation of the discretely monitored price for loopback options, as used by Fusai et al. (2016), are based on a distribution of the maximum (or minimum) of the process at t > 0. Similarly to the schemes for barrier options also described in that paper, the first date is taken out of the Spitzer-based scheme and the result for N − 1 dates is multiplied by the characteristic function. This gives a smooth probability distribution, i.e. p X M ∈ C ∞ , as illustrated on the left hand plot of Figure 1, where X M is used to denote that we are using the maximum for t > 0. As described by Ruijter et al. (2015), Fourier based pricing methods using this PDF will therefore not be negatively affected by the Gibbs phenomenon and can achieve exponential error convergence with the number of log-price grid points.
However, for the α-quantile options we price in this paper, we require the distribution for the maximum (minimum) for t ≥ 0. As all Lévy processes have the property that X(0) = 0, the value of the maximum for t ≥ 0 cannot go below 0. Therefore obtaining p X M (x) using the Spitzer-based scheme with the full number of dates alters the PDF so that it now has an abrupt discontinuity and a large spike at t = 0 as shown on the right hand panel of Figure 1. The large spike corresponds to a single probability mass equal to 0 −∞ p X M (x)dx. This can be understood by seeing that the cumulative distribution function (CDF) of X M is If P (X M ≤ 0) = 0 as is generally the case for Lévy processes, then P X M (x) has a jump at Thus when the CDF is differentiated to give the PDF p X M (x), this will contain a probability mass at x = 0 of size 0 −∞ p X M (x)dx. We also note that the introduction of the discontinuity and spike has caused oscillations in the plot of the pdf via the Gibbs phenomenon.
The existence of the discontinuity means that, as described by Boyd (2001) and Gottlieb and Shu (1997) for example, we would no longer obtain exponential error convergence with grid size using these distributions to price options. We therefore use spectral filtering, as successfully implemented for Fourier based option-pricing methods by Ruijter et al. (2015), Phelan et al. (2018b) and Cui et al. (2017). Following it's use in previous literature we use an exponential filter of order 12 (Gottlieb and Shu, 1997).

Results for discretely monitored α-quantile options
We implemented the pricing method described above using the numerical procedure described step-by-step in Appendix A. Figure 2 shows results with α = 0.75 for monitoring dates N of 52, 252 and 1008 for the Gaussian, VG and Merton jump-diffusion processes. The error convergence for this method is extremely fast; we were able to achieve a CPU time of 10 −2 seconds or less for an error of 10 −8 using the Gaussian and Merton processes and 10 −6 using the VG process. There are too few points before reaching the error floor of 10 −11 , caused by the inverse z-transform, to assess whether the error convergence is exponential or high order polynomial. This is not a concern for using the method in practise but a more accurate z-transform would allow us to understand the convergence better.
In Table 1 we show the variation of the price with α for 252 monitoring dates, as expected the option price increases with α, and compare the results with those from Monte-Carlo pricing methods. We use two Monte-Carlo methods, the first simulates N monitoring date paths and selects the αN th smallest value. The other method uses the same approach as Ballotta and Kyprianou (2001) by combining the Monte-Carlo simulation with the Dassios-Port-Wendel identity. Two independent paths of length αN and (1 − α)N dates are simulated and the sum of their respective infimum and supremum are calculated to provide an estimate of the α-quantile. Further details about this Monte-Carlo scheme are included in Appendix C.

Perpetual Bermudan options
Bermudan options have the same payoff as European options but they can be exercised at a discrete set of intermediate dates rather than simply at a final expiry date. They can also be thought of as a discrete version of American options and, indeed, the prices for Bermudan options are often used as a approximation for the value of American options, see e.g. Feng and Lin (2013). Perpetual Bermudan and American options have no expiry date and are therefore "live" until they are exercised. Pricing perpetual options is an easier problem than the valuation of those with finite expiry as the infinite time horizon means the optimal exercise boundary is constant rather than a function of time. Indeed, closed-form formulas exist for perpetual American options with simple processes, whereas there are none for finite expiry options. We look at two methods for pricing perpetual Bermudan options. Firstly we implement the method by Green (2009) which uses residue calculus. We also implement a new method which uses Spitzer identities and show a novel way to calculate the optimal exercise barrier. New numerical truncation bounds for the log-price domain are specified; these are required due to the infinite time horizon.
Both methods require the probability distribution of the current value of a process, subject to the path not having crossed a lower barrier l, i.e.
Using techniques from complex analysis, Green (2009) and Green et al. (2010) showed that this can be expressed using the Spitzer identities as x ≥ l, For pricing single-barrier options, Fusai et al. (2016) used the version of the Spitzer identity for x ≥ l. In contrast, both the methods for pricing Bermudan options described here require the version of the identity for x < l. That is we require the distribution of X(t n ) subject to t n being the first time it has crossed the discretely monitored barrier l. This is illustrated in Figure 3 for a process with 10 monitoring dates and a barrier of −0.2. The value of X(t 10 ) for paths 1 and 2 would contribute to the distribution as the paths stay above l for monitoring dates 0-9 but go below l at the 10 th monitoring date. Note that path 1 is acceptable as, even though its path does go below l between dates 3 and 4, it is above the barrier at the actual monitoring dates. Path 3 does not count towards the distribution as it is above l at the 10 th monitoring date and path 4 does not count as it goes below l at an earlier monitoring date.

Green's residue method
We here briefly recap the method described by Green (2009) to provide a background to the results from our numerical implementation and also because we re-use some of the same techniques in deriving the new Spitzer based method described in Section 4.2. The scheme is based on a combination of a first-touch option paying K − D, where K is the strike, and an overshoot option which pays the difference between the barrier D and the underlying asset S(t) = e X(t) the first time the barrier is crossed. A first-touch option requires the probability that the first time the underlying asset crosses the barrier l is time n∆t, i.e.
where ∆t is the time step between monitoring dates. The next insight by Green (2009) is that the summation and inverse z-transform cancel each other if we use q = e −r∆t because the summation equates to a forward z-transform. Then the price for a first-touch option is The value of the payoff of an overshoot option at time n is the expected overshoot D − S(t n ) conditional on n being the first time that the underlying asset process falls below l, τ l is We can calculate the option value by using tower property to take the expectation over all discrete monitoring and substitute p l (x, n) : The final line uses the Plancherel relation to move into the Fourier domain and we have q = e −r∆t to cancel the inverse z-transform cancel as before. Using Eq. (5) for φ(ξ) for a put option with a = −∞, b = l, k = l, α d > 0 and solving via residue method As the price of a perpetual Bermudan option v B (0, 0) is the sum of a first-touch option with payoff (K − D) from Eq. (24) and the price of an overshoot option from Eq. (28) then Green (2009) showed that, as the optimal exercise barrier gives a maximum price, solving ∂v B (0, 0)/∂D = 0 for D gives the optimal exercise barrier

New formulation based on Spitzer identities
The method by Green described in Section 4.1 is very elegant mathematically but does not particularly lend itself to an intuitive understanding. Therefore we also devised an alternate method for pricing perpetual Bermudan options, including a new way of calculating the optimal exercise boundary. We first recognise that the expected return from exercising a perpetual Bermudan put (subject to us being below the optimal exercise barrier) at monitoring date n can be expressed as To obtain the value of the option we sum over all monitoring dates and we can also substitute p l (ξ, q) := P l− (ξ, q)Φ (ξ, q) for the probability distribution We can again use the trick by Green (2009) of using q = e −r∆t so that the summation cancels with the inverse z-transform to give This integral can then be expressed in the Fourier domain using the Plancherel relation where φ(ξ) is the damped payoff for a put option from Eq. (4) with l being the optimal exercise boundary in the log-price domain. For this method, the calculation of the optimal exercise boundary is based on the idea that if the underlying asset is exactly at the optimal exercise boundary, i.e. S 0 = S(0) = D opt , then the value of the payoff from exercising the option is equal to the expected value from continuing to hold the option. Furthermore, the boundary used to calculate p l (x, n) via p l (ξ, q) := P l− (ξ, q)Φ (ξ, q) is l = log(D opt /S 0 )) and so l = 0 and therefore we can rewrite Eq. (32) as If we differentiate the expression on the first line of Eq. (36) with respect to S 0 , we obtain which is constant. Therefore, if we were to plot Eq. (36) against S 0 we obtain a straight line and the point where that line crosses the payoff function (K − S 0 ) represents the optimal exercise barrier. This is illustrated in Figure 4. Therefore, by calculating Eq. (36) for two values of S 0 (D 1 and D 2 ) we can obtain the corresponding straight line equation with gradient m and y-axis intercept c. We can then calculate the optimal exercise boundary D opt , corresponding to the point where the lines cross, as We can also speed up the computational time by noting that if we set D 1 = 0 in Eq. (36) we obtain the price which means that we only need to perform the inverse Fourier transform for the other calibration point(D 2 ). To avoid computational errors, rather than calculating the Spitzer identity with l = 0 we select l = l , where l is the log-price domain grid step size spacing. This value does not depend on S 0 and so the calculation of the gradient in Eq. (37) still returns a constant. The option price is then calculated using Eq. (35) with l = D opt /S 0 .

Calculation of truncation limits
For the Fourier based methods used for finite expiry boundary options, the range of the log-price domain grid for the numerical calculations was based on the cumulants of the distribution over a single time step and, for the parameters used by Fusai et al. (2016) and Phelan et al. (2018b), were calculated as approximately ±2. However, for perpetual options we must consider the shape of the probability distribution far in the future. This is especially true when the risk-free rate r is low as the contribution from future dates is discounted away extremely slowly. Therefore, the truncation limits used for finite expiry options are far too narrow for this application. We base our calculation of the new limits on the idea that the we should truncate the log-price domain at the value where the discount factor means any distortion of the distribution of the underlying asset process will have negligible effect on the final price calculation. We select a value that we wish the error to be below, i.e. 10 −λ , and calculate the time it will take for the discount factor to be below this value T bound = λ log(10)/r.
We then approximate the standard deviation of the underlying process at this time with where σ is the underlying process volatility, normalised to 1 year. The bounds of the log-price domain [−b, b] are now given by Figures 5 and 6 show results for the two methods for pricing perpetual Bermudan options. The results labelled "Green" are from the residue method described in Section 4.1 and those labelled "Spitzer" are from the new method described in Section 4.2. Results are presented for error vs. CPU time for risk-free rates of 0.02 and 0.05. Using a Gaussian process for the underlying asset means the results as ∆t → 0 can be compared with closed-form calculations for perpetual American options. The convergence of the price for Bermudan options to the continuous case is shown in Tables 2-3 and we discuss the further verification of these results in Section 4.5.

Results for perpetual Bermudan options with the Gaussian process
Both methods perform well, and as ∆t → 0, the results approach those for the corresponding perpetual American option. The new Spitzer based method outperformed Green's residue method, achieving an error of 10 −7 in approximately one tenth of the CPU time. In contrast, it is interesting to note that the convergence of the barrier is slower for the Spitzer based method than Green's residue method. Moreover the barriers of both methods converge at the same rate or slower than the price so we can see that the barrier error has a limited effect on the price error. This can be understood by considering Figure 4. In this plot (K − S 0 ) represents the value of exercising the option, whereas v(0, 0) represents the value keeping the option. Close to D opt this difference is extremely small and so a minor error in D opt has minimal effect on the price. Absolute price error  Figure 5: Results for the price and barrier error convergence with CPU time with an underlying asset driven by a Gaussian process with a risk-free rate r = 0.02. Notice that the pricing error convergence for the new method described in Section 4.2 labelled "Spitzer" is faster than for the residue method described in Section 4.1 labelled "Green", whereas the barrier error convergence is worse.   Figures 7-10 show results for the price and optimal exercise barrier error vs. the number of grid points and CPU time, with the underlying asset driven by the VG and Merton processes. Once again, the Spitzer-based method performs better, achieving an error of 10 −7 about 10 times quicker than Green's method. The error is calculated as the precision compared to the result with the maximum number of grid points and we discuss the verification of these results below. Absolute barrier error Spitzer, t=0.1 Green,t=0.1 Spitzer,t=0.25 Green,t=0.25 Spitzer, t=1 Green, t=1 Figure 6: Results for the price and barrier error convergence with CPU time with an underlying asset driven by a Gaussian process with a risk-free rate r = 0.05. Notice that the pricing error convergence for the new method described in Section 4.2 labelled "Spitzer" is faster than for the residue method described in Section 4.1 labelled "Green", whereas the barrier error convergence is worse.  Table 3: Barrier calculated using both methods for perpetual Bermudan options with 2 20 price grid points showing the convergence to the barrier for a perpetual American option.

Spitzer
In Tables 2 and 3 we showed that that the results with the Gaussian process converge to the closed-form solution by Merton (1973) as ∆t → 0. However, we do not have a general closedform solution for other Lévy processes, so a Monte-Carlo method was used as an approximation.
We wrote Monte-Carlo pricing procedures with underlying assets driven by Gaussian, VG and Merton jump-diffusion processes. Although the discrete nature of perpetual Bermudan options is appropriate for a Monte-Carlo simulation, the absence of an expiry date means that a Monte-Carlo scheme with a finite number of dates will not be a true representation of the contract. However, in our simulation we truncate the Monte-Carlo simulation far enough in the future that the effect of disregarding these future dates is less than the standard deviation of the Monte-Carlo method itself (clearly this method is more feasible for high discount factors and large time steps, as the effect of future dates is discounted away more rapidly). The calculation of the optimal exercise barrier uses the same philosophy as the new method described in Section 4.2, i.e. the price was calculated with S 0 = D 1 and S 0 = D 2 and we found the intersection between the line through these points and the straight line for the payoff K − S 0 .
We calculated an approximate 95% confidence interval (corresponding to 2 standard deviations) for the Monte-Carlo methods and Table 4 shows that the results for our new methods described in Sections 4.1 and 4.2 are within this range for all cases tested. Furthermore, we can see that the results for Green's method and the new Spitzer based method are the same. More details about this Monte-Carlo scheme are included in Appendix C. Absolute barrier error Spitzer, t=0.1 Green,t=0.1 Spitzer,t=0.25 Green,t=0.25 Spitzer, t=1 Green, t=1 Figure 7: Results for the price and barrier error convergence with CPU time with an underlying asset driven by a VG process with a risk-free rate r = 0.02. Notice that the pricing error convergence for the new method described in Section 4.2 labelled "Spitzer" is faster than for the residue method described in Section 4.1 labelled "Green", whereas the barrier error convergence is worse. Absolute barrier error Spitzer, t=0.1 Green,t=0.1 Spitzer,t=0.25 Green,t=0.25 Spitzer, t=1 Green, t=1 Figure 8: Results for the price and barrier error convergence with CPU time with an underlying asset driven by a VG process with a risk-free rate r = 0.05. Notice that the pricing error convergence for the new method described in Section 4.2 labelled "Spitzer" is faster than for the residue method described in Section 4.1 labelled "Green", whereas the barrier error convergence is worse. Absolute barrier error Spitzer, t=0.1 Green,t=0.1 Spitzer,t=0.25 Green,t=0.25 Spitzer, t=1 Green, t=1 Figure 9: Results for the price and barrier error convergence with CPU time with an underlying asset driven by a Merton process with a risk-free rate r = 0.02. Notice that the pricing error convergence for the new method described in Section 4.2 labelled "Spitzer" is faster than for the residue method described in Section 4.1 labelled "Green", whereas the barrier error convergence is worse. Absolute barrier error Spitzer, t=0.1 Green,t=0.1 Spitzer,t=0.25 Green,t=0.25 Spitzer, t=1 Green, t=1 Figure 10: Results for the price and barrier error convergence with CPU time with an underlying asset driven by a Merton process with a risk-free rate r = 0.05. Notice that the pricing error convergence for the new method described in Section 4.2 labelled "Spitzer" is faster than for the residue method described in Section 4.1 labelled "Green", whereas the barrier error convergence is worse.  Table 4: Comparison between the value of a perpetual Bermudan option with 2 20 price grid points and the value for the same contract using a Monte-Carlo approximation. Notice that the prices calculated using the new Spitzter based method and Green's method are the same and within two standard deviations of the Monte-Carlo price.

Perpetual American options
As described by Green et al. (2010), and implemented in Phelan et al. (2018a), we can exploit the relationship between Laplace and z-transforms described in Eq. (9) to extend the pricing methods from discrete to continuous monitoring (i.e. to perpetual American options). Unlike the previous examples for option pricing with continuous monitoring, where the application of option pricing was used as a motivating example for techniques which have relevance for other fields, the continuous (i.e. American) case is commonly used in financial contracts. By defining q = e −s∆t , we can write the continuously monitored equivalent to Φ(ξ, q) = 1 − qΨ(ξ, ∆t) as Φ c (ξ, s) = s − ψ(ξ), where ψ(ξ) is the characteristic exponent of the characteristic function Ψ(ξ, ∆t). Here, as q = e −r∆t then s = r. Both methods described in Sections 4.1 and 4.2 are converted to continuous monitoring and we compare the results below. For assessing the accuracy of the methods with a Gaussian process we have the advantage over Bermudan options that closed-form formulas exist for the calculation of both the barrier and option price and so the results in Figure 11 use the closed-form calculation from Merton (1973) to calculate the error. For the other process we do not have a closed-form result, and the continuous nature of American options means that they cannot be accurately represented using Monte-Carlo methods which are inherently discrete. Therefore the absolute errors displayed in Figures 12 and 13 are calculated against the result for the same method and process with the maximum number of FFT points.

Results for perpetual American options
In contrast to the results for perpetual Bermudan options, the performance of the two methods are very different. Our new Spitzer-based method is far superior, giving errors of approximately 10 −6 in 10 −2 seconds or less. For the Gaussian and Merton processes, Green's method fails to reach an error level of 10 −6 , and for the VG process it reaches this level about 100 times slower than our new method. Absolute barrier error Spitzer, r=0.02 Green,r=0.02 Spitzer,r=0.05 Green,r=0.05 Spitzer, r=0.1 Green, r=0.1 Figure 11: Results for the price and barrier error convergence with CPU time with an underlying asset driven by a Gaussian process with a range of risk-free rates. The error is calculated using the closed-form expression by Merton as a reference. Notice that the price and barrier error convergence for the new method described in Section 4.2 labelled "Spitzer" is faster than for the residue method described in Section 4.1 labelled "Green".   Figure 12: Results for the price and barrier error convergence with CPU time with an underlying asset driven by a VG process with a range of risk-free rates. The error is calculated using the numerical result with the maximum grid size as a reference. Notice that the price and barrier error convergence for the new method described in Section 4.2 labelled "Spitzer" is faster than for the residue method described in Section 4.1 labelled "Green".  Figure 13: Results for the price and barrier error convergence with CPU time with an underlying asset driven by a Merton process and a range of risk-free rates. The error is calculated using the numerical result with the maximum grid size as a reference. Notice that the price and barrier error convergence for the new method described in Section 4.2 labelled "Spitzer" is faster than for the residue method described in Section 4.1 labelled "Green".

Comparison between American and Bermudan option prices
The performance for the direct calculation of the price of American options is sufficiently good for practical purposes, with the Spitzer method having an error of 10 −6 for a CPU time of 10 −2 s. However, it is of academic interest to study the use of the price for Bermudan options as an approximation to the price for American options. The left hand plot in Figure 14 shows the price error compared to American options plotted against ∆t with an underlying asset driven by a Gaussian process and it is clear that the relationship is linear. By extrapolating this line, we can see that in order to achieve an error of 10 −8 for r = 0.02 a step size of ∆t =0.3E-06 is required. The price convergence with this step size is shown in Figure 14 and we can see that reducing the step size this low destroys the monotonicity of the convergence and the excellent error performance that we were seeing for more realistic step sizes.

Conclusion
We implemented three new pricing methods based on the Spitzer identities for pricing exotic options. We saw very fast error convergence with FFT grid size and, by extension, CPU time. Due to the low errors and the presence of the 10 −11 error floor in the inverse z-transform, it was not possible to identify whether this error was exponentially or high order polynomially convergent. The extremely low errors and computational time make this a highly attractive method as it stands. However, it would be of academic interest to implement this method using an inverse z-transform with a lower error floor in order to determine the exact order of convergence which we achieve. We presented a novel method for pricing perpetual Bermudan options which includes a new technique for directly computing the optimal exercise barrier and was based on the Spitzer identities and the sinc-based fast Hilbert transform. We also provide the first numerical implementation of the method designed by Green (2009) which uses residue calculus to remove the requirement for the final inverse Fourier-z transform. Both methods performed well, but our new method showed significantly lower errors and CPU times, with a computational speed ten times faster for an error of 10 −7 .
These methods were extended to perpetual American options (i.e. with continuous monitoring) and very different results were observed for the two methods. For Green's method, the factorisation error was very high and dominated the error convergence. However, for the new Spitzer based method, the factorisation error was much lower and therefore the effect only became dominant for errors below approximately 10 −7 .
For errors greater than 10 −7 , the new Spitzer-based method for perpetual American options has errors at least as low as the method for perpetual Bermudan options so we conclude that, for practical purposes, there is no advantage in using discrete monitoring as an approximation for continuous monitoring. Comparing the two as an academic exercise showed that there is as a linear relationship between the size of the time step ∆t with discrete monitoring and the error compared to continuous monitoring. However, reducing ∆t to a sufficiently low value that the discrete method would be predicted to have a lower error than the continuous method significantly degraded the error convergence of the discretely monitored method such that there was no advantage gained.

A.4 Pricing procedure for perpetual American options
For both methods we adapt the pricing procedure for discretely monitored options to continuous monitoring by replacing Steps 1-2 in the procedures described in Sections A.2 and A.3 with 1. Compute the characteristic exponent ψ(ξ + iα d ) of the underlying transition density.
We then continue with the calculations in Section A.2 or A.3 as before but with Φ c (ξ, r) in place of Φ(ξ, e −r∆t ). Table 5 contains the process parameters used for the numerical tests.

C Monte-Carlo pricing procedures
We describe in more detail the Monte-Carlo pricing procedures which were used to verify the prices calculated by the numerical pricing procedures in this paper.  Table 5: Parameters for the numerical tests; Ψ(ξ, t) is the characteristic function of the process that models the log return of the underlying asset.

C.1 α-quantile options
For underlying assets modelled by Gaussian, VG and Merton jump-diffusion processes, these have paths which can be simply constructed using the built-in Matlab commands for normal, gamma and Poisson random variables. In all cases the risk neutral drift was calculated using µ RN = ψ(−i) as described by Feng and Linetsky (2008), where ψ(ξ) is the characteristic exponent of the Lévy process.
VG process. A single step over time ∆t was simulated using the following procedure: 1. Calculate the gamma-distributed random variable ∆G which has a probability distribution Γ( 1 ν∆t , ν).

Calculate the process step
where ζ is a standard normal random variable.
MJD process. A single step over time ∆t was simulated using the following procedure: 1. Calculate the Poisson distributed random variable with parameter λ∆t.
2. Calculate the process step: where ζ 1 and ζ 2 are independent standard normal random variables.
We used two methods to calculate the value for the α-quantile. For the first we simulated a path of N points and then found the j th smallest value. For the second we combined the Dassios-Port-Wendel identity with the Monte-Carlo method as in Ballotta and Kyprianou (2001). Two independent paths of length αN and (1−α)N dates are simulated and the sum of their respective infimum and supremum are calculated to provide an estimate of the α-quantile.

C.2 Perpetual Bermudan options
We wrote Monte-Carlo pricing procedures with underlying assets driven by Gaussian, VG and Merton jump-diffusion processes using the same techniques to simulate the paths as described in Section C.1.
Although the discrete nature of perpetual Bermudan options is appropriate for a Monte-Carlo simulation, the absence of an expiry date means that a Monte-Carlo scheme with a finite number of dates will not represent the contract exactly. However, we truncate the Monte-Carlo simulation so far in the future that the effect of disregarding these dates is less than the standard deviation of the Monte-Carlo method itself. (Clearly this is more feasible for high discount factors and large time steps as the effect of future dates is discounted away more rapidly.) Finding the optimal exercise barrier uses the same philosophy as the new method described in Section 4.2, i.e. the price was calculated with S 0 = D 1 and S 0 = D 2 and we found the intersection between the line through these points and the straight line for the payoff K −S 0 .