Pricing double barrier options on homogeneous diffusions: a Neumann series of Bessel functions representation

This paper develops a novel analytically tractable Neumann series of Bessel functions representation for pricing (and hedging) European-style double barrier knock-out options, which can be applied to the whole class of one-dimensional time-homogeneous diffusions even for the cases where the corresponding transition density is not known. The proposed numerical method is shown to be efficient and simple to implement. To illustrate the flexibility and computational power of the algorithm, we develop an extended jump to default model that is able to capture several empirical regularities commonly observed in the literature.


Introduction
In this paper we develop a novel option pricing methodology based on an analytically tractable Neumann series of Bessel functions (hereafter, NSBF) representation. The new representation is derived by applying the NSBF expansion to the arising Sturm-Liouville problem. To highlight the potential of our method, we derive a new analytical tractable representation of the price (and Greeks) of double barrier European-style knock-out options (henceforth, DBKO options), though applications to other similar problems may be designed using this conceptual framework.
Barrier option contracts are path-dependent exotic options traded at over-the-counter markets on several underlying assets, e.g., stocks, stock indexes, currencies, commodities, and interest rates. They have been actively traded mainly because they are cheaper than the corresponding vanilla options and offer an important tool for risk managers and traders to better express their market views without paying for outcomes that they may find unlikely. Moreover, they are also used as building blocks of many structured products.
Given their popularity in the market, a vast literature dedicated to their valuation has been developed. For instance, alternative pricing (and hedging) schemes for DBKO options under the classic geometric Brownian motion (hereafter, GBM) assumption have been proposed by [25], [16], [43], [36], [41], or [5]. Such modeling framework assumes the volatility is constant throughout the option's life, several attempts have been made to overcome this unrealistic assumption implicit in the GBM diffusion.
It is well-known that the constant elasticity of variance (hereafter, CEV) diffusion model of [9], where the volatility is a function of the underlying asset price, is able to better reproduce two empirical regularities commonly observed in the literature, namely the existence of a negative correlation between stock returns and realized volatility (leverage effect) and the inverse relation between the implied volatility and the strike price of an option contract (implied volatility skew ). To accommodate these observations, the valuation of DBKO options under the CEV model have been performed by [4] through a trinomial scheme, by [11] using a pricing framework based on the numerical inversion of Laplace transforms, by [12] via an eigenfunction expansion approach, and by [30] whose approach rests on the construction of an approximation based on continuous-time Markov chains, amongst others.
More recently, [14] tackle the valuation of DBKO options (using the stopping time approach as well as the static hedging approach) under the so-called jump to default CEV (hereafter, JDCEV) model proposed by [7], which is known to be able to capture the empirical evidence of a positive correlation between default probabilities (or credit default swap spreads) and equity volatility. 1 Moreover, it nests the GBM and CEV models as special cases and, therefore, it also accommodates the aforementioned leverage effect and implied volatility skew stylized facts. The importance of linking equity derivatives markets and credit markets has thus generated a new class of hybrid credit-equity models with the aim of pricing derivatives subject to the risk of default-for other applications of jump to default models, see, for instance, [32], [29], [28], [40], [34], [33], and the references contained therein. Moreover, the new algorithms recently provided by [13] for computing truncated and raw moments of a noncentral χ 2 random variable can be also used on the pricing of barrier options under the JDCEV model considered in [14].
The main purpose of this paper is the development of a new analytically tractable NSBF representation for pricing (and hedging) European-style DBKO options which can be applied to the whole class of one-dimensional time-homogeneous diffusions, independently of knowing the corresponding transition density. Similarly to [12], we solve the boundary value problem for the parabolic partial differential equation using a classical separation of the variables method. This technique reduces the problem to the determination of the eigenvalues and eigenfunctions of the associated Sturm-Liouville problem. The approach based on the NSBF representation allows one to compute large sets of eigendata with non-deteriorating accuracy. Hence, we are able to calculate the prices for the general time homogeneous diffusion models, not relying on the knowledge of the exact solutions as for example is done in [12] for the CEV model and [7] and [14] for the JDCEV model. Therefore, the novel NSBF representation permits the construction of a fast and accurate algorithm for pricing DBKO options and opens its application to other similar problems. 2 We note that [7] are able to obtain closed-form solutions for European-style plain-vanilla options, survival probabilities, credit default swap spreads, and corporate bonds in the JDCEV model by exploring the powerful link between CEV and Bessel processes. By adopting the hybrid creditequity JDCEV architecture modeling framework, [14] are restricted to the volatility and default intensity specifications that are implicit in the JDCEV model. By contrast, since we do not need to be restricted to such specific modeling assumptions we are able to quickly and accurately price DBKO options for a larger class of models. We illustrate our numerical method on extended jump to default constant elasticity of volatility (hereafter, EJDCEV) model, which nests the JDCEV model as a special case.
In summary, our method can be considered as an alternative powerful computational tool for pricing DBKO options. Since we are able to quickly construct the whole value function and not just the price, we can easily observe the behavior of the option price through time and different initial values. Moreover, the NSBF representation also presents an easy way to calculate the derivatives of the value function and consequently the 'Greeks' of the option, and thus it can be used for the design of hedging strategies. Given its accuracy, efficiency, and easy implementation, the novel valuation method can be used also for analyzing the empirical performance of models for barrier option valuation under alternative underlying asset pricing dynamics, e.g. [18].
The remainder of the paper is structured as follows. Section 2 sets the general financial framework and defines the boundary value problem. Section 3 provides the main result of the paper (in Proposition 1): the representation of the solution to the boundary value problem and the price for a DBKO option. Section 4 illustrates the calculation of the so-called 'Greeks'. Next sections are dedicated to the algorithm implementation and numerical examples. First, we present, in section 5, some recurrence formulas which are more robust and efficient for computation of the coefficients that appear in the direct representations of the value function and its derivatives presented in previous sections. Then, section 6 offers the conceptual algorithm for the computation of the price of DBKO options. Section 7 presents numerical experiments for the EJDCEV framework. For illustrative purposes, the analysis is separated into two different horizons: medium (six months) and short (one day). The final section presents the concluding remarks and the possible directions for further research.

