Option Pricing by Willow Tree Method for Generalized Hyperbolic L´evy Processes

In this paper, a new approach is proposed to construct willow tree (WT) for generalized hyperbolic (GH) L´evy processes. Tere are two advantages of our proposed approach compared to the classical WTmethods. Firstly, it avoids the moments matching from Johnson curve in the known WTconstruction. Secondly, the error of European option pricing is only determined by time partition ∆ t under some conditions. Since the moments of L´evy 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´evy processes if their characteristic functions are expressed by explicit forms.


Introduction
To provide more modeling tools on jump types, the Lévy process is a generalization of the difusion processes by allowing infnite jump activity.Te standard form of the Lévy process assumes stationary increments, hence resulting in nice analytic tractability.Te 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 fnancial feld.Variance gamma (VG, see [2,3]) and normal inverse gamma (NIG, see [4][5][6][7]) processes are two special cases of the GH model.VG and NIG models can be obtained from a randomized time-changed clock.Te 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 fve popular numerical methods to price options under Lévy processes, binomial tree methods (BTMs, see [3,8,9]), fnite diference methods (FDMs, see [10][11][12][13]), Monte Carlo methods (MCM, see [14]), FFTbased transformation methods (see [15][16][17][18]), 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 efciency of them are needed to be enhanced.To use FDMs to valuate options, corresponding partial integral-diferential equations (PIDEs) should be established.Generally, PIDEs governed by Lévy processes are very complicated such that PIDEs are limited in several options.Te Monte Carlo method is straightforward, but it is quite time-consuming to generate random samples for the Lévy process.Te FFT-based transformation method is the most popular one in option pricing, such as the COS method (see [16]) and PROJ method (see [17]).Te COS method employs a cosine series expansion on the riskneutral return density and estimates the European option price based on the numerical integration on (− ∞, +∞).However, it is hard to determine a proper fnite interval [a, b] to truncate (− ∞, +∞) for the integration (see [15,18,20]) and is hard to be extended to path-dependent options.Te 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.
Te 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.Te 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.Tis algorithm is unifed and suitable for all GH subclass models.(2) In FFT-based transformation methods (see [15,18,20,21]), it is hard to fnd a proper fnite interval [a, b] to truncate the integration over (− ∞, +∞).While in the WT method, we propose an adaptive integration method in which the appropriate integral interval [a, b] is automatically found. In

