Abstract

In this paper, a new approach is proposed to construct willow tree (WT) for generalized hyperbolic (GH) Lévy processes. There are two advantages of our proposed approach compared to the classical WT methods. Firstly, it avoids the moments matching from Johnson curve in the known WT construction. Secondly, the error of European option pricing is only determined by time partition under some conditions. Since the moments of Lévy measure are removed from this algorithm, our approach improves the stability and accuracy of WT in option pricing. Numerical experiments support our claims. Moreover, the new approach can be extended to other Lévy processes if their characteristic functions are expressed by explicit forms.

1. Introduction

To provide more modeling tools on jump types, the Lévy process is a generalization of the diffusion processes by allowing infinite jump activity. The standard form of the Lévy process assumes stationary increments, hence resulting in nice analytic tractability. The assumption of Lévy processes makes it a good choice for pricing equity derivatives. Generalized hyperbolic (GH, see [1]) process is a class of Lévy processes with wide application in financial field. Variance gamma (VG, see [2, 3]) and normal inverse gamma (NIG, see [47]) processes are two special cases of the GH model. VG and NIG models can be obtained from a randomized time-changed clock. The time-changed processes and general Lévy processes exhibit stochastic volatility such that they become capable of capturing volatility smile, smile skew, and term structure of the smile. Also, the hyperbolic (HYP) distribution and the Gauss normal (GN) distribution are viewed as subclasses of the GH model.

Currently, there are five popular numerical methods to price options under Lévy processes, binomial tree methods (BTMs, see [3, 8, 9]), finite difference methods (FDMs, see [1013]), Monte Carlo methods (MCM, see [14]), FFT-based transformation methods (see [1518]), and cosine-willow tree methods (see [19]). For BTMs and WTMs, computing transition probabilities (TPs) is an insurmountable barrier. Although a few literature studies discuss BTMs and WTMs for Lévy processes, the accuracy and computational efficiency of them are needed to be enhanced. To use FDMs to valuate options, corresponding partial integral-differential equations (PIDEs) should be established. Generally, PIDEs governed by Lévy processes are very complicated such that PIDEs are limited in several options. The Monte Carlo method is straightforward, but it is quite time-consuming to generate random samples for the Lévy process. The FFT-based transformation method is the most popular one in option pricing, such as the COS method (see [16]) and PROJ method (see [17]). The COS method employs a cosine series expansion on the risk-neutral return density and estimates the European option price based on the numerical integration on . However, it is hard to determine a proper finite interval to truncate for the integration (see [15, 18, 20]) and is hard to be extended to path-dependent options. The PROJ method overcomes these shortcomings and is extendable to Asian options, variance swaps, and American options but its extendability is still limited compared to the Cosine-willow tree method.

The key of the WT method for option pricing is to construct willow tree structure (see [13]). In this paper, we propose WT algorithms for GH Lévy processes. The main contributions of this paper are fourfolds.(1)An numerical algorithm is proposed to compute the probability density functions (PDFs) and cumulative distribution probabilities (CDPs) from the characteristic functions (CDs) of GH processes. This algorithm is unified and suitable for all GH subclass models.(2)In FFT-based transformation methods (see [15, 18, 20, 21]), it is hard to find a proper finite interval to truncate the integration over . While in the WT method, we propose an adaptive integration method in which the appropriate integral interval is automatically found. In determining , the WT algorithm does not consume too much extra computational effort.(3)By setting appropriate discrete stock prices at each time and calculating transform probabilities from to , WT structure is constructed. Unlike Ma et al. [19], it is not needed to estimate the -order moments of the GH Lévy model when selecting nodes at each time . The determination of nodes at each time only relies on PDFs or CDFs of GH processes, which makes our WT programming run fast than those developed by Ma et al. (see [19]).(4)The convergence rate of European option on the WT structure is proved. Numerical experiments show that American options computed by our WT algorithm are also convergent for time partition and the underlying partition number .

The remaining parts of this paper are arranged as follows. Basic conceptions are reviewed in Section 2. Calculation of PDFs and CDFs is illustrated in Section 3. The convergence analysis for European options is discussed in Section 4. Numerical examples of GH Lévy processes are carried out in Section 5. Some conclusions and remarks are given in the final section.

2. GH Lévy Processes

In option pricing with non-Gaussian processes, the asset price process is defined as an exponential Lévy process , i.e.,where is the risk-free interest rate and is a martingal adjustment parameter under risk-neutral measure . Under the GH Lévy processes (1), option pricing becomes more complicated than the classical BS model. For FDMs, partial integral-differential equations (PIDEs) are needed, and then numerical schemes should be designed to solve these PIDEs. For BTMs and WTMs, calculating transform probabilities (TPs) cannot be avoided. Formulating PIDEs, designing numerical schemes for PIDEs, and calculating TPs for BTMs and WTMs are not easy tasks.