Modelling framework
This section presents the general financial model for our pricing method and describes the boundary value problem associated to an option contract containing two barrier (knock-out) provisions. We recall that the holder of such a DBKO option has the right to receive, at the expiration date T , a payoff f (y) = (y − K) + in the case of a call option, or f (y) = (K − y) + for the case of a put option, if the underlying asset price (Y ) remains in the range (L, U ). The real constants U > L > 0 are designated by the upper and lower bounds (i.e., the knock-out trigger barrier levels), whereas K ∈ R : L < K < U is the strike price. 3

The general financial model setup
Hereafter, and during the trading interval [0, T ], for some fixed time T (> 0), uncertainty is generated by a probability space (Ω, G, Q), where the equivalent martingale measure Q associated to the numéraire "money market account" is taken as given. The price dynamics of the underlying asset is assumed to be governed, under the risk-neutral measure Q, by the time-homogeneous (or time invariant) diffusion with initial value Y 0 = y 0 , and where the functions µ (y) and σ (y) are, respectively, the (state dependent) instantaneous drift and instantaneous volatility (whose regularity properties will be formally defined later in Assumption 1), whereas B t ∈ R is a standard Wiener process defined under measure Q.
To explain the development of our pricing methodology we consider a European-style DBKO option contract whose payoff at the expiry date T is a function of the single state variable Y . Given the contractual clauses of such derivative security, the process is considered on the interval [L, U ], where L and U are, respectively, the lower and upper bounds of the DBKO option contract. The end-points L and U are set to be knock-out boundaries, because if at any time between the initial date of the contract and its expiration either the upper barrier or the lower barrier is hit, then the option contract is canceled (i.e. it is knocked out). We are considering the case with no rebate for the simplicity. However, it is possible to incorporate a rebate value, as shown in the note in the appendix.
As in [12], if any of these end-points is a regular boundary, we adjoin a killing boundary condition at that end-point, sending the process to a cemetery state {∆} at the first hitting time of the end-point. Consequently, the hitting time for our problem (with two knock-out provisions) is defined by We also consider the possibility that the process may be killed by a sudden jump to {∆}, i.e. the spot price is allowed to jump to zero from the interior of the interval. This implies that the default event forces the option knock-out. There is no recovery value (in the case of a DBKO put) upon default. This is accomplished by imposing a killing time defined as where h (y) ≥ 0 is the default intensity or the hazard rate, whereas E is an exponential random variable with unit mean, independent of Y and time. Therefore, the combined stopping time of the process is denoted by Remark 1. The abuse of notation is in place; to simplify notation, the stopped process set on the domain [L, U ] ∪ {∆} is denoted by the same letter Y as the original process. This avoids bloating the paper with many technical details that are standard in the literature-see, for instance, [ [29]. We notice that the defaultable asset price process is adapted not to the filtration F = {F t , t ≥ 0} generated by the predefault process, but rather to the enlarged filtration As usual, to ensure that the constructed (killed) process remains a martingale it is necessary to set the drift of equation (1) as wherer (y) andq (y) are the (time-homogeneous) continuously compounded interest rate and dividend yield. In summary, the main purpose of this paper is to develop an efficient and flexible pricing methodology for computing risk-neutral expectations of the form This will be accomplished by applying the NSBF expansion to the associated Sturm-Liouville problem.

The boundary value problem
Let us introduce the differential operator Then, the value function (4) is the solution of the following boundary value problem The pricing of the DBKO option will be performed by solving problem (6). For convenience, we rewrite the operator A in the following form At this point, we can set the needed conditions for the coefficients of the process Y t through the problem (6). In this illustration paper we are looking only at the regular case, so we will need the following assumptions: Assumption 1. The functions p, p ′ , q, w and w ′ are real valued and continuous on [L, U ]. Additionally, it is assumed that p ′ and w ′ are absolutely continuous and that p > 0, σ > 0 and w > 0.

Assumption 2.
The payoff function f is square integrable.
Application of Fourier' separation of variables method to the partial differential equation in (6)-see, for example, [31] and [37] for a general exposition of the method and [12] for financial applications-, leads to the eigenfunction expansion of the value function (4) as where the pairs (λ n , ϕ n ) are solutions of the Sturm-Liouville problem The functions ϕ n form a complete orthogonal basis for the space . The coefficients f n are the Fourier coefficients of the function f relative to the basis {ϕ n } n∈N with scalar product The Hilbert space L 2 w ([L, U ]) with the above defined scalar product is denoted by H w . The function f can be explicitly decomposed as where f n are the Fourier coefficients of the payoff function f defined as We recall that Assumption 2 guarantees that the series converges to the function f in L 2 w norm. We note also that Assumption 1 ensures that problem (9) is a regular Sturm-Liouville problem. The eigenvalues are real, positive and can be listed as λ 1 < λ 2 < .. < λ n < ..., with lim n→∞ λ n = +∞. 4 Remark 2. We further notice that for the put and call barrier options under consideration, problem (6) possesses non-consistent (discontinuous) boundary conditions. For the barrier call option the value function is discontinuous at the point A similar observation can be made for the point (L, T ) in the barrier put case. Nevertheless, in both cases the value function is still in the space H w and can be represented by its Fourier series (11).

An analytical representation through NSBF of the value function
This section presents the pricing formula for the DBKO option using the NSBF representation for the Sturm-Liouville problem (9) recently proposed by [22] for the one-dimensional Schrödinger equation-i.e. the case with w (y) = 1-and then extended to a more general function w (y) in [24]. In a nutshell, this powerful technique consists in the representation of the solutions of the Sturm-Liouville problem (9) and their derivatives in terms of NSBF with explicit formulas for the coefficients.
To extend this approach to our option pricing problem, let us first introduce the space H 1,0 w . This is the subspace of the functions u ∈  (4) is a solution to problem (6) and can be represented as The series converges in the norm of the space H 1,0 w . 5 Moreover the series converges uniformly with Before providing a formal proof of Proposition 1 and for the sake of completeness, let us first highlight some important details aiming to offer a better exposition. Consider the following identities: • The eigenfunctions, solutions to the Sturm-Liouville problem (9) are 6 ϕ n (y) = sin (ω n l (y)) ρ (y) • The spherical Bessel functions of the first kind, j ν (y), are given by where J µ (y) are the Bessel functions of the first kind shown in [17, 8.461.1].
• The function l (y) is defined by 7 • The function ρ (y) is defined by • The roots of the eigenvalues λ n , denoted as ω n , are solutions of the equation • The functions α n (y), n ≥ 0 are defined as The efficient recursive method for computing α n will be presented in Section 5.
• l k,n is the coefficient of x k in the Legendre polynomial of order n-see, for instance, [1,Chapter 8].
• Φ k (y) are the formal powers that will be defined in Definition 1.
The formal powers Φ k (y) are constructed on the basis of one non-vanishing solution g of the with an initial condition set as .
6 Note that these functions are not normalized. 7 A detailed analysis of the role of this transformation in this decomposition and in the transmutation operators theory can be found in [21]. 8  Definition 1. Let p, q, w satisfy Assumption 1 and let g be a non-vanishing solution of equation (16) that satisfies condition (17). Then, the associated formal powers are defined, for k = 0, 1, 2, ..., as where two families of the auxiliary functions are defined as Remark 3. We note that these formal powers arise in the Spectral Parameter Power Series (SPPS) representation for the solution of the Sturm-Liouville problem (9). The SPPS method was introduced in [20]-see also [23] and [19].
Next we provide the formal proof of Proposition 1.
Proof (Proposition 1). The application of the Fourier separation of the variables method to problem (6) leads to representation (8). It is a weak solution of problem (6) (12) and guarantees the approximation of the eigenfunctions uniformly in ω. Let us denote by f N (y) = N n=1 f n ϕ n (y) the approximation of the function f of the order N .
The uniform convergence of the series is due to majorization by decreasing sequence e −λnT 0 .
In summary, Proposition 1 provides a powerful computational technique with the potential to be applied in a wide range of finance applications due to the fact that the NSBF representation can be used as a simple and efficient numerical method. Furthermore, the proposed novel representation is applicable to a large class of option pricing models and it represents not only the price but also the entire value function. This feature allows us to view the behavior of the option price under different initial values for the asset (i.e., to construct the value surface as will be shown in Figure  1).

The analytical representation of 'Greeks'
Since Proposition 1 presents an analytical representation of the value function, we are then able to offer a similar representation for its derivatives, commonly known as 'Greeks' in the option pricing literature. This should be a useful computational tool for both academics and practitioners, since numerical differentiation is known to be problematic in this kind of problems.

Vega
For the calculation of the Vega, we assume that the instantaneous volatility σ is differentiable and σ ′ (y 0 ) = 0. Then by the application of the chain rule and the derivative of the inverse function theorem, we have For the constant σ we cannot apply this formula. 9

Theta
The direct differentiation with respect to t of (12) provides us with a formula for the Theta As in the case of the Delta, it is necessary to assume the continuity of ∂v ∂t at (y 0 , 0).

Recurrence formulas for the coefficients α n (y) and β n (y)
For the (efficient and robust) computation of the coefficient functions α n (y) and β n (y) it is convenient to use recurrence formulas borrowed from [24]. These formulas increase the robustness of the calculations by solving the numerical issue in (15) related to the fast growth of the Legendre coefficients.
We first introduce A n (y) = l n (y) α n (y) and B n (y) = l n (y) β n (y) .
A n (y) = 2n + 1 2n − 3 l 2 (y) A n−2 (y) + (2n − 1) g (y)θ n (y) and B n (y) = 2n + 1 2n − 3 l 2 (y)B n−2 (y) + 2(2n − 1) p(y) w(y) g ′ (y)ρ(y) + g(y)ρ ′ (y) θ n (y) The initial values A 0 , A 1 , B 0 and B 1 can be calculated from There is a useful practical test for the verification of the coefficients α n and β n -its details may be consulted in [24,Equations 7 and ∞ m=0 β m (y) = l(y) q(y) 4ρ(y)w(y) The relations (25) - (28) can also be used as an indicator for the optimal choice of the number K of coefficients in the truncated series (12) and (13) to include in computations, monitoring the differences between the right sides of equations (25) -(28) and the partial sums of these. Let also k 0 be the index ofỹ (i.e. y k 0 =ỹ). Hence, we can set α n (y) = 0 for all n < k 0 . They are larger only due to the numerical error. A similar construction can be performed for the coefficients β n (y).

Implementation of the pricing algorithm
For the sake of completeness and to better describe important details of our pricing methodology, let us now provide the conceptual steps for implementing our algorithm: 1. Compute the coefficients p, q and w of the associated Sturm-Liouville problem using (7). 2. Create or choose an indefinite integration scheme. In this paper, we have used the Newton-Cotes six point integration rule-see [22] for discussions on the use of other possible methods.
3. Construct or find any non-vanishing solution g to equation (16) that satisfies (17). In our implementation, we have used the SPPS method presented in [23]. For example, if q (y) ≡ 0 (as in the case of the standard CEV model), we can choose g (y) = 1 ρ(L) as a particular solution. 4. Construct the formal power Φ (1) using Definition 1, compute the constanth and the functions G 1 (y) and G 2 (y) using formulas (24), (22) and (23), respectively. 5. Compute recursively the coefficients A m (y) and B m (y) using the representation highlighted in Section 5. 6. Compute coefficients α m (y) and β m (y) using equations (21) and verify them using relations (25) - (28) and Remark 4. We notice that this procedure can incorporate a test for estimating an optimal M (truncation parameter for the series (13) and for the second sum in the series (12)) to be used. 7. Find the eigenvalues λ n = ω 2 n from equation (14). Note that the values of the spherical Bessel functions j 2m+1 for varying indices m = 0, 1, . . . , M at the same point x can be computed efficiently using backward-recursive formula, see [1,Equation 10.
8. Construct the eigenfunctions of the problem (9) given by (13). 9. Decompose function f into the Fourier series (11) using the eigenfunctions ϕ n . 10. Construct the function v through a truncated expression (8). By N we denote the number of terms in the truncated series. 11. Calculate the Greeks via expressions (18) - (20).

Remark 5.
Notice that the proposed algorithm can be significantly simplified if we are interested only in the price of the option v (y 0 , 0). If this is the case, then in steps 4, 5 and 6 we only need terms relative to A n and α n (y). Moreover, after calculating f n we do not need to keep the eigenfunctions, but only values at the point y 0 . Further, at step 10, we construct only v (y 0 , t) and step 11 is not necessary.

Computational experiments and particular examples
In this section we apply the algorithm described in the previous section to the EJDCEV model, whose details will be described next. For illustrative purposes, we have separated the examples in two different time horizons, the medium (six months) and the short (one day). This particular choice will highlight the eigendata needed for the accurate computations and the sensibility of the model to the chosen parameters.
We note that for the regular Sturm-Liouville problems that we are considering in this paper, the asymptotics for the eigenvalue growth is of the order of n 2 , e.g. [38,Section 2.13]. In the case of the long horizon the exponential term e −λn(T −t) , with T − t of order 1 2 , decays rapidly and the representation (12) converges quickly. Hence, few eigenvalues and eigenfunctions are needed to secure a good approximation. For the short horizon case, with T − t of order 1 360 , that is analyzed in the second set of numerical experiments, we need more eigenvalues to have an accurate approximation for the option value. We further note that the NSBF method calculates the required eigendata with the same accuracy. This NSBF's important property, of not loosing accuracy for the highest order eigenvalues, makes it an exciting tool for applications to problems requiring large sets of eigenvalues.
Another computational advantage of our method is that there is no need in any two-dimensional grid for computation. The formulas for the steps 1-9 are one-dimensional. In order to make the integration errors negligible and to concentrate mainly on the numerical performance we have used an overwhelming number of mesh points (10000) on the interval [L, U ] to represent all the functions involved, moreover, even 3000 mesh points produced close results. Once all the coefficients are obtained, the computation of the value function and Delta may be performed only for the arguments (y, t) ∈ [L, U ] × [0, T ] required by application. E.g., option price can be obtained as the value of v at one point (y 0 , 0); the value surface requires calculation of the function v on a mesh of about 101 × 101 points, etc. Even though the main purpose here is to present the ability of the algorithm to be used in a wide variety of modeling contexts and not the optimization speed, the very small computational burden that is required is remarkable.
All the calculations where done in Matlab R2015a.

The EJDCEV model
The volatility specification under the time-homogeneous version of the JDCEV model is given by with δ > 0 and β ∈ R. The drift is given by expression (3), withr > 0,q ≥ 0 and hazard rate with b > 0 and c > 0. The properties of the constructed diffusion with different parameterizations can be consulted in [7] and [29]. The nice feature of the JDCEV model is its analytical tractability, due to the special form of the assumed hazard rate h 1 (y). The advantage of the NSBF representation is that it allows us to consider different default intensity specifications without any additional effort. Following [6], we have also considered a default intensity specification guaranteeing a positive relationship between the default probability and volatility. Hence, in our variant of the JDCEV model, that we name as the EJDCEV, the default intensity is assumed to be dependent of the constant parameter γ ≥ 0 as h (y) := h alt 1 (y) = b + cσ γ (y) .
It is important to point out that we are not restricted to functions of the form h alt 1 ; we can choose any positive continuously differentiable function. This feature allows the possibility of obtaining a wide alternative of default rates when calibrating the model to market prices.

Medium horizon (6 month)
For this set-up, we adopt the parameter configuration considered in [14, Table 2, Panel C], that is y 0 = 100, L = 90, U = 120, T = 0.5,r = 0.1,q = 0, b = 0.02, c = 0.5, and σ 0 = 0.25. As usual, the scale parameter δ is calculated, for each β value 11 , through the relation while keeping the initial instantaneous volatility σ 0 = 0.25. We notice that the determination of the spectral parameters in step 7 was performed by interpolation with a grid of equally distributed 100 points on the interval (0, 15) for the practical purposes and the grid of 1000 points on the interval (0, 50) for the construction of the illustration Table 2 and graphs. The practical reasoning is to cut out the eigenvalues λ n > 15 2 due to the term e λn(T −t) in our formulas, this indirectly sets the parameter N somewhere around 6, as can be observed in Table 2. Figure 1 illustrates the value function under the JDCEV model using the aforementioned parameters coupled with K = 100, β = −1 and γ = 2. Using the same set of parameters, Figure 2 shows the detail of the approximation of the function f (y) = (y − K) + at the boundary for t = T . It is possible to observe a sharp decline at the boundary U , this is the illustration of the Remark 2.    Table 1 shows the prices of European-style DBKO call and put options and the corresponding Greeks under the EJDCEV model for different moneyness levels with K ∈ {95, 100, 105}, β ∈ {0.5, 0.0, −1.0, −2.0} and γ ∈ {0, 1, 2}. We further note that the six cases with γ = 2 and β ∈ {−1.0, −2.0} originate the values of DBKO put options shown in [14, Table 2, Panel C]. A direct comparison reveals that the results are exactly the same (rounded to four decimal places of accuracy), which gives further evidence on the robustness of our algorithm. More importantly, this also allows us to test our methodology under a larger set of volatility and default intensity specifications that until now were not possible to be tackled in the literature.

Short horizon (1 day)
In the case of the short horizon, the time is of the order 1 360 and, hence, the term e −λn(T −t) decays much slower as n grows. For this case, in step 7, we have used 1000 points for the interval (0, 100). Some values are presented in Table 2. In order to illustrate the convergence and the necessity of calculating accurately a significant number of eigenvalues, we introduce the contribution of the partial sum from equation (8) We observe, in Table 3, that the value of the parameter β under the short horizon does not have much influence on the price. However, the γ parameter associated with default intensity is relevant. It is important to note that although the prices of the options for different β values differ slightly, the corresponding Sturm-Liouville problems are very different. This can be observed in Tables 2  and 4. The observation of the Table 4 induces the choice of N around 45. This table shows the eigenvalues for different γ and β parameters, with y0 = 100, K = 100, L = 90, U = 120, T = 0.5,r = 0.1,q = 0, b = 0.02, c = 0.5, and σ0 = 0.25.   This table shows the value of the contribution defined in equation (32) for one-day to maturity European-style DBKO call options for different γ and β parameters, with y0 = 100, K = 100, L = 90, U = 120, T = 1/360,r = 0.1,q = 0, b = 0.02, c = 0.5, and σ0 = 0.25.

Concluding remarks and future research
This paper provides a new methodology for pricing (and hedging) European-style DBKO options via the application of the NSBF decomposition of the Sturm-Liouville equation associated to the corresponding boundary value problem. The illustration of the method was done through the EJDCEV model. The modeling techniques applied in this paper open several avenues for future research. For instance, it should be possible to apply the NSBF decomposition and similar constructions to other singular problems that naturally appear in many financial applications, e.g. plain-vanilla options (unbounded domains), default cases (singularities in the coefficients) and others. It would also be interesting to apply the method to calibrate a parametric curve of parameters to real market data. Finally, it has also the potential to be applied to stopping time problems and related subjects. scholarship granted by the Mexican Government via the Ministry of Foreign Affairs which gave him the opportunity to develop this work during his stay in the CINVESTAV, Mexico.

Appendix A. Pricing rebates
Consider the case of the call option with rebate R > 0. The upper boundary condition for the problem (6) changes to v (U, t) = R.