GH Lévy Processes
In option pricing with non-Gaussian processes, the asset price process S t is defned as an exponential Lévy process X t , i.e., where r is the risk-free interest rate and ω GH is a martingal adjustment parameter under risk-neutral measure Q.Under the GH Lévy processes (1), option pricing becomes more complicated than the classical BS model.For FDMs, partial integral-diferential 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 fve parameters (α, β, δ, μ, λ) is defned by the following equation: where K λ (•) is the λ th order-modifed Bessel function of the second kind.Te density function of the GH model ρ GH (x) can be derived as follows: In CFs (2), μ ∈ R and the other parameters satisfy the following constraints: Te GH distribution embeds various distributions under special choices of the parameters.Te parameter λ determines the subclass of the GH distribution.When λ � 1, the GH distribution reduces to the hyperbolic distribution (HYP) whose logarithm of density is a hyperbolic.In addition, when δ ⟶ ∞ and δ/α ⟶ σ 2 , the GH distribution reduces to the Gauss normal (GN) distribution.Furthermore, when δ � 0 and μ � 0, the GH distribution becomes the variance gamma (VG) distribution; when λ � − 1/2, it becomes the normal inverse Gaussian (NIG) distribution.Table 1 gives the four diferent categories of the GH model.
Figure 1 plots the graphs of ρ GH (x) with diferent 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 infnitely divisible, one can construct a GH Lévy process X t whose distribution at fxed time t has characteristic function ϕ X t .Te characteristic function is described by the following equation (see [22]): ( where is the density function of GH process at time t, and ϕ GH (u) is defned as (2).Te function  ϕ GH (u) � ϕ GH (u)e − iμu is the part in which the oscillating factor e iμu is removed.Te characteristic exponent of GH process is as follows: Since a Lévy process is an infnitely divisible distribution, the Lévy-Khintchine theorem (see Teorem 8.1 in [24]) for an infnitely divisible distribution can be applied to establish the characteristic exponent ψ X (u) of a Lévy process X t .ψ X (u) admits the Lévy-Khintchine representation, i.e., 1 − e iux + iux1 |x|<1  Π(dx).(7) Te Lévy measure Π(dx) is defned on the real domain excluding zero with Π( 0 { }) � 0. Te triplet (μ, σ 2 , Π) is called the Lévy characteristic of X t , where μ ∈ R is the constant drift, σ > 0 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.Te 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 Te density function ρ GH t (x) can be restored from characteristic function ϕ X t (u) by numerical inverse integral.Once the characteristic function ( 5) is given, the density function ρ GH t (x) is expressed as with defnition  ϕ X (u) in (5).We note that the integral factor e − iu(x− μt) is the part with high frequency oscillation, whereas  ϕ X (u) is decaying quickly as |u| ⟶ ∞.Tis observation is very important for computing the density function ρ GH t (x), numerically.is underlying prices, and p (n)  ij 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,

Willow Tree Algorithm
Firstly, the discrete pairs (X n i , q i )  (i � 1, 2, • • • , m; n � 1, 2, . . ., N) are selected to approximate the distribution of X t at time t n .Te cumulative distribution functions (CDFs) of X t � X n i at t n are computed by After setting discrete probabilities, we can numerically determine X n i by solving equation ( 9) if the integral can be computed efciently and accurately.Very diferent from the classical WT method (see [19]), the determination of nodes X n i avoids computing k-order moments of X t .Tis modifcation makes our algorithms more easily and widely available for Lévy processes.
Table 1: Special cases of GH distribution.

Parameters Type of distribution
Tere are two most important aspects in willow tree construction.One is to calculate PDF ρ GH t n (x) and then generate willow tree nodes X n i   and another is to compute transition probability p (n)  ij via ρ GH ∆t (x) for given ∆t.Once the PDFs ρ GH t n (x) and ρ GH ∆t (x) are computed, the willow tree algorithm can be applied in some types of option pricing, such as European, American, Asian options, and so on.

Numerical Computing of PDFs and CDFs for GH Processes.
Firstly, we numerically compute the PDFs ρ GH t (x) defned by (8).Tere is a challenge for the high frequency oscillation of integral kernel e − iu(x− μt) when the value of (x − μt) is large.We consider the following truncation: In the above expression, A is a sufciently large negative number and B is a large enough positive constant.Let partition nodes  4 Journal of Mathematics with notation ). (iii) Transition probabilities p (n)  ij (see formula ( 11) can be calculated by numerical approximation with t � ∆t.(iv) Using formula (14), we can expand the range of [A, B] until the result ρ GH t (x) does not change much, which can be seen as an adaptive algorithm.
Remark 2. Since the absolute values of characteristic function ϕ X (u) are decreasing to zero as u tends to ± ∞, the nodes u k   can be selected as nonuniform, for example, Using nonuniform nodes u k  , we can achieve more accurate  ρ GH t (x) with less nodes.
i with sufciently large negative number x 0 and sufciently small ∆x ≔ � max 1≤j≤i [x j − x j− 1 ], CDFs F n (X n i ) are approximated by with being denoted by expression (12).Using ( 9) and ( 14)-( 18), Algorithm 1 describes the details of constructing willow tree X n i , S n i , p (n) ij   for GH Lévy processes.

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 WTand we omit the details.
Defne by V n i the option values at time t n with underlying S n i � S 0 e (r+ω GH )t n +X n i .Option values V at t � 0 with parameters (S 0 , K, T, r) can be computed from willow tree as follows: where payof function f(ξ) � (ξ − K) + for call options, whereas f(x) � (K − ξ) + for put options.American option value V at time zero with parameters (S 0 , K, T, r) can be determined backward on willow tree.
American options via willow tree are described in Algorithm 2.

A Simple Method to Determine m Nodes at Each Time t n .
To generate m nodes X n i at each time t n , it is much timeconsuming according to Step 6 in Algorithm 1.We give a simple algorithm to generate m nodes only depending on CDFs F n (•) (see (17)).Assume X n i   as given, the nodes X n+1 i are generated as follows.Find X n+1 i such that Journal of Mathematics i can be obtained by interpolation: Now, we describe the method of generating nodes X n i , S n i   as in Algorithm 3. Compared with the method proposed by this algorithm avoids to compute the k-moments of X t .Te modifcation of our algorithm saves much CPU time in willow tree construction.

Convergence Analysis
Errors of option pricing from willow tree have two aspects: the errors ε(ρ) � ρ GH t (x) −  ρ GH t (x) (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.

Errors of PDFs and CDFs.
We know that the absolute value of characteristic function ϕ X (u) is decreasing to zero as u tends to ± ∞.So, the error ε(ρ) can be estimated as follows.

Theorem 3. Assume that the characteristic function satisfes |ϕ
where PDFs ρ GH t (x) and estimated PDFs  ρ GH t (x) are defned by (14) for the GH model is constructed.
%% Input parameters and prepare willow tree Step 1. Set willow tree parameters (N, m) and computational parameters M. N represents the number of time discretization and m represents the number of X-nodes at each time t n .M is the number in NI formula (14) for PDF ρ GH t (x).
Step 2. Generate willow tree X n i , ALGORITHM 2: American option pricing on willow tree.

6
Journal of Mathematics where From error estimation of composite trapezoidal rule (see Richard and Dougals Faires [26]), we have which is the result of (23).

□
Remark 4. Given x ∈ R, to ensure the error ε(ρ) be small enough, |A|, B, and M should be taken as large as enough.
For given ε(ρ), a simple choice is to select |A| and B such that min and then fnd M such that As an example, Table 2 gives the choices of A, B, M for diferent values of x and fxed ε(ρ) � 10 − 6 .

Convergence of Willow Tree for European Option.
Denoted by E (k)  ∆t the conditional k th moments of increment ∆X � X n+1 − X n with given X n � X n i , i.e., where ρ GH ∆t (u) is the PDFs of the increment of X t at t n+1 given X n � X n i .For convenience, we give an assumption as follows.

Lemma . Assume X t as a GH Lévy process, then E[X t ] and Var[X t ] have the following results:
Step 1. Give the frst node X 0 1 � 0 and S 0 1 � S 0 at time zero.%% generate CDFs F on nodes x Step 2. Based on (17), generate discrete distribution F j � F n�1 (x j ) for j � 0, 2, • • • , L. Here, x 0 is a sufcient small number, x L is a sufcient larger number and max j�1: L |x j − x j− 1 | sufcient small.%% generate

Journal of Mathematics
Proof.Te proof can be seen in Cont and Tankov [27].□ Lemma 7. Given X t following GH Lévy processes, the k th conditional moment of increment ∆X � X n+1 − X n i has estimation E (k)  ∆t � O((∆t) k ) for k ≥ 1.
Proof.Given the characteristic exponent ψ X (u) of X t and characteristic function ϕ X Δt (u), it implies that Since where Te following lemma gives an estimation of R (k)  ∆t,m , the error between k th moments, and their discrete approximations.

□ Lemma 8. Given X t following GH Lévy processes, under Assumption 5, i.e., h is bounded by H/m with constant H, the errors between E (k)
∆t and its discrete approximations Proof.Using the following integral operator: the errors for k th moments E (k) ∆t can be written as Tanks to Cauchy-Schwarz inequality and Lemma 6, we have

8
Journal of Mathematics which is the result of (33).

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

Theorem 9. Given the asset price S t governed by the exponential Lévy model
i being generated by (9).If h satisfes conditions in Assumption 5, the error between the true value V(x, t) of the European option and the computed value V(x, t) by the backward induction (19) in the willow tree is Proof.It is known that the European option V(y, t) with y � log S(t) is the solution of the following partial integrodiferential equation (PIDE): with the terminal condition V(y, T) � f(y), (f(y) � (e y − K) + for call option or f(y) � (K − e y ) + for put option), and Π(dθ) being the Lévy measure (see Lévy representation Teorem 2.7 in [7]).Given willow tree, the European option V(y n i , t n ) with y n i � logS n i � (r + ω GH )t n + X n i is computed by the backward induction as in (19), i.e., Expanding V(y n+1 j , t n+1 ) at (y n i , t n ) by the Taylor series, we have where ∆Y n ij � y n+1 j − y n i and ξ ∈ (y n i , y n i + ∆Y n ij ).Tus, the backward induction (38) can be written as From Lemma 8, the discrete approximation of the frstand second-order moments of the GH process can be estimated as From we have On the other hand, using the Taylor expansion of V(y n i + θ, t n ) at (y n i , t n ), we have with ξ ∈ (y n i , y n i + ∆Y n ij ).Applying the properties of Lévy measure, it is obtained that 10 Journal of Mathematics Terefore, using (45) and (46), it is yielded that (47)    3.

□
Remark 10.To obtain option values with good accuracy, the number of m should be taken as O(N 3/2 ).Since European options are path-independent, the number of N should be selected as small as possible, whereas the number m should be taken as large enough since N should be taken as large as possible for American options.3.

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 r � 0.03, S 0 � 10, T � 1. (N � 1, m � 200) are set in the WT algorithm for European options, whereas (N � 100, m � 200) for American options.In Monte Carlo (MC) simulation, (N mc � 100, M mc � 5, 000, 000) are set with N mc representing the number of time partition and M mc represents the number of simulated paths.In numerical formula (14), we set numerical partition M � 50, 000 for variable u 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 t � 1. Te fgure shows that the PDFs computed from numerical formula (14) are very close to explicit expression (3).So, we believe that PDFs ρ GH t (x) computed from ( 14) with time t � ∆t are also accurate enough. Figure 4 plots the shape of cumulative distribution functions (CDFs) with diferent numbers of M, from which we see the computed CDFs are convergent as M becoming larger.Figure 5 plots the trajectories of nodes X n t and S n t 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 10 − 3 .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.Te results in Table 4 illustrate the efectiveness of the proposed WT method.Te analytical formula of European options under NIG and VG processes can be seen in literatures [2,6,28,29].Tose 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 m of space nodes and time partition number N, some experiments are carried out.Figure 7    ], from which we see the errors are decreasing as N (and so m) increases.Table 5 lists the numerical convergent rates with respect to the values of N. Tese results in Figures 7 and 8 and Table 5 support the theoretical conclusion in Teorem 9.

Conclusions
In this paper, a unifed and robust approach is proposed to construct the willow tree structure for GH Lévy processes.
Tere 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 Δt.Te ffthmoment 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., [30][31][32][33]).We will discuss those models in the future.

