Fourier-Cosine Method for Pricing and Hedging Insurance Derivatives

Ludovic Goudenège¹, Andrea Molent², Xiao Wei³, Antonino Zanette⁴

¹Fédération de Mathématiques de CentraleSupélec, Université Paris-Saclay, Paris, France.
²Dipartimento di Management, Università Politecnica delle Marche, Ancona, Italy.
³China Institute for Actuarial Science, Central University of Finance and Economics, Beijing, China.
⁴Dipartimento di Scienze Economiche e Statistiche, Università di Udine, Udine, Italy.

DOI: 10.4236/tel.[phone]                

Abstract

We introduce the Fourier-Cosine method for pricing and hedging insurance derivatives. We implement this method for a particular problem of variable annuities under the Black-Scholes model for the investment account. The numerical results show the reliability of the Fourier-Cosine method for pricing and hedging insurance derivatives.

Keywords

Fourier-Cosine, Pricing, Hedging, Insurance, Derivatives

1. Introduction

Life insurance products with guarantees dependent on stock market indicator are sold in insurance markets of different countries. For instance we can cite united-linked contracts sold in the United Kingdom since the late 1960s, the variable annuities in the United States which form a major part of the sells in annuities market, the segregated fund contracts in Canada since late 1990s and the variable annuities starting to sell from 2011 in China. This kind of insurance products offer the opportunity to earn money from the bull financial market, while providing the protection against bear financial market. The opportunity to earn money and the protection from losing money become an ideal management tool for financing the post-retirement income of the policyholder, this leads to the popularity of variable annuities in many insurance markets. Concerning the population aging problem in the world and the potential inflation risk in economic, the demand of this kind of products will continue to grow.

As far as the financial market is considered, the insurer of the equity-linked insurance (see Hardy [1] ) faces the financial risk in addition to mortality risk and lapse risk. Thus the valuation of this kind of insurance contracts corresponds to find hedging against the risk of the contracts. The financial risk stems from the guaranteed benefits, it can be treated as an embed option in the contract. Properly evaluate the financial risk is essential to the insurer. This kind of contracts caused inefficient financing to their providers during the financial crisis in 2008, forcing some insurers to reconsider the design and pricing of these products. As an insurance risk, the mortality risk is usually hedged by diversification, which means that we can assume that the mortality rate is to be deterministic when there are a sufficiently large number of contracts. In some cases, the lapse risk is also treated as diversifiable, thus it is also assumed as a deterministic probability of the contracts in force at maturity, however for the contract provided guarantees of withdrawal benefit, the lapse risk depends on the guaranteed benefits, then it will be considered in the pricing of the embed option. In order to price the equity-linked insurance it is essential to evaluate the guaranteed benefit as an embed option while considering the lapse risk.

The research about the valuation of equity-linked insurance began with the work of Boyle and Brennan, Schwartz in [2] and [3] . They assumed both mortality and lapse risk to be diversifiable and the stock indicator to follow the Black-Scholes models, thus the embed guaranteed benefits can be priced by applying the famous Black-Scholes formulas in option pricing. More recently a general framework to price different guaranteed benefits has been proposed in [4] using Monte Carlo simulation, where the close formula is not available due to the complexity of some guaranteed benefits, especially when the lapse risk cannot be treated as diversifiable, the guaranteed benefits with the lapse risk is a path dependent option. Other authors use Monte Carlo methods. Bacinello et al. [5] evaluate variable annuities (including GLWBs) are priced using a Monte Carlo approach. Holz et al. [6] and Kling et al. [7] used a Monte Carlo approach to price variables annuities products. As Monte Carlo simulation is time consuming for path dependent option, other methods for option pricing would help to increase the computation efficiency of valuation of the guaranteed benefits. PDE methods for these types of guarantees are proposed in Chen et al. [8] , Forsyth et al. [9] and Belanger et al. [10] .

The Fourier-Cosine method, which was first introduced to option pricing by Fang and Oosterlee in [11] , is based on the risk-neutral option valuation formula (discounted expected payoff approach). Classically the option value of a financial product could be known without the knowledge of the transitional density function. But this density can be recovered from its characteristic function by a truncated Fourier-Cosine expansion. Thus this approach permits to price approximately the option value, by computing its Fourier-Cosine coefficients. Here we propose to introduce this Fourier-Cosine method for pricing and hedging the guaranteed benefits, and to prove its efficiency with a comparison with already known numerical methods.