We consider the generalized hyperbolic (GH) model (see [1, 22, 23]) whose characteristic function with five parameters is defined by the following equation:where is the order-modified Bessel function of the second kind. The density function of the GH model can be derived as follows:

In CFs (2), and the other parameters satisfy the following constraints:

The GH distribution embeds various distributions under special choices of the parameters. The parameter determines the subclass of the GH distribution. When , the GH distribution reduces to the hyperbolic distribution (HYP) whose logarithm of density is a hyperbolic. In addition, when and , the GH distribution reduces to the Gauss normal (GN) distribution. Furthermore, when and , the GH distribution becomes the variance gamma (VG) distribution; when , it becomes the normal inverse Gaussian (NIG) distribution. Table 1 gives the four different categories of the GH model.

Figure 1 plots the graphs of with different values of parameters. From the Figure, we see the shape of tail, skewness, and kurtosis of the GH distribution which are controlled by , and , respectively.

Because the GH law is infinitely divisible, one can construct a GH Lévy process whose distribution at fixed time has characteristic function . The characteristic function is described by the following equation (see [22]):where , is the density function of GH process at time , and is defined as (2). The function is the part in which the oscillating factor is removed. The characteristic exponent of GH process is as follows:

Since a Lévy process is an infinitely divisible distribution, the Lévy-Khintchine theorem (see Theorem 8.1 in [24]) for an infinitely divisible distribution can be applied to establish the characteristic exponent of a Lévy process . admits the Lévy–Khintchine representation, i.e.,

The Lévy measure is defined on the real domain excluding zero with . The triplet is called the Lévy characteristic of , where is the constant drift, is the constant volatility of the continuous component, and is the Lévy measure that represents the expected number of jumps per unit time. In GH model (1), the parameter is set as zero. The Lévy measure of the GH model has no explicit analytic form and it can be expressed only in terms of integrals (see [25]). To satisfy the martingale condition in SDE (1), the adjustment parameter is chosen as .

The density function can be restored from characteristic function by numerical inverse integral. Once the characteristic function (5) is given, the density function is expressed aswith definition in (5). We note that the integral factor is the part with high frequency oscillation, whereas is decaying quickly as . This observation is very important for computing the density function , numerically.

3. Willow Tree Algorithm

3.1. Willow Tree Structure

A willow tree for the GH model is represented by , where is discrete nodes of Lévy processes at time , is underlying prices, and is the transition probability. As shown in Figure 2, there are two main stages to construct a willow tree: (i) selecting the discrete tree nodes, ( and ), for at each time and (ii) determining the transition probability, , from at to at on a discrete time points with and , .

Firstly, the discrete pairs are selected to approximate the distribution of at time . The cumulative distribution functions (CDFs) of at are computed by

After setting discrete probabilities,we can numerically determine by solving equation (9) if the integral can be computed efficiently and accurately. Very different from the classical WT method (see [19]), the determination of nodes avoids computing -order moments of . This modification makes our algorithms more easily and widely available for Lévy processes.

Secondly, the transition probability between two consecutive discrete times and is computed aswith increments of Lévy processes. Increments and are defined as

In computation (11), and , and an exception is that with .

There are two most important aspects in willow tree construction. One is to calculate PDF and then generate willow tree nodes and another is to compute transition probability via for given . Once the PDFs and are computed, the willow tree algorithm can be applied in some types of option pricing, such as European, American, Asian options, and so on.

3.2. Numerical Computing of PDFs and CDFs for GH Processes

Firstly, we numerically compute the PDFs defined by (8). There is a challenge for the high frequency oscillation of integral kernel when the value of is large. We consider the following truncation:

In the above expression, is a sufficiently large negative number and is a large enough positive constant. Let partition nodes for with . Then, can be numerically integrated (NI), i.e.,with notation