Figure 2 :
Figure 2: Graphical depiction of the willow tree lattice with 5 possible asset prices and 4 discrete times.

Figure 3 :
Figure 3: Density functions under (a) GH, (b) HYP, (c) NIG, and (d) VG models.All parameters of these models are listed in Table3.

Figure 4 :
Figure 4: Cumulative distribution functions of (a) GH model and (b) NIG model with ∆t � 0.25 and diferent number of M. All parameters of two models are listed in Table3.

Figure 5 :
Figure 5: Trajectories of X t and S t for the NIG model with m � 10 and N � 30: (a) nodes X n t and (b) nodes S n t .

Figure 7 :
Figure 7: Convergence of willow tree with respect to m.All parameters are listed in Table 3: (a) GH model, (b) HYP model, (c) NIG model, and (d) VG model.

Figure 8 :
Figure 8: Convergence of willow tree with respect to N with fxed m � 50.All parameters are listed in Table 3: (a) GH model, (b) HYP model, (c) NIG model, and (d) VG model.

Figure 8
Figure8plots the errors with diferent N and corresponding m � [N 3/2 ], from which we see the errors are decreasing as N (and so m) increases.Table5lists the numerical convergent rates with respect to the values of N. Tese results in Figures7 and 8and Table5support the theoretical conclusion in Teorem 9.
determining [a, b], the WT algorithm does not consume too much extra computational efort.(3) By setting appropriate m discrete stock prices at each time t n and calculating transform probabilities p (n) Basic conceptions are reviewed in Section 2. Calculation of PDFs and CDFs is illustrated in Section 3. Te 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 fnal section.
[19]rom t n to t n+1 , WT structure is constructed.Unlike Ma et al.[19], it is not needed to estimate the k-order moments of the GH Lévy model when selecting m nodes at each time t n .Te determination of m nodes at each time t n 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)Teconvergence rate O(∆t) 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 ∆t and the underlying partition number m.Te remaining parts of this paper are arranged as follows.
and M is the number of nodes u k  .Proof.It is obvious that %% Input parameters Step 1. Set GH-model parameters (α, β, μ, δ, λ, T). %% Compute probability density function Step 2. Set willow tree parameters (N, m) and computational parameters M. N represents the number of time discretization and m represents the number of X-nodes at each time t n .M is the number in NI formula (14) for PDFs ρ GH t (x).%% Compute probability density function for t � t n and t � Δt.Step 3. Set discrete space nodes u k � A + kΔu with Δu � (B − A)/M for k � 1, 2, • • • , M. Step 4. Compute θ t k and C t k for k � 0, 1, • • • , M − 1 using formulas (15); Step 5.According to (14), compute PDFs ρ GH t (x) on given points x i   i�1: m .%% Generate WT nodes and compute transition probabilities.Step 6.According to (17), compute F n (X n i ) at each time t � t n .By solving (3.1), determine nodes

Table 3 :
Parameters of the GH model for some special cases.

Table 4 :
Option values calculated by WT method, MC method, and analytical formulas.
European call options by WT, MC methods, and analytical formula are labeled by "ECWT," "ECMC," and "ECA," respectively.American put options by WT and MC methods are labeled by "APWT" and "APMC," respectively.ERRs are the corresponding errors between the WT method and Monte Carlo method.
plots the errors for diferent parameters m with fxed N � 20.
In the above formula, K c (x) is the MacDonald function, β(x, y) is the beta function, and Φ(c 1 , c 2 , c 3 ; x, y) is the degenerate Appell hypergeometric function.