In order to compare with already known results, we present the implementation of this methodology for the valuation of Guaranteed Lifelong Withdrawal Benefits (GLWB) fully described in [9] for the Black-Scholes model and in [12] for Black-Scholes model with stochastic interest rates.

The paper is organized as follows. In Section 1 we describe the Fourier-Cosine methodology, in Section 2 the treated insurance derivative, and in Section 3 we propose the numerical results. Finally we give some conclusion and remarks about all the methodology and perspectives in the last section.

2. Fourier-Cosine Methodology

Let us give a brief introduction of the COS option pricing method. The concept behind COS pricing method is to recover the conditional density by its characteristic function through Fourier-Cosine expansions. It can be applied for all processes as soon as the characteristic function is available, which includes all affine processes. The method performs impressively, especially when the underlying follows a Lévy processes.

Everything starts from the risk-neutral valuation formula

$v\left(x,t_0\right)={\text{e}}^{-r\Delta t}{\displaystyle {\int }_{-\infty }^{\infty }}v\left(y,T\right)f\left(y|x\right)\text{d}y,$

where $v\left(x\mathrm{,}t\right)$ is the option value, and $x\mathrm{,}y$ can be any increasing functions of the underlying at $t_0$ and T respectively. In order to numerically evaluate this integral, we truncate the integration range, aiming at

$v\left(x\mathrm{,}{t}_0\right)\approx {\text{e}}^{-r\Delta t}{\displaystyle {\int }_a^b}v\left(y\mathrm{,}T\right)f\left(y\mathrm{<mo>\; |\; </mo>}x\right)\text{d}y\mathrm{.}$ (2.1)

The integration range has to been chosen such that the condition

$\left|{\displaystyle {\int }_{ℝ}}f\left(y\mathrm{<mo>\; |\; </mo>}x\right)\text{d}y-{\displaystyle {\int }_a^b}f\left(y\mathrm{<mo>\; |\; </mo>}x\right)\text{d}y\right|<TOL$

is satisfied, where TOL is a pre-specified tolerance value An error analysis of the various approximations can be found in [11] [13] .

The conditional density function of the underlying is then approximated by means of the characteristic function via a truncated Fourier-Cosine expansion of order N, as follows:

$f\left(y\mathrm{<mo>\; |\; </mo>}x\right)\approx \frac{2}{b-a}{\displaystyle \sum }{\text{'}}_{k=0}^{N-1}\Re \left(\varphi \left(\frac{k\text{π}}{b-a}\mathrm{<mo>\; ;\; </mo>}x\right)exp\left(-i\frac{ak\text{π}}{b-a}\right)\right)cos\left(k\text{π}\frac{y-a}{b-a}\right)\mathrm{,}$ (2.2)

where $\Re$ means taking the real part of the expression in brackets, and $\varphi \left(\omega \mathrm{<mo>\; ;\; </mo>}x\right)$ is the characteristic function of $f\left(y\mathrm{<mo>\; |\; </mo>}x\right)$ defined as:

$\varphi \left(\omega \mathrm{<mo>\; ;\; </mo>}x\right)=\mathbb{E}\left[{\text{e}}^{i\omega y}\mathrm{<mo>\; |\; </mo>}x\right]\mathrm{.}$ (2.3)

The prime at the sum symbol in (2.2) indicates that the first term in the expansion is multiplied by one-half. Replacing $f\left(y\mathrm{<mo>\; |\; </mo>}x\right)$ in (2.1) by its approximation (2.2) and interchanging integration and summation, we obtain the COS algorithm to approximate the value of a European option:

$v\left(x,t_0\right)={\text{e}}^{-r\Delta t}{\displaystyle \sum }{\text{'}}_{k=0}^{N-1}\Re \left(\varphi \left(\frac{k\text{π}}{b-a};x\right)\exp \left(-i\frac{ak\text{π}}{b-a}\right)\right)V_k,$ (2.4)

where

$V_k=\frac{2}{b-a}{\displaystyle {\int }_a^b}v\left(y\mathrm{,}T\right)cos\left(k\text{π}\frac{y-a}{b-a}\right)\text{d}y$ (2.5)

is the Fourier-Cosine coefficient of $v\left(y\mathrm{,}T\right)$ , which is available in closed form for several payoff functions of European options.

The COS algorithm exhibits an exponential convergence rate for all processes whose conditional density $f\left(y\mathrm{<mo>\; |\; </mo>}x\right)$ is in $C^{\infty }\left(\left(a\mathrm{,}b\right)\subset ℝ\right)$ . The size of the integration interval $\left[a\mathrm{,}b\right]$ can be determined with help of the cumulants [11] .

3. Insurance Derivatives

In this section, divided into two subsections, we present the Guaranteed Lifelong Withdrawal Benefits contract description. In the first subsection we describe the dynamics of the proposed insurance product between monitoring dates, and in the second subsection we describe the events occurring at monitoring dates.

3.1. Guaranteed Lifelong Withdrawal Benefits Model

In this section we will consider the index $\hat{S}$ which follows a Geometric Brownian motion with initial value ${\hat{S}}_0$ , constant drift r and volatility $\sigma$ (Black-Scholes model). It satisfies the following stochastic differential equation:

$\text{d}{\hat{S}}_t=r{\hat{S}}_t\text{d}t+\sigma {\hat{S}}_t\text{d}{W}_t,$ (3.1)

where $W_t$ is a Brownian motion.

Let the mortality function $M\left(t\right)$ be defined as the fraction of the original owners of the Guaranteed Lifelong Withdrawal Benefits (GLWB) contract who die in the interval $\left[t\mathrm{,}t+\text{d}t\right]$ ( $M\left(t\right)$ is assumed to be constant in this fixed interval of size $\text{d}t$ ). This number can be computed from a mortality table, given the probability $\mathbb{Q}\left(x\right)$ that an owner which is x years old will die the next year. The fraction of the original owners still alive at time t, denoted by $R\left(t\right)$ is

$R\left(t\right)=1-{\displaystyle {\int }_0^t}M\left(s\right)\text{d}s\mathrm{.}$

Time t is measured in years from the contract inception date, so the maturity corresponds to the moment when everyone has died (and of course the end of the mortality table).

Let S be the amount in the investment account (i.e. mutual fund) of any holder of the GLWB contract still alive at time t. We suppose that this investment account tracks an index $\hat{S}$ previously described. Let A be the guarantee account balance. We suppose that percentage fees based on the value of the investment account S are charged to the policy holder at the annual rate ${\alpha }_{tot}$ and withdrawn continuously from that account. These fees include mutual fund management fees ${\alpha }_m$ and a fee charged to fund the guarantee (also known as the rider) ${\alpha }_g$ , so that ${\alpha }_{tot}={\alpha }_g+{\alpha }_m$ . Let $V\left(S\mathrm{,}A\mathrm{,}t\right)$ be the value of the entire contract (express in backward time) as sum of the no-arbitrage value of the guarantee only portion of the contract (the GLWB rider) and the amount in the investment accounts of those remaining alive. We have the following dynamics for $V\left(S\mathrm{,}A\mathrm{,}t\right)$ :

$V_t=\frac{{\sigma }^2{S}^2}{2}{V}_{SS}+\left(r-{\alpha }_{tot}\right)S{V}_S-rV+{\alpha }_{m}R\left(t\right)S+M\left(t\right)S\mathrm{,}$ (3.2)

with initial value $V\left(S,A,0\right)=0$ (which is in fact the terminal condition in real life time).

Between two annual dates $t_i$ and $t_{i+1}$ , the contract follows this partial differential equation. At date $t_{i+1}$ there are jumps depending on the specification of the contract.

3.2. Jumps Events

At each monitoring date $t_i$ there are jumps events to take into account the product features, e.g. withdrawal, bonus, ratchet, penalties. We assume that the order of event occurring at an event time $t_i$ is first the ratchet event then the withdrawal events containing possibly bonus or penalties.

Ratchet Event. If the contract specifies a ratchet (step-up) feature, then the value of the guarantee account A is increased if the investment account has increased. The guarantee account A can never decrease, unless the contract is partially or fully surrendered. At a ratchet event time $t_i$ we have the following relation

$V\left(S\mathrm{,}A\mathrm{,}{t}_i^{+}\right)=V\left(S\mathrm{,}max\left(S\mathrm{,}A\right)\mathrm{,}{t}_i^{-}\right)\mathrm{.}$

General Withdrawal Event. The contract will typically specify a withdrawal rate $G_r$ . Given a time interval of $t_i-t_{i-1}$ between withdrawals, the contract withdrawal amount at $t=t_i$ is $G_r\left(t_i-t_{i-1}\right)A$ . At this point we do not make any particular assumptions about the withdrawal strategy of the policy holder. In general terms, the policy holder’s actions at $t_i$ can be represented by a policy parameter ${\gamma }_i$ , where $0\le {\gamma }_i\le 2$ . Withdrawals of amounts less than or equal to the contract withdrawal amount $G_r\left(t_i-t_{i-1}\right)A$ are represented by ${\gamma }_i\in \left[\mathrm{0,1}\right]$ . Withdrawals in excess of the contract amount are indicated by ${\gamma }_i\in \left(\mathrm{1,2}\right]$ , with ${\gamma }_i=2$ corresponding to full surrender. Withdrawal events can be written in the general form

$V\left(S\mathrm{,}A\mathrm{,}{t}_i^{+}\right)=V\left(S^{{\gamma }_i}\mathrm{,}{A}^{{\gamma }_i}\mathrm{,}{t}_i^{-}\right)+\text{CashFlow}\left(S\mathrm{,}A\mathrm{,}{t}_i\mathrm{,}{\gamma }_i\right)\mathrm{,}$

where $S^{{\gamma }_i}$ and $A^{{\gamma }_i}$ are particular values depending on the withdrawal event, and where “Cash Flow” is the cash flow from the event depending on the withdrawal event represented by the value ${\gamma }_i$ .

Bonus Event ( ${\gamma }_i=0$ ). If the contract holder chooses not to withdraw at $t=t_i$ , this is indicated by ${\gamma }_i=0$ . Let the bonus fraction be denoted by $B\left(t_i\right)$ (if no bonus is possible at $t=t_i$ , then $B\left(t_i\right)=0$ ). Thus we define

$S^{{\gamma }_i}\mathrm{<mo>\; :\; </mo>}=S\mathrm{,}{A}^{{\gamma }_i}\mathrm{<mo>\; :\; </mo>}=A\left(1+B\left(t_i\right)\right)\text{and}\text{CashFlow}=0.$

Withdrawal Not Exceeding Contract Amount ( ${\gamma }_i\in \left(\mathrm{0,1}\right]$ ). Note that withdrawals at the contract rate (or less) are allowed even if the amount in the investment account is $S=0$ . The withdrawal amount is ${\gamma }_i{G}_r\left(t_i-t_{i-1}\right)A$ and we obtain the following definitions

$\begin{array}{l}{S}^{{\gamma }_i}\mathrm{<mo>\; :\; </mo>}=max\left(S-{\gamma }_i{G}_r\left(t_i-t_{i-1}\right)A\mathrm{,0}\right)\mathrm{,}\\ A^{{\gamma }_i}\mathrm{<mo>\; :\; </mo>}=A\text{and}\text{CashFlow}=R\left(t_i\right){\gamma }_i{G}_r\left(t_i-t_{i-1}\right)A\mathrm{.}\end{array}$

Partial or Full Surrender ( ${\gamma }_i\in \left(\mathrm{1,2}\right]$ ). Finally consider the case of a withdrawal of an amount greater than the contract amount $G_r\left(t_i-t_{i-1}\right)A$ , i.e. the withdrawal amount is given by

$G_r\left(t_i-t_{i-1}\right)A+\left({\gamma }_i-1\right)max\left(S-G_r\left(t_i-t_{i-1}\right)A\mathrm{,0}\right)\left(1-\kappa \left(t_i\right)\right)$

where $\kappa \left(t_i\right)\in \left[\mathrm{0,1}\right]$ is a penalty for withdrawal above the contract amount. In this case we have to define the following quantities

$S^{{\gamma }_i}:=\left(2-{\gamma }_i\right)\max \left(S-{\gamma }_i{G}_r\left(t_i-t_{i-1}\right)A,0\right),A^{{\gamma }_i}:=A\left(2-{\gamma }_i\right)$

and

$\text{CashFlow}=R\left(t_i\right)G_r\left(t_i-t_{i-1}\right)A+\left({\gamma }_i-1\right)max\left(S-G_r\left(t_i-t_{i-1}\right)A\mathrm{,0}\right)\left(1-\kappa \left(t_i\right)\right)\mathrm{.}$

Note that it is assumed that the guarantee account value A is reduced proportionately for any withdrawal above the contract rate.

Now we consider the withdrawal strategy. The risk neutral price is the cost of hedging. If we consider that the insurer should charge a price which ensures that no losses can occur (assuming that the claim is hedged), then the withdrawal strategy is assumed to verify

${\gamma }_i=\underset{\gamma \in \left[0,2\right]}{{\displaystyle \arg \max }}\left(V\left(S^{\gamma },A^{\gamma },t_i^{-}\right)+\text{CashFlow}\left(S,A,t_i,\gamma \right)\right).$

Assuming such a strategy by policy holders, and hedging against it, is obviously very conservative from the standpoint of the insurer, since it seeks to provide complete protection against policy holder withdrawal behavior (given assumptions about parameter values such as volatilities). In other words, if investors follow this strategy, and if the insurer hedges continuously, the balance in the insurer’s overall hedged portfolio will be zero. On the other hand, if the investor deviates from this strategy, then the insurer’s portfolio will have a positive balance.

4. Numerical Results

At each date $t_i$ , we need to compute the withdrawal strategy, so we look for the maxima by a linear search on values $\gamma \in \left[\mathrm{0,2}\right]$ . From the previous experiments on this model, the authors can say that only values 0, 1 and 2 are actually selected by this linear search. That is to say: when withdrawal is chosen, only full withdrawal is selected, and in case of surrender event, only a full surrender event is selected.

With this optimal dynamical approach, or in a static approach where all ${\gamma }_i$ values are fixed at the beginning of the contract, it is possible to compute (backward in real life time, but forward here) the value $V\left(S\mathrm{,}A\mathrm{,}t\right)$ for all variables S, A and $t\in \left[\mathrm{0,}T\right]$ .

The pricing problem now reduces to find the rider fee ${\alpha }_g^{\mathrm{<mo>\; *\; </mo>}}$ (driving the dynamic between time events) such that

$V\left[{\alpha }_g^{\mathrm{<mo>\; *\; </mo>}}\right]\left(S_0\mathrm{,}{A}_0\mathrm{,}T\right)=S_0\mathrm{,}$ (4.1)

when $V\left[{\alpha }_g^{\mathrm{<mo>\; *\; </mo>}}\right]$ is the solution of the partial differential Equation (3.2) when ${\alpha }_g={\alpha }_g^{*}$ . Viewing $V\left[{\alpha }_g\right]$ as being parametrized by the rider fee ${\alpha }_g$ , we solve the Equation (4.1) using a classical secant method. Typically, only 5 or 7 iterations are necessary to obtain convergence of the algorithm under a fixed tolerance of 10⁻⁸.

4.1. Numerical Parameters

In all our numerical experiments we have used the following parameters which are compatible with actual literature on GLWB. Table 1 refers to the Black-Scholes parameters, Table 2 refers to the product features and Table 3 refers to the policy holder parameters describing its behavior.

We consider that there are no initial fees and that all the fees and the death benefits are paid continuously. There is no ratchet feature. The bonus is of 5% if there is no withdrawal. If there is a partial surrender then there are penalties given in Table 4.

Moreover we used the DAV 2004R mortality Table, 65 year old German male (see [9] for the table).

We have used the following parameters for the Fourier-Cosine method. The range of integration is given by $a=\left(r-{\alpha }_g\right)-L\sigma +log\left(S_0\right)$ and

$b=\left(r-{\alpha }_g\right)+L\sigma +log\left(S_0\right)$ with $L=17.5$ . We have remarked that a number

Table 1. The Black-Scholes parameters.

Table 2. The product features.

Table 3. The policy holder parameters.

Table 4. The surrender penalties.

$N=20$ of Fourier-Cosine coefficients provides a fast method, but the choice of $N=90$ seems better and more stable.

For the partial differential equation method (PDE), we use an implicit finite difference scheme with a time step ${\delta }_t=0.01$ and a non-uniform mesh of 750 points in the space variable.

4.2. Pricing and Hedging the Insurance Product

We compute the fee ${\alpha }_g$ such that (4.1) is verified with a tolerance of 10⁻⁵ i.e.

$\left|V\left[{\alpha }_g^{*}\right]\left(S_0,A_0,T\right)-S_0\right|<{10}^{-5}.$

We perform all these computations with parameters detailed in Tables1-4. We only change the value of the volatility parameter, respectively $\sigma =0.10$ and $\sigma =0.15$ .The numerical results are presented in Table 5.

With this method it is possible to numerically compute the classical greeks. But in the insurance market, and for this kind of product, an interesting value could be given by the sensibility with respect to the mortality table. The authors have implemented a ±10% shock on the probability $\mathbb{Q}\left(x\right)$ given by the mortality table. Table 6 gives the sensibility with respect to this shock when $\sigma =0.15$ .

The numerical results provide that Fourier-Cosine method may be a reliable method in order to price variables annuities. The parameters defining the Fourier-Cosine method have been chosen in order to obtain fast and accurate results.

5. Conclusions

In this article we have introduced the Fourier-Cosine method for pricing and hedging insurance derivatives. We have investigated the efficiency of such a method on a particular insurance product which is the Guaranteed Long life Withdrawal Benefits model of variable annuities. We have assumed that the index tracked by the investment account follows a Geometric Brownian Motion, but this method can be extended to various models as soon as we have an explicit formula for the density function for the underlying process.

First numerical results say us that Fourier-Cosine method seems to be a promising technique in order to price and hedge variable annuities.

Future research can be extended to take into account a model for the index

Table 5. Values ${\alpha }_g$ of the fee.

Table 6. Sensibility with respect to ±10% shocks on mortality probabilities.

with stochastic volatility and/or stochastic interest rate. For instance in the Heston model we have an explicit formula for the density function, then the extension to this case could be straightforward.

Conflicts of Interest

The authors declare no conflicts of interest.

References

[1]	Hardy, M. (2003) Investment Guarantees: Modeling and Risk Management for Equity Linked Life Insurance. Wiley, New York.
[2]	Boyle, P.P. and Schwartz, E. (1977) Equilibrium Prices of Guarantees under Equity-Linked Contracts. Journal of Risk and Insurance, 44, 639-680. https://doi.org/10.2307/251725
[3]	Brennan, M.J. and Schwartz, E. (1976) The Pricing of Equity-Linked Life Insurance Policies with an Asset Value Guarantee. Journal of Financial Economics, 3, 195-213. https://doi.org/10.1016/0304-405X(76)90003-9
[4]	Bauer, D., Kling, A. and Russ, J. (2008) A Universal Pricing Framework for Guaranteed Minimum Benefits in Variable Annuities. Astin Bulletin, 38, 621-651. https://doi.org/10.1017/S0515036100015312
[5]	Bacinello, A.R., Millossovich, P., Olivieri, A. and Pitacco, E. (2011) Variable Annuities: A Unifying 953 Valuation Approach. Insurance: Mathematics and Economics, 49, 285-297. https://doi.org/10.1016/j.insmatheco.2011.05.003
[6]	Holz, D., Kling, A. and Ruβ, J. (2007) GMWB for Life: An Analysis of Lifelong Withdrawal Guarantees. Zeitschrift für die gesamte Versicherungswissenschaft, 101, 305-325. https://doi.org/10.1007/s12297-012-0193-3
[7]	Kling, A., Ruez, F. and Ruβ, J. (2014) The Impact of Stochastic Volatility on Pricing, Hedging, and Hedge Efficiency of Variable Annuity Guarantees. European Actuarial Journal, 4, 281-314. https://doi.org/10.1007/s13385-014-0093-0
[8]	Chen, Z., Vetzal, K. and Forsyth P. (2008) The Effect of Modelling Parameters on the Value of GMWB Guarantees. Insurance: Mathematics and Economics, 43, 165-173. https://doi.org/10.1016/j.insmatheco.2008.04.003
[9]	Forsyth, P. and Vetzal, K. (2014) An Optimal Stochastic Control Framework for Determining the Cost of Hedging of Variable Annuities. Journal of Economic Dynamics and Control, 44, 29-53. https://doi.org/10.1016/j.jedc.2014.04.005
[10]	Belanger, A., Forsyth, P. and Labahn, G. (2009) Valuing the Guaranteed Minimum Death Benet 957 Clause with Partial Withdrawals. Applied Mathematical Finance, 16, 451-496. https://doi.org/10.1080/13504860903075464
[11]	Fang, F. and Oosterlee, C.W. (2008) A Novel Pricing Method for European Options Based on Fourier-Cosine Series Expansions. SIAM Journal on Scientific Computing, 32, 826-848.
[12]	Fang, F. and Oosterlee, C.W. (2010) A Fourier-Based Valuation Method for Bermudan and Barrier Options under Heston Model. SIAM Journal on Financial Mathematics, 2, 439-463. https://doi.org/10.1137/100794158
[13]	Goudenège, L., Molent, A. and Zanette, A. (2016) Pricing and Hedging GLWB in the Heston and in the Black-Scholes with Stochastic Interest Rate Models. Insurance: Mathematics and Economics, 70, 38-57. https://doi.org/10.1016/j.insmatheco.2016.05.018