Remark 1. The numerical formula (14) plays a key role for both and . (i) Let , we can calculate CDFs at each time . (ii) Nodes can be determined by solving for (see expression (9)). (iii) Transition probabilities (see formula (11) can be calculated by numerical approximation with . (iv) Using formula (14), we can expand the range of until the result does not change much, which can be seen as an adaptive algorithm.

Remark 2. Since the absolute values of characteristic function are decreasing to zero as tends to , the nodes can be selected as nonuniform, for example,with for . Using nonuniform nodes , we can achieve more accurate with less nodes.
Next, we compute transition probabilities with definition (11). Taking discrete values with sufficiently large negative number and sufficiently small , CDFs are approximated byfor and . The transition probability is approximated bywith and being denoted by expression (12). Using (9) and (14)–(18), Algorithm 1 describes the details of constructing willow tree for GH Lévy processes.

%% Input parameters
Step 1. Set GH-model parameters ().
%% Compute probability density function
Step 2. Set willow tree parameters () and computational parameters . represents the number of time discretization and represents the number of -nodes at each time . is the number in NI formula (14) for PDFs .
%% Compute probability density function for and .
Step 3. Set discrete space nodes with for .
Step 4. Compute and for using formulas (15);
Step 5. According to (14), compute PDFs on given points .
  %% Generate WT nodes and compute transition probabilities.
Step 6. According to (17), compute at each time . By solving (3.1), determine nodes for and .
Step 7. Compute discrete underlying for and .
Step 8. Based on (18), compute transition probability for and .
Step 9. Willow tree for the GH model is constructed.
3.3. European and American Options

After the willow tree is constructed, European and American options can be valued on WT backward. Other options (Asian, Lookback, and so on) also can be computed from WT and we omit the details.

Define by the option values at time with underlying . Option values at with parameters can be computed from willow tree as follows:where payoff function for call options, whereas for put options.

American option value at time zero with parameters can be determined backward on willow tree.

American options via willow tree are described in Algorithm 2.

%% Input parameters and prepare willow tree
Step 1. Set willow tree parameters () and computational parameters . represents the number of time discretization and represents the number of -nodes at each time . is the number in NI formula (14) for PDF .
Step 2. Generate willow tree by Algorithm 1.
%% Compute options backward from time to
Step 3. Compute payoff function at time , i.e., for , where for call options, whereas for put options.
Step 4. For to 1
—— Compute European values for
—— Compute exercise values for
—— Take American values for
EndFor
%% Output American option
Step 5. Calculate option value .
3.4. A Simple Method to Determine Nodes at Each Time

To generate nodes at each time , it is much time-consuming according to Step 6 in Algorithm 1. We give a simple algorithm to generate nodes only depending on CDFs (see (17)). Assume as given, the nodes are generated as follows. Find such thatfor , where is the CDFs of the GH model (defined in (17)) and . If the discrete probability distribution is given with probabilities and corresponding increments , then can be obtained by interpolation:

Now, we describe the method of generating nodes as in Algorithm 3. Compared with the method proposed by this algorithm avoids to compute the -moments of . The modification of our algorithm saves much CPU time in willow tree construction.

Step 1. Give the first node and at time zero.
%% generate CDFs on nodes
Step 2. Based on (17), generate discrete distribution for . Here, is a sufficient small number, is a sufficient larger number and sufficient small.
%% generate for .
Step 3. For , generate by interpolation,
, .
  END
Step 4. Compute discrete underlying and then the nodes of willow tree are generated.

4. Convergence Analysis

Errors of option pricing from willow tree have two aspects: the errors (see (14)) coming from the calculation for PDFs (or CDFs) and the errors from backward computation on willow tree (see (19) and (20)). We give some detailed analysis in this section.

4.1. Errors of PDFs and CDFs

We know that the absolute value of characteristic function is decreasing to zero as tends to . So, the error can be estimated as follows.

Theorem 3. Assume that the characteristic function satisfies for , thenwhere PDFs and estimated PDFs are defined by (14) and is the number of nodes .

Proof. It is obvious thatwhere and for . From error estimation of composite trapezoidal rule (see Richard and Dougals Faires [26]), we havewhich is the result of (23).

Remark 4. Given , to ensure the error be small enough, , , and should be taken as large as enough. For given , a simple choice is to select and such thatand then find such thatAs an example, Table 2 gives the choices of for different values of and fixed .

4.2. Convergence of Willow Tree for European Option

Denoted by the conditional moments of increment with given , i.e.,where is the PDFs of the increment of at given . For convenience, we give an assumption as follows.

Assumption 5. There exists a positive number such thatfor and . Here, and are denoted by expression (12).

Lemma 6. Assume as a GH Lévy process, then and have the following results:with .

Proof. The proof can be seen in Cont and Tankov [27].

Lemma 7. Given following GH Lévy processes, the conditional moment of increment has estimation for .

Proof. Given the characteristic exponent of and characteristic function , it implies thatSince , we havewhere is the cumulants of . Thus, we obtain for .
The following lemma gives an estimation of , the error between moments, and their discrete approximations.

Lemma 8. Given following GH Lévy processes, under Assumption 5, i.e., is bounded by with constant , the errors between and its discrete approximations with can be estimated by

Proof. Using the following integral operator:the errors for moments can be written asThanks to Cauchy–Schwarz inequality and Lemma 6, we havewhich is the result of (33).

The following theorem gives the error estimation of the willow tree algorithm for European options.

Theorem 9. Given the asset price governed by the exponential Lévy model , on discrete times with and discrete values being generated by (9). If satisfies conditions in Assumption 5, the error between the true value of the European option and the computed value by the backward induction (19) in the willow tree is when is in .

Proof. It is known that the European option with is the solution of the following partial integrodifferential equation (PIDE):with the terminal condition , ( for call option or for put option), and being the Lévy measure (see Lévy representation Theorem 2.7 in [7]). Given willow tree, the European option with is computed by the backward induction as in (19), i.e.,Expanding at by the Taylor series, we havewhere and . Thus, the backward induction (38) can be written asFrom Lemma 8, the discrete approximation of the first- and second-order moments of the GH process can be estimated asFrom (40)–(43), we haveOn the other hand, using the Taylor expansion of at , we havewith . Applying the properties of Lévy measure, it is obtained thatTherefore, using (45) and (46), it is yielded thatCombining (44) and (47), we haveat all nodes . Comparing (48) with (37), we see any European option value computed by the backward induction (38) satisfies PIDE (37) with error term . When is in , the error term is emerged as .

Remark 10. To obtain option values with good accuracy, the number of should be taken as . Since European options are path-independent, the number of should be selected as small as possible, whereas the number should be taken as large enough since should be taken as large as possible for American options.

5. Numerical Examples

To test the performance of the willow tree (WT) method for option valuating, we consider four classes of choices (GH, HYP, NIG, and VG models) with parameters being listed in Table 3. Parameters of risk-free interest, initial stock price, and maturity time are set as . are set in the WT algorithm for European options, whereas for American options. In Monte Carlo (MC) simulation, are set with representing the number of time partition and represents the number of simulated paths. In numerical formula (14), we set numerical partition for variable of characteristic functions.

All experiments are carried out by MATLAB R2012b running on a machine with Intel(R) Core (TM) i7-8550U CPU @ 1.80 GHz, 8 GB RAM under Windows 10.

Figure 3 plots probability density functions (PDFs) of GH, HYP, NIG, and VG models with . The figure shows that the PDFs computed from numerical formula (14) are very close to explicit expression (3). So, we believe that PDFs computed from (14) with time are also accurate enough. Figure 4 plots the shape of cumulative distribution functions (CDFs) with different numbers of , from which we see the computed CDFs are convergent as becoming larger. Figure 5 plots the trajectories of nodes and for the NIG model.

Option values computed from WT algorithms, MC simulations, and analytical solutions (labeled by “ECA”) are listed in Table 4. Figure 6 plots values of European and American options, from which we see WT solutions are very close to those obtained by MC simulation. We see that the errors between WT solution and analytical solutions (or MC solutions) are about . In European options computation, the CPU time consumed from WT is less than 3 seconds whereas more than 18 s for MC simulation. In American option computation, the CPU time consumed from WT is less than 10 s, whereas more than 160 s for the MC method. The results in Table 4 illustrate the effectiveness of the proposed WT method. The analytical formula of European options under NIG and VG processes can be seen in literatures [2, 6, 28, 29]. Those pricing formulas are also listed in Appendix A and Appendix B.

To test the convergence of the willow tree method with respect to the number of space nodes and time partition number , some experiments are carried out. Figure 7 plots the errors for different parameters with fixed . Figure 8 plots the errors with different and corresponding , from which we see the errors are decreasing as (and so ) increases. Table 5 lists the numerical convergent rates with respect to the values of . These results in Figures 7 and 8 and Table 5 support the theoretical conclusion in Theorem 9.

6. Conclusions

In this paper, a unified and robust approach is proposed to construct the willow tree structure for GH Lévy processes. There are two advantages of our proposed approach compared to that in [19]. First, it avoids the moment matching failure by the Johnson curve under some circumstances in the willow tree construction. Second, the error of European call option pricing is only determined by . The fifth-moment term of Lévy measure is removed from the error bound, so our approach improves the stability and accuracy of the willow tree in option pricing. Numerical experiments support our claims.

Moreover, we believe the proposed willow tree method can be extended to other option models, such as variance and volatility swaps with stochastic volatility and stochastic volatility model with regime switching stochastic mean reversion (see, e.g., [3033]). We will discuss those models in the future.

Appendix

A. Analytical Solution under VG Model

When the risk-neutral dynamics of the stock price is given by the VG process (for risk-neutral parameters ), the European call option price on a stock is (see page 88 in [2]):where ,and is defined in terms of the modified Bessel function of the second kind.

B. Analytical Solution under NIG Model

When the risk-neutral dynamics of the stock price is given by NIG, the European call option price is (see Theorem 2.1, Theorem 2.2, and Corollary 2.1 in Ivanov [6])with definitionsand

In the above formula, is the MacDonald function, is the beta function, and is the degenerate Appell hypergeometric function.

Data Availability

The authors confirm that the data supporting the findings of this study are available within the article.

Conflicts of Interest

The authors declare that there are no conflicts of interest regarding the publication of this paper.

Acknowledgments

This work was supported by the National Natural Science Foundation of China (grant no. 12171409) and Key Project of Hunan Education Department (grant no. 238945).