Notes on the SWIFT method based on Shannon Wavelets for Option Pricing
aa r X i v : . [ q -f i n . C P ] M a y May 28, 2020 0:34 swift˙notes
Notes on the SWIFT method based on ShannonWavelets for Option Pricing
Fabien Le Floc’h (v0.1 released January 2018)
Abstract
This note shows that the cosine expansion based on the Vieta formula is equivalent to a discretization ofthe Parseval identity. We then evaluate the use of simple direct algorithms to compute the Shannon coefficients for thepayoff. Finally, we explore the efficiency of a Filon quadrature instead of the Vieta formula for the coefficients relatedto the probability density function.
Key Words:
SWIFT method, Wavelets, Heston, stochastic volatility, characteristic function, quantitative finance
1. Introduction
Ortiz-Gracia and Oosterlee (2016) describe a novel approach to the pricing of European op-tions under models with a known characteristic function, based on Shannon Wavelets, referredto as the SWIFT method hereafter. This note shows that the cosine expansion based on Vi-eta’s formula is equivalent to a discretization of Parseval’s identity. We then evaluate theuse of simple direct algorithms to compute the Shannon coefficients for the payoff. Finally,we explore the efficiency of a Filon quadrature instead of Vieta’s formula for the coefficientsrelated to the probability density function.The equivalence with Parseval’s identity is also stated in (Maree et al. , 2017).
2. Equivalence with Parseval’s identity
With the SWIFT method, the price at time t of a Vanilla Put option of maturity T andlog-moneyness x = ln FK , with K the strike and forward F is v ( x, t ) = B ( t, T ) k X k = k c m,k V m,k (1)where c m,k = h f | φ m,k i = 2 m Z R f ( x ) φ (2 m x − k ) dx, , (2) V m,k = Z I m v ( y, T ) φ m,k ( y ) dy , (3) φ m,k ( x ) = 2 m φ (2 m x − k ) , (4) φ ( x ) = sin πxπx (5) ISSN: Print/ISSN Online/YY/issnum1stp–ppcnt © ay 28, 2020 0:34 swift˙notes Fabien Le Floc’h and k , k , m ∈ Z , m > = 1 suitably chosen, f the probability density function and v ( y, T )is the payoff at maturity with y = ln F ( T,T ) K , that is v ( y, T ) = K | − e y | + for a vanilla Putoption.In (Ortiz-Gracia and Oosterlee, 2016), the coefficients c m,k and V m,k are computed using anapproximation based on Vieta formula for the cardinal sinus: φ ( x ) ≈ J − J − X j =1 cos (cid:18) j − J πx (cid:19) (6)where J is chosen sufficiently large.As mentioned in paragraph 3.1.2 of their paper, c m,k can also be computed by Parseval’sidentity: c m,k = h f | φ m,k i = 12 π D ˆ f | ˆ φ m,k E (7)where ˆ f , ˆ φ m,k are the Fourier transforms of f and φ m,k . In particular ˆ f ( z ) = ψ ( − z ) where ψ is the characteristic function andˆ φ m,k ( w ) = e − i k m w m rect (cid:16) w m +1 π (cid:17) (8)where rect is the rectangular function, that is rect( x ) = 1 for | x | < , rect( x ) = for | x | = ,rect( x ) = 0 for | x | > .Via a the change of variable t = w m +1 π , we obtain c m,k = 2 m Z − h ˆ f (2 m +1 πt ) e i πkt i dt (9)= 2 m +1 ℜ "Z ˆ f (2 m +1 πt ) e i πkt dt (10)as ℜ ( ψ ( x )) = ψ ( x )+ ψ ( x )2 = ψ ( x )+ ψ ( − x )2 .Let us now discretize in 2 J − equidistant steps equation (10) of size J at the mid-points t j = j − J for j = 1 , , ..., J − , we obtain c ⋆m,k = 2 m J − J − X j =1 ℜ h ˆ f (2 m π t j ) e iπk t j i (11)= 2 m J − ℜ J − X j =1 ˆ f (cid:18) m π (2 j − J (cid:19) e iπk j − J (12)= 2 m J − ℜ e iπ k J J − − X j =0 ˆ f (cid:18) m π (2 j − J (cid:19) e iπk j J (13)This is exactly equation (24) of (Ortiz-Gracia and Oosterlee, 2016, p. B127) which correspondsto their expansion based on Vieta’s formula. Their expansion is thus equivalent to the mid-point quadrature applied to Parseval’s identity.A particularly important property of equation 13 is that it can be computed by fast Fouriertransform (FFT). The typical FFT algorithm computes the tranform (or inverse transform)from index 0 to n . Here, we start with a negative index k . The coefficients can be obtained ay 28, 2020 0:34 swift˙notes Draft with the relation c = 2 m J − ℜ (cid:2) e T F − { f } (cid:3) (14)where F − is the unscaled inverse discrete Fourier transform of size 2 J , the vector f haselements f j = ˆ f (cid:16) m π (2 j +1)2 J (cid:17) e iπ k j J and the vector e has elements e l = e iπ l + k J . We alsoassumed that ˆ f (cid:16) m π (2 j +1)2 J (cid:17) = 0 for j ≥ J − . This leads to a very efficient way to computethe coefficients c m,k , for all k , together. In practice, this means that the bounds k ≤ k < k are chosen so that k − k < J .In particular, if we center the interval around zero, that is for k = − J − , we can savea bit of computation by directly using f j = ˆ f (cid:16) m π (2 j +1)2 J (cid:17) and swapping ( g , ..., g J − ) with( g J − , ...g J − ) where g = F − { f } .
3. Alternative quadratures
Instead of the mid-point method, we could have considered the trapezoidal method, this wouldresult in c ⋆m,k = 2 m J − ℜ J − − X j =0 w j ˆ f (cid:18) m π (2 j )2 J (cid:19) e iπk j J (15)where w j = 1 for j ≥ w = .The fast inverse discrete Fourier transform of length 2 J can be directly used to compute c m,k by using f j = ˆ f (cid:16) m π (2 j )2 J (cid:17) for 1 ≤ j < J − , f j = 0 for 2 J − ≤ j , and f = ˆ f (0).We will see in the numerical examples that it can be much more accurate than the mid-pointmethod.In the same framework, we could also explore other quadratures, such as the Simpson’squadrature. The problem is that those tend to behave worse than the midpoint or trapezoidalrules on oscillatory functions. In fact, the trapezoidal rule can achieve exponential convergenceon oscillatory functions (Johnson, 2011; Trefethen and Weideman, 2014). In the case of theprobability density transform function ˆ f , this can be also be seen from the Euler-Maclaurinformula where as all the derivatives ˆ f (2 l +1) (2 m π ) will be small if the characteristic functiondecreases exponentially. Instead of quadrature with a fixed number of steps, we can use an adaptive Filon quadratureto compute the coefficients c m,k by equation (10). This is particularly interesting since thecost of computing the characteristic function is relatively high.Is it more important to reduce its number of evaluations than to use Fast Fourier Transformtricks to compute c m,k ? We will explore this in the numerical examples (section 8).An alternative adaptive quadrature, close in spirit, is to use an adaptive cubic-Hermitequadrature to integrate ˆ f , and use the integration nodes to compute the piecewise cubicHermite interpolant of ˆ f . Then we can use the trapezoidal-FFT approach on a dense dis-cretization. This saves explicit computations of the characteristic function while still allowingthe use of the FFT algorithm. ay 28, 2020 0:34 swift˙notes Fabien Le Floc’h
4. Sine and Exponential integrals for the payoff
For a Vanilla Put option, the payoff at maturity is V ( y, T ) = K (1 − e y ) + . According toequation (3), the payoffs coefficients are then V m,k = K m Z a (1 − e y ) sin ( π (2 m y − k )) π (2 m y − k ) dy (16)= K m Z a sin ( π (2 m y − k )) π (2 m y − k ) dy − K m Z a e y sin ( π (2 m y − k )) π (2 m y − k ) dy (17)= K m π Z − πkπ (2 m a − k ) sin tt dt − Ke k m m π Z − πkπ (2 m a − k ) e tπ m sin tt dt (18)The first integral corresponds the sine integral Si ( x ) = R x tt dt . Many efficient algorithmsexist to compute it (MacLeod, 1996; Jin and Jjie, 1996). Most mathematical software (for ex-ample Octave, Matlab) or libraries (for example netlib) include the function. It can effectivelybe considered as a closed form function.The second integral can be reduced to evaluations of the complementary exponential integral Ein ( z ) = R z − e − t t dt in the complex plane. In deed, it can be verified that we have the identity Z e − at sin( bt ) t dt = ℑ Ein ( a + ib ) (19)The complementary exponential integral is related to the exponential integral Ei ( z ) = − R ∞ z e − t t dt by the relation Ein ( z ) = γ + ln | z | + i ℑ ( − z ) | arg( − z ) | − Ei ( − z ). Again manyefficient algorithms exist to compute the complementary exponential integral (Amos, 1990;Jin and Jjie, 1996; Pegoraro and Slusallek, 2011).In terms of those special functions, the coefficients are: V m,k = K m π e k m ℑ (cid:20) Ein (cid:18) − t a π m + it a (cid:19) − Ein (cid:18) − t π m + it (cid:19)(cid:21) − K m π [ Si ( t a ) − Si ( t )] (20)with t a = π (2 m a − k ) and t = − πk .The expansion based on Vieta’s formula might require thousands of terms to reach anacceptable accuracy (Table 1). With the same number of terms, a Simpson 3/8 quadratureis more accurate and faster to compute. Our simple implementation of the algorithm fromPegoraro and Slusallek (2011) is much faster and achieves machine epsilon accuracy while thealgorithm from the CERN libary Mathlib (K¨olbig, 1990) is even faster for a close to machineepsilon accuracy as it relies on simple rational and pad´e expansions in the zone of interest.In practice, the implementation of the SWIFT method will still benefit from a cache table of V m,k for example for m ∈ { , ..., } and k ∈ {− , ..., } .Table 1.: V m,k for m = 6, k = − a = −
1. Vieta’s formula or Simpson’s quadrature use 2 J − terms.Method Value Time(ns)Vieta J = 5 -0.0555195115435162 600Simpson J = 5 -0.0020905045216672 520Vieta J = 10 0.0020428901436639 17300Simpson J = 10 0.0020420973936057 15300CERN 0.0020420954069492 420Pegoraro 0.0020420954069488 2500 ay 28, 2020 0:34 swift˙notes Draft
5. Alternative payoff coefficients
The interval [ a, b ] is centered along the spot F (0 , T ), we can express the payoff in terms ofthe spot F (0 , T ) instead of the strike K . This leads to V m,k = 2 m Z ba F (cid:12)(cid:12)(cid:12)(cid:12) KF − e y (cid:12)(cid:12)(cid:12)(cid:12) + sin ( π (2 m y − k )) π (2 m y − k ) dy (21)= F m Z za ( e z − e y ) sin ( π (2 m y − k )) π (2 m y − k ) dy (22)= Ke − z m Z za ( e z − e y ) sin ( π (2 m y − k )) π (2 m y − k ) dy (23)where z = ln KF and y = ln S T F .In terms of those special functions, the coefficients are: V m,k ( z ) = Ke k m − z m π ℑ (cid:20) Ein (cid:18) − t a π m + it a (cid:19) − Ein (cid:18) − t z π m + it z (cid:19)(cid:21) − K m π [ Si ( t a ) − Si ( t z )] (24)with t a = π (2 m a − k ) and t z = π (2 m z − k ).The price of the option of strike K corresponds then to v (0 , t ). The coefficients c m,k become independent of x . This is nearly equivalent to the Levy based equation (33) in(Ortiz-Gracia and Oosterlee, 2016) that defines the coefficients V αm,k ( x ). The difference lies inthe interval considered. In their paper, V αm,k ( x ) = R ba K | − e u | + φ m,k ( u + z ) du with u = ln S T K .This can be rewritten as V αm,k ( x ) = K Z za + z (cid:0) − e t − z (cid:1) φ m,k ( t ) dt (25)with the change of variable t = u + z . The interval [ a, b ] is thus shifted from z upwards.Our choice of interval is more accurate as it corresponds directly to the Levy characteristicfunction, while their interval is based on the shifted Levy characteristic function. Also theirLevy formulation (as well as ours) leads to options prices different from the classic formulation:for the two to be equivalent, the integers k and k should be adjusted to k = 2 m ( a + z ) and k = 2 m ( b + z ). But then some of the density coefficients need to be recomputed at each strikeas the window [ k , k ] moves forward as z increases and the Levy approach loses in efficiency.Another advantage of having c m,k independent of the strike is that the integers k and k canalso be determined in a strike independent manner from the value of the density coefficients c m,k and the area under the curve defined by the probability density (which should sum to oneminus a user-defined tolerance) as explained by Ortiz-Gracia and Oosterlee (2016), instead ofrelying on the relatively rough guess given by the cumulants (the fixed interval [ a, b ]). Withthe cumulants approach, it is not always obvious how large the truncation level L should bechosen to achieve a desired accuracy.
6. Alternative FFT-compatible payoff coefficients
In a similar fashion to Maree et al. (2017), we start from the definitionsin( πx ) πx = Z cos( πxw ) dw . (26)We can then choose an appropriate discretization that has good convergence, andis allows computation of the payoff coefficients V m,k by the FFT. The choice from ay 28, 2020 0:34 swift˙notes Fabien Le Floc’h
Ortiz-Gracia and Oosterlee (2016) is equivalent to the mid-point quadrature. On this prob-lem, the Trapezoidal rule would not lead to an increase in accuracy . A particularly simple aneffective choice is the second Euler-Maclaurin summation formula, that is the Euler-Maclaurinextension to the mid-point rule.sin( πx ) πx ≈ N N − X n =0 cos (cid:18) πx n + 12 N (cid:19) + πx N sin( πx ) . (27)Using the above in equation (23) leads to V m,k = Ke − z m N N − X n =1 Z za ( e z − e y ) cos (cid:18) π (2 m y − k ) 2 n + 12 N (cid:19) dy + π N Ke − z m Z za (2 m y − k ) ( e z − e y ) sin ( π (2 m y − k )) dy (28)Let C n ( a, z ) = R za ( e z − e y ) cos (cid:0) π m y nN (cid:1) dy and S n ( a, z ) = R za ( e z − e y ) sin (cid:0) π m y nN (cid:1) dy .Using the trigonometric cos and sin identities, we obtain V m,k = Ke − z m N N − X n =0 C n + ( a, z ) cos (cid:18) πk n + 12 N (cid:19) + S n + ( a, z ) sin (cid:18) πk n + 12 N (cid:19) − ( − k πk N Ke − z m S N ( a, z ) − ( − k N Ke − z m D ( a, z ) (29)with D ( a, z ) = e z (cid:0)(cid:0) p m + p m (cid:1) z − p m − (cid:1) sin ( p m z ) + (cid:0)(cid:0) p m + p m (cid:1) z + 2 p m (cid:1) cos ( p m z ) p m ( p m + 2 p m + 1)+ e z (cid:0) p m + 2 p m + 1 (cid:1) sin ( a p m ) + (cid:0) − a p m − a p m − a p m (cid:1) cos ( a p m ) p m ( p m + 2 p m + 1)+ e a (cid:0) ( − a − p m + (1 − a ) p m (cid:1) sin ( a p m ) + (cid:0) a p m + ( a − p m (cid:1) cos ( a p m ) p m ( p m + 2 p m + 1) , (30) C n ( a, z ) = e z sin ( q n,m z ) − q n,m cos ( q n,m z ) q n,m (1 + q n,m ) − e z sin ( q n,m a ) q n,m + e a cos ( q n,m a ) + q n,m sin ( q n,m a )1 + q n,m , (31) S n ( a, z ) = − e z cos ( q n,m z ) + q n,m sin ( q n,m z ) q n,m (1 + q n,m ) + e z cos ( q n,m a ) q n,m + e a sin ( q n,m a ) − q n,m cos ( q n,m a )1 + q n,m , (32) q n,m = nN p m , (33) p m = π m . (34)In particular, S N and D are independent of k . Computing V m,k with the Euler-Maclauringcorrection for all k requires only k more multiplications than the mid-point quadrature. The It can be shown that the mid-point is actually more accurate by a factor of two.ay 28, 2020 0:34 swift˙notes
Draft sum over n corresponds to the mid-point quadrature, and can be computed with two fastFourier transforms of size N (see appendix A).
7. Choice of m and k , k The SWIFT method accuracy is fully determined by the choice of the scale m and the trun-cation k , k . There is some interplay between those since the scale m also determines thetruncation of the characteristic function: the characteristic function will not be evaluatedbeyond 2 m π .If we want to use the radix-2 FFT algorithm to compute the payoff coefficients V m,k , thereis little reason not to use k − k = 2 J , centered on zero, where a reasonably good guess for J can be obtained from the model characteristic function cumulants. In the evaluation of asingle option strike, the cost of computing the payoff coefficients will dominate the cost ofevaluating the price based on the sum of the V m,k multiplied by the (precomputed) densitycoefficients c m,k . Furthermore, the number of coefficients must be a power of two and mustinclude [ k , k ).The scale m is more challenging to guess. It can be guessed from the rule used to truncate theintegral of the more standard Fourier based approach from Andersen and Piterbarg (2010),refined in Le Floc’h (2013). It then directly depends on the asymptotic behaviour of thecharacteristic function. Maree et al. (2017) propose a simple iterative method to determine n automatically (with very few iterations on m ).
8. Numerical examples V m,k and the FFT Vieta’s formula is not very efficient to compute a single coefficient V m,k but as we computeclose to 2 J coefficients the FFT improves its performance significantly. For 2 V m,k for k = 0 , ..., J − − m = 6, a = − J FFT CERN5 1.7 10.710 56 360While the raw difference in performance is impressive. It is more interesting to look at theactual performance difference when pricing vanilla Put options under the Heston stochasticvolatility model. We consider two different Heston parameter sets for two distinct optionmaturities. This leads to two vastly different truncation ranges [ a, b ], computed according tothe Heston cumulants. As a result 2 J = 2 = k − k for the first set and 2 J = 2 = k − k for the second set. Ignoring the initialization time where the c m,k are computed, which needsto be done only once per option expiry, the direct CERN algorithm is between five to eighttimes slower. ay 28, 2020 0:34 swift˙notes Fabien Le Floc’h
Table 3.: Heston parameter sets.Name v κ θ σ ρ F T Set 1 0.1 1.0 0.1 1.0 -0.9 1.0 2 daysSet 2 0.0225 0.1 0.01 ,2.0 0.5 1000000 1 yearTable 4.: Time in milliseconds taken to compute the Put option price under two differentHeston parameter sets.Heston Method Price Error Time (ms)Set 1 (J=12) FFT 117.9149 -1.4704 0.250CERN 117.9144 -1.4708 1.370Set 2 (J=8) FFT 0.006361 -3.49e-15 0.016CERN 0.006361 6.42e-13 0.101
We consider options of maturity 2 days (short) in order to make the issue more visible andwe consider the Heston parameters s κ = 1 . , θ = 0 . , σ = 1 . , ρ = − . , v = 0 .
1, along witha forward price at valuation time F = 1 .
0. Those parameters are not extreme, and are in thetypical range of a Heston fits to market option prices.In Figure 1, we look at the absolute error in price for a scale m = 8 and a truncation L = 12 based on the Heston cumulants. This truncation corresponds to an interval [ a, b ] =[ − . , . J , so that the overall error is dominatedby error in the payoff formula. We stop at strike K = 1 .
32 since then ln KF > b . Figure 1Figure 1.: Error in Vanilla option prices of maturity 2 days with Heston parameters κ = 1 . , θ = 0 . , σ = 1 . , ρ = − . , v = 0 . , F = 1 . L = 12 andscale m = 8. Strike ab s ( E rr o r) Method
ClassicNew shows that the error of the new formula stays below 10 − , close to machine epsilon while the ay 28, 2020 0:34 swift˙notes Draft error of the classic formula can be as high as 1 . · − when the strike approaches the upperboundary F e b . c m,k and quadratures We consider the same Heston model parameters as in the previous section. The trapezoidalrule is three to six times more accurate than the mid-point rule (or equivalently the formulafrom Ortiz-Gracia and Oosterlee (2016) based Vieta’s formula) across strikes and on bothHeston sets. Both rules use exactly the same number of points.Table 5.: Price of an out-of-the-money option under two different Heston parameter sets.Heston Method Strike Price ErrorSet 1 ( m = 8 , J = 12) Midpoint 250000 114.51 -4.87Trapezoidal 250000 117.91 -1.47Midpoint 4000000 3866.59 -85.33Trapezoidal 4000000 3931.09 -20.82Set 2 ( m = 6 , J = 5) Midpoint 1.0064 0.0063611 3.97e-07Trapezoidal 1.0064 0.0063606 -7.39e-08Midpoint 1.064 4.77e-06 5.09e-07Trapezoidal 1.064 4.18e-06 -8.22e-08We now look at the time to initialize the SWIFT method for a given option maturity. Thiscorresponds to the calculation of the density coefficients c m,k , either with the FFT appliedon the trapezoidal quadrature, or with the direct adaptive Filon quadrature on a relativetolerance of 10 − (which leads to a similar accuracy as the FFT approach). For a similarTable 6.: Initialization time of the SWIFT method for two different Heston parameter setsand different quadratures.Heston Method Points Time (microseconds)Set 1 FFT 4096 433Filon 585 76000Set 2 FFT 32 16Filon 497 588accuracy, the initialization based on the adaptive Filon quadrature is slower by a factor ofmore than 32 although the characteristic function is evaluated 585 times compared to 2048times for the FFT calculation. There is then a lot of room if we were to make the FFT densitycalculation adaptive by doubling successively the interval [ k , k ].
9. Conclusion
The use of the fast Fourier transform (FFT) to compute the payoff coefficients is particu-larly important and makes the SWIFT method competitive with some of the fastest pricingmethods such as COS method of Fang and Oosterlee (2008). Our alternative formula centeredon the forward is more accurate in general than the original payoff coefficients formula fromOrtiz-Gracia and Oosterlee (2016) while being of equivalent computational cost. ay 28, 2020 0:34 swift˙notes REFERENCES
The calculation of the density coefficients also benefits from the FFT, even though the re-lated characteristic function is relatively expensive to compute. The FFT based on the trape-zoidal rule is much more accurate than the original formula from Ortiz-Gracia and Oosterlee(2016) for a slightly lower computational cost. Using more fancy adaptive quadratures is no souseful. A simple adaptive scheme based successively doubling the truncation interval [ k , k ]according to the accuracy of the area underneath the curve is good enough. References
Amos, D. E. (1990) Algorithms 683: a portable FORTRAN subroutine for exponential integrals of a complex argument,
ACM Transactions on Mathematical Software (TOMS) , 16(2), pp. 178–182.Andersen, L. B. and Piterbarg, V. V. (2010)
Interest Rate Modeling, Volume I: Foundations and Vanilla Models , (AtlanticFinancial Press London).Fang, F. and Oosterlee, C. W. (2008) A novel pricing method for European options based on Fourier-cosine seriesexpansions,
SIAM Journal on Scientific Computing , 31(2), pp. 826–848.Jin, J. and Jjie, Z. S. (1996)
Computation of special functions , (Wiley).Johnson, S. G., Numerical integration and the redemption of the trapezoidal rule. (2011) , Technical report, MIT AppliedMath.K¨olbig, K., Exponential Integral for Complex Argument. (1990) , Technical report, CERN.Le Floc’h, F. (2013) Fourier Integration and Stochastic Volatility Calibration,
Available at SSRN 2362968 .Lord, R. and Kahl, C. (2007) Optimal Fourier inversion in semi-analytical option pricing,
SSRN pa-pers.ssrn.com/abstract=921336 .MacLeod, A. J. (1996) Rational approximations, software and test methods for sine and cosine integrals,
NumericalAlgorithms , 12(2), pp. 259–272.Makhoul, J. (1980) A fast cosine transform in one and two dimensions,
IEEE Transactions on Acoustics, Speech, andSignal Processing , 28(1), pp. 27–34.Maree, S. C., Ortiz-Gracia, L. and Oosterlee, C. W. (2017) Pricing early-exercise and discrete barrier options by Shannonwavelet expansions,
Numerische Mathematik , 136(4), pp. 1035–1070.Ortiz-Gracia, L. and Oosterlee, C. W. (2016) A highly efficient Shannon wavelet inverse Fourier technique for pricingEuropean options,
SIAM Journal on Scientific Computing , 38(1), pp. B118–B143.Pegoraro, V. and Slusallek, P. (2011) On the evaluation of the complex-valued exponential integral,
Journal of Graphics,GPU, and Game Tools , 15(3), pp. 183–198.Trefethen, L. N. and Weideman, J. (2014) The exponentially convergent trapezoidal rule,
SIAM Review , 56(3), pp.385–458.
Appendix A. Computing the discrete Cosine and Sine transforms together from theFFT
The calculation of the V m,k by the formula described in Appendix A ofOrtiz-Gracia and Oosterlee (2016) is the sum of a type 2 discrete cosine transform(DCT) and a type 2 discrete sine transform (DST). It can be summarized by the followingequation V m,k = N − X j =0 a j cos πk j + N ! + b j sin πk j + N ! (A1)with N = 2 ¯ J − for some positive integer ¯ J . Makhoul (1980) gives a simple algorithm tocompute the DCT of size N with one FFT of size N . We simply initialize the FFT coefficients c j with: c j = a j , c N − − j = a j +1 for j = 0 , ..., N − c , the DCT coefficients ˆ a areˆ a k = ℜ h ˆ c j e − iπ k N i (A3)Makhoul does not specify the equivalent formula for the DST, but we can do something ay 28, 2020 0:34 swift˙notes REFERENCES similar. We first initialize the FFT coefficients c j with: c j = b j , c N − − j = − b j +1 for j = 0 , ..., N − c , the DST coefficients ˆ b areˆ b k = −ℑ h ˆ c j e − iπ k N ii