A High Order Method for Pricing of Financial Derivatives using Radial Basis Function generated Finite Differences
AA High Order Method for Pricing of FinancialDerivatives using Radial Basis Functiongenerated Finite Differences
Slobodan Milovanovi´c and Lina von Sydow
Department of Information TechnologyUppsala UniversitySweden
Abstract
In this paper, we consider the numerical pricing of financial deriva-tives using Radial Basis Function generated Finite Differences in space.Such discretization methods have the advantage of not requiring Carte-sian grids. Instead, the nodes can be placed with higher density in areaswhere there is a need for higher accuracy. Still, the discretization matrixis fairly sparse. As a model problem, we consider the pricing of Europeanoptions in 2D. Since such options have a discontinuity in the first deriva-tive of the payoff function which prohibits high order convergence, wesmooth this function using an established technique for Cartesian grids.Numerical experiments show that we acquire a fourth order scheme inspace, both for the uniform and the nonuniform node layouts that we use.The high order method with the nonuniform node layout achieves veryhigh accuracy with relatively few nodes. This renders the potential forsolving pricing problems in higher spatial dimensions since the computa-tional memory and time demand become much smaller with this methodcompared to standard techniques.
Keywords:
Pricing of Financial Derivatives; Radial Basis Function generatedFinite Differences; High Order Methods; Smoothing of Initial Data.
In this paper, we are concerned with the numerical pricing of financial deriva-tives. A financial derivative is a contract whose value depends on an underlyingasset such as a stock, an interest rate, or a commodity. The trading in financial1 a r X i v : . [ q -f i n . C P ] A ug erivatives has increased tremendously during the last decades, mostly due tothe possibility to hedge positions in the underlying asset. Another importantfeature of financial derivatives is the potential for leverage, since a small move-ment in the underlying asset can cause a large movement in the value of thefinancial derivative.Due to the large traded volume of financial derivatives, efficient and accu-rate pricing of such contracts is of utmost importance. In most cases, there isno analytical formula available, and it becomes necessary to use a numericalmethod to compute the prices of the contracts. When a financial derivativedepends on several underlying assets, the problem becomes multi-dimensional.Traditionally, the only way to price such financial derivatives is to use MonteCarlo methods for a stochastic differential equation (SDE) formulation of theproblem. However, due to their arguably slow convergence, considerable effortsin the research community have been devoted to deriving efficient methods fora partial differential equation (PDE) formulation of the pricing problem. Themain problem with these methods is the so-called curse of dimensionality — thenumber of degrees of freedom in the problem grows exponentially in the numberof dimensions.Numerical methods for the PDE formulation include adaptive Finite Differ-ences (FD) [15, 12, 11, 16, 25], high-order compact schemes [3, 4], AlternatingDirection Implicit (ADI) schemes [8, 6], Radial Basis Function (RBF) approx-imation [17, 10], Radial Basis Function Partition of Unity (RBF-PU) method[20, 22, 21], and Radial Basis Function generated Finite Differences (RBF-FD)[14, 13]. In [24] and [26] several methods for pricing of options are implementedand evaluated.As a numerical example, we consider pricing of a European two-dimensionaloption. An option is a financial derivative which gives the holder the right,but not the obligation, to buy (for call options) or sell (for put options) anunderlying asset at a specified strike price K at or before the time of maturity T . The method that we employ is RBF-FD. The main idea behind it is tocombine the desirable features from FD (sparsity of the discretization matrices— as opposed to RBF) and RBF (meshfree — as opposed to FD). Such methodshave the potential to be of high order, depending on the number of nodes usedin the discretization stencil. However, for many option pricing problems, thepayoff function has a discontinuity in the function itself or its derivatives, whichlimits the order of convergence obtained in numerical simulations. For thisreason, we smooth the payoff function according to [9], before employing thenumerical method. This smoothing increases the order of convergence to theexpected one from the discretization that is used.In Section 2 we define the discretization in space and time, while Section 3 isdevoted to the model problems that we solve, as well as node layouts, stencils,2oundary conditions, and smoothing of the initial data. Finally, the results arepresented in Section 4, and conclusions are drawn in Section 5. We consider pricing of financial derivatives where the problem can be formulatedas a PDE in D spatial dimensions and time ∂u∂t + L u = 0 ,u ( s , . . . , s D ,
0) = g ( s , . . . , s D ) , (2.1) s i ≥ , i = 1 , . . . , D ; 0 ≤ t ≤ T. Here the solution u ( s , . . . , s D , t ) denotes the price of the financial derivative, t denotes time, s i , the value of the underlying asset with index i , and g the payofffunction of the financial derivative. In many pricing problems the original PDEis a final value problem solved backward in time. We consider problems inforward time as in (2.1), i.e., when necessary the problem is transformed intoan initial value problem.In Sections 2.1 and 2.2 we define the spatial and temporal discretization of(2.1), respectively. In RBF-FD the spatial operator L u in (2.1) at a location s c = ( s c , s c , . . . , s cD ),is approximated as a linear combination of the solution at the m closest nodes s k (possibly including s c ), k = 1 , . . . , m L u | s c ≈ m (cid:88) k =1 w k u | s k . (2.2)The weights w k are calculated by enforcing (2.2) to be exact for an RBF φ ( r ) φ ( (cid:107) s − s (cid:107) ) . . . φ ( (cid:107) s − s m (cid:107) )... . . . ... φ ( (cid:107) s m − s (cid:107) ) . . . φ ( (cid:107) s m − s m (cid:107) ) w ... w m = L φ ( (cid:107) s − s (cid:107) ) | s c ... L φ ( (cid:107) s − s m (cid:107) ) | s c . (2.3)It is given from RBF interpolation that (2.3) is a nonsingular system, and hencea unique set of weights w k , k = 1 , . . . , m can be computed.Typical choices of RBFs are listed in Table 1. For the first four examples,the parameter ε ∈ R is the shape parameter of the RBF. For polyharmonicsplines (PHSs), the parameter q ∈ N .In this paper, we follow [5], [1], and [13] and use PHSs as basis functionstogether with polynomials of degree p in the interpolation. With that approach,3he polynomial degree (instead of the RBF) controls the rate of convergence,while the RBFs contribute to reduction of approximation errors and are neces-sary in order to have a stable approximation. φ ( r )Gaussian e − ( εr ) Inverse quadratic 1 / (1 + ( εr ) )Multiquadric (cid:112) εr ) Inverse multiquadric 1 / (cid:112) εr ) Polyharmonic splines r q − Table 1: A list of commonly used RBFs φ ( r ).In (2.4), we augment (2.3) with monomials of degree one s . . . s D B s m . . . s mD . . . . . . s . . . s m . . . s D . . . s mD . . . w ... w m γ γ ... γ D = L φ ( (cid:107) s − s (cid:107) ) | s c ... L φ ( (cid:107) s − s m (cid:107) ) | s c L | s c L s | s c ... L s D | s c , (2.4)where B is the coefficient-matrix in (2.3).Now, we place N computational nodes s ci , i = 1 , . . . , N at the locations wherewe want to approximate the solution. The weights for each computational nodefrom solving (2.4) are assembled row-wise into the sparse differentiation matrix W ∈ R N × N , with m nonzero elements per row. This leads to the followingsemi-discretization of (2.1) dd t ¯ u ( t ) + W ¯ u ( t ) = ¯0 , ¯ u (0) = ¯ g, (2.5)where ¯ u ( t ) ∈ R N × is the vector of unknowns at time t , with approximationsof u in the computational nodes s ci , i = 1 , . . . , N , ¯0 ∈ R N × is a vector withonly zeros, and ¯ g ∈ R N × is the vector with the function g evaluated in thecomputational nodes s ci , i = 1 , . . . , N .Equation (2.5) forms a system of linear ordinary differential equations (ODEs)in time. In the next section, we describe how to solve it.4 .2 Temporal discretization For the time discretization of (2.5), we use the Backward Differentiation Formulaof order two (BDF2) [7]. This time-stepping scheme requires the solution at twoprevious time steps, and we therefore employ Backward Euler (BDF1) for thefirst time step. It is convenient to have the same coefficient matrix in all timesteps, so we use non-equidistant time steps as described in [10] and later usedin e.g., [14, 13]. This is accomplished by discretizing the time interval with M steps of length ∆ t (cid:96) = t (cid:96) − t (cid:96) − , where (cid:96) = 1 , . . . , M . We define ω (cid:96) = ∆ t (cid:96) / ∆ t (cid:96) − for (cid:96) = 2 , . . . , M and arrive at¯ u − ¯ u = ∆ t W ¯ u , (2.6)¯ u (cid:96) − β (cid:96) ¯ u (cid:96) − + β (cid:96) ¯ u (cid:96) − = β (cid:96) W ¯ u (cid:96) , (cid:96) = 2 , . . . , M, (2.7)where β (cid:96) = ∆ t (cid:96) ω (cid:96) ω (cid:96) , β (cid:96) = (1 + ω (cid:96) ) ω (cid:96) , β (cid:96) = ω (cid:96) ω (cid:96) . (2.8)We compute the values for ω (cid:96) using the recursive condition β (cid:96) = β (cid:96) − , whichkeeps the coefficient matrix constant throughout all time steps. Since our timeinterval has the length T , we chose the initial time step length ∆ t from M (cid:88) (cid:96) =1 ∆ t (cid:96) = T = ∆ t (1 + M (cid:88) (cid:96) =2 (cid:96) (cid:89) (cid:96) (cid:48) =2 ω (cid:96) ) . (2.9)Finally, we start the time integration by setting ¯ u = ¯ g .From the temporal discretization we get the following linear system of equa-tions to solve in each time step A ¯ u (cid:96) = ¯ b (cid:96) , (2.10)where A = I − ∆ t W , ∆ t is given by (2.9), ¯ b (cid:96) = β (cid:96) ¯ u (cid:96) − − β (cid:96) ¯ u (cid:96) − , (cid:96) = 2 , . . . , M ,and ¯ b = ¯ u = ¯ g . As a model problem we consider a European call option issued on two underlyingassets s and s ∂u∂t + L u = 0 , (3.11) s ≥ , s ≥ , ≤ t ≤ T, with L u = 12 (cid:16) σ s ∂ u∂s + σ s ∂ u∂s (cid:17) + ρσ σ s s ∂ u∂s ∂s + r (cid:16) s ∂u∂s + s ∂u∂s (cid:17) − ru, (3.12)5nd u ( s , s , T ) = g ( s , s ) = (cid:18)
12 ( s + s ) − K (cid:19) + . (3.13)Here ( f ( x )) + = max( f ( x ) , r denotes the risk-free interest rate in the market, σ i the volatility of asset i , and ρ the correlation between the assets. As a close-field boundary condition in s = s = 0 we set u (0 , , t ) = 0 , ≤ t ≤ T, (3.14)and as a far-field boundary condition we set u ( s , s , t ) = (cid:18)
12 ( s + s ) − Ke − rt (cid:19) , ≤ t ≤ T, (3.15)for s and s large enough. The parameters used are given in Table 2. r σ σ ρ , K T We consider both a uniform and a nonuniform node layout, presented in Figure1. Unlike classical grid-based methods (e.g., standard FD methods) we do notneed to use a rectangular domain. Instead, we only use the lower-triangular halfof the rectangle which reduces the number of computational nodes by a factorof two, and hence the computational complexity significantly.The reason for introducing a nonuniform node layout is that we can clusternodes where we are most interested in having an accurate solution. In general,we are most interested in having an accurate solution in the neighborhood of s + s = 2 K , which is also where the truncation error is largest due to largederivatives in the solution from the discontinuity in the first derivative of thepayoff function.We start to present a nonuniform node distribution in 1D that is generatedas introduced in [8] and later used for RBF-FD and option pricing in [14].6onsider N equidistant nodes x (1)1 < . . . < x ( i )1 < . . . < x ( N )1 constructed by x ( i )1 = arcsinh (cid:18) − Kc (cid:19) + ( i − x , i = 1 , . . . , N , (3.16)where c is a positive real constant which specifies how dense the node distribu-tion becomes around the strike price K ,∆ x = 1 N (cid:20) arcsinh (cid:18) s max − Kc (cid:19) − arcsinh (cid:18) − Kc (cid:19)(cid:21) ,and s max denotes the far-field boundary. Then, the nonuniform node distribu-tion s is generated pointwise as s ( i )1 = K + c · sinh( x ( i ) ), i = 1 , . . . , N . (3.17)The nonuniform node layout is generated by using the one-dimensional nodelayouts from (3.16) and (3.17), along the axes s and s , and then uniformlyplacing the internal points in the diagonal direction. The number of nodes alongeach diagonal is increased by one for each diagonal. The far-field boundary islocated at s + s = s max = 8 K . The density tuning parameter used in Figure 1for the nonuniform node layout and in the numerical experiments presented inSection 4, is c = 0 .
8. It should be noted that a too small value of c eventuallyleads to an ill-conditioned problem.0 K K K K K K s s (a) Uniform node layout. K K K K K K s s (b) Nonuniform node layout. Figure 1: Uniform and nonuniform computational node layouts in 2D. Theboundary conditions are employed in the blue triangle node (the close-fieldboundary condition) and in the red square nodes (the far-field boundary condi-tion).We also introduce the notation N s for the number of nodes along one of the7xes, i.e., N s ( N s + 1)2 = N. (3.18)The nearest neighbors for constructing the stencils are efficiently determinedusing the k -D tree algorithm, [2]. In Figure 2 we show examples of stencils atdifferent locations in the domain. The polynomial space is of size ν = (cid:32) p + Dp (cid:33) , which we use to set the size of the stencils to m = 5 ν , following [5, 1, 13]. Weare aiming for a fourth order scheme and use p = 4 and q = 5 which gives ν = (cid:32) (cid:33) = 15 . Hence, we use a stencil size that is m = 5 ·
15 = 75. s s Figure 2: Examples of nearest neighbor based stencils used for approximatingthe differential operator on a nonuniform node layout. The central node of eachdisplayed stencil is denoted by a white cross mark. All stencils are of the samesize m = 75.For the boundary nodes we use different treatments depending on where thenode is located. For the node s = s = 0 (the blue triangle in Figure 1), we setthe close-field boundary condition from (3.14). For the nodes s + s = 8 K (the8ed squares in Figure 1), we set the far-field boundary conditions from (3.15).For the boundary nodes along the axes, i.e., s = 0, s > s = 0, s > k -Dtree algorithm generates one-sided stencils for those nodes. Since the initial data g ( s , s ) in (3.13) has a discontinuity in the first derivative,the obtained spatial order of convergence for a finite difference scheme is limitedto two, regardless of the formal order of the scheme. The formal spatial orderof the scheme in (2.4) is p , i.e., for p > defined by its Fourier transformˆΦ ( ω ) = (cid:18) sin( ω/ ω/ (cid:19) (cid:18) ( ω/ (cid:19) . (3.19)Using Wolfram Mathematica to compute the inverse Fourier transform of(3.19) givesΦ ( s ) = 172 (cid:16) − ( s − · sgn( s − − ( s + 3) · sgn( s + 3)+ 12( s − · sgn( s −
2) + 12( s + 2) · sgn( s + 2) − s − · sgn( s − − s + 1) · sgn( s + 1)+ 56 s · sgn( s ) (cid:17) , (3.20)where sgn( x ) = | x | x . Following [9], [18], and [4], we get the smoothed final condition on a uniformnode layout as˜ g ( s , s ) = 1∆ s s (cid:90) − s s (cid:90) − s Φ (cid:18) ˜ s ∆ s (cid:19) Φ (cid:18) ˜ s ∆ s (cid:19) g ( s − ˜ s , s − ˜ s ) d˜ s d˜ s . (3.21)Since g ( s , s ) is smooth in a large part of the computational domain, we onlyneed to compute (3.21) in the nodes that are close enough to s + s = 2 K to be affected from the smoothing. Also, since the nodes along a diagonal allhave the same distance to s + s = 2 K , we only need to compute one value of˜ g ( s , s ) for each diagonal and use that value for all nodes on that diagonal.9he theory in [9] shows that replacing the final condition g ( s , s ) with˜ g ( s , s ) defined in (3.21) gives a fourth order scheme for Cartesian grids, i.e.,the node layout that we here refer to as uniform. Here, we want to use thissmoothing also for our nonuniform node layout defined in Section 3.1. Thislayout can be seen as a slightly skewed Cartesian grid and the nodes are equidis-tantly distributed along the diagonals. For this node layout we replace ∆ s in(3.21) with ∆ s i = k (cid:54) = c min k =1 ,...,m (cid:107) s ci − s ki (cid:107) , i = 1 , . . . , N. The numerical method described in Section 2 applied to the model problemsdescribed in Section 3 is implemented in
Matlab . In all experiments, we startby scaling the original problem such that s max = 1 and time runs forward inthe PDE. After the integration, the solution is transformed back to the originalproblem.The linear system defined in (2.10) is solved using GMRES [19], with anincomplete LU factorization as the preconditioner using nofill . The conver-gence tolerance for the iterations is set to 10 − , and as the initial condition foreach iteration we use the computed solution from the previous time step.The numerical experiments are performed on a laptop equipped with a 2.3GHz Intel Core i7 CPU and 16 GB of RAM. The computation of the RBF-FDweights is performed in parallel using the parallel toolbox command parfor with four workers.In Figure 3 we plot the error ∆ u max as a function of ˆ h ≡ / √ N as well asof CPU-time for the model problem. The error is defined as∆ u ( s , s ) = | u c ( s , s , − u ∗ ( s , s , | , (4.22)where u c is the computed solution and u ∗ is a reference solution computed witha second order finite difference method on a very fine grid. We use (4.22) todefine ∆ u max = max [ s ,s ] ∈ ˆΩ ∆ u ( s , s ) , (4.23)where ˆΩ = (cid:2) K, K (cid:3) × (cid:2) K, K (cid:3) . We denote standard second order finitedifferences [23] by FD. RBF-FD-GS is an RBF-FD method with Gaussian RBFs,stencil size m = 25 and a node density dependent shape parameter that ispresented in detail in [14]. Abbreviation RBF-FD-PHS is used for the methodthat is presented in this paper. Moreover, we use designation smoothed in thesuperscript for the computations performed with the smoothing of the initialdata, and uniform and nonuniform to specify the node layouts.10ndependent of spatial discretization that is used, we employ BDF2 with M = N s time steps in all experiments. For the RBF-FD methods N s is definedin (3.18) and for FD it is defined by N s = √ N . With this number of time steps,the temporal discretization error is not visible in the plots.10 − − − − − − ˆ h ∆ u m a x Convergence O (ˆ h ) RBF-FD-GS uniform RBF-FD-PHS nonuniform O (ˆ h ) RBF-FD-GS nonuniform RBF-FD-PHS smootheduniform
FD RBF-FD-PHS uniform
RBF-FD-PHS smoothednonuniform − timePerformanceFigure 3: ∆ u max as a function of ˆ h and CPU-time in seconds for the Europeancall option.In Figure 3 we see that all methods, but the two that are using a smoothedfinal condition, exhibit second order convergence. Among those five secondorder methods, RBF-FD with PHS exhibits the smallest error for a given N ,and a nonuniform node layout gives a smaller error than the uniform one usingthe same N . RBF-FD-PHS with a smoothed final condition exhibits fourthorder spatial convergence whether we are using a uniform or nonuniform nodelayout (apart from a small deviation for the nonuniform node layout). Whenit comes to computational time to reach a certain ∆ u max , FD is competitivefor the larger errors displayed. This makes sense since the RBF-FD methodsall have to compute the weights w k , k = 1 , . . . , m before the time-stepping.Moreover, our model problem has a fairly short time to maturity T = 0 .
2. Forlonger times to maturity, FD does not perform equally well compared to theRBF methods, see [14, 13]. We also establish that the fourth order methods11uickly become superior when it comes to CPU-time to reach a certain ∆ u max .That is especially true for RBF-FD-PHS smoothednonuniform . Even though this methodhas a computational prephase that includes both computation of weights w k , k = 1 , . . . , m and smoothing of the final condition, the method requires a muchsmaller CPU-time than the other methods for ∆ u max < − .0 1 2 3 401234 s s RBF-FD-PHS uniform s RBF-FD-PHS smootheduniform − − − ∆ u Figure 4: Heat maps of ∆ u for the European call basket option on uniform nodelayouts. The boundary of ˆΩ is marked with a white dash-dotted line.0 1 2 3 401234 s s RBF-FD-PHS nonuniform s RBF-FD-PHS smoothednonuniform − − − ∆ u Figure 5: Heat maps of ∆ u for the European call basket option on nonuniformnode layouts. The boundary of ˆΩ is marked with a white dash-dotted line.12n Figures 4–5 we show a heat map of the error ∆ u defined in (4.22) forthe model problem using the method RBF-FD-PHS. In Figure 4 we present theerror on the uniform node layout, and in Figure 5 the error on the nonuniformnode layout, with N = 6105 for both node layouts. To the right in both figures,we have used the smoothed final condition, while the original one is used inthe plots on the left. The errors in Figures 4–5 are presented for 0 ≤ s j ≤ j = 1 ,
2, in order to have a better view of the error profile around the smoothedarea. The color scale is the same in all four plots.From Figures 4–5 we conclude that the smoothing of the final conditionrenders a ∆ u that has a smaller magnitude compared to the original final con-dition. Moreover, ∆ u obtained from the smoothed final condition has threemaxima along a line s = s + const. while the corresponding number of max-ima is 1 for the nonsmoothed final condition. We also note that the magnitudeof ∆ u is smaller for the nonuniform node layout compared to the uniform one.That is due to the fact that for the nonuniform node layout, the number of nodesis larger in the area where the solution has large derivatives, i.e., around thestrike price. We end this section by concluding that in this particular example,∆ u max is more than one order smaller using smoothing of the final condition ona nonuniform grid than without the smoothing on a uniform grid for the samenumber of nodes. In this paper, we have implemented a solver to price financial derivatives basedon RBF-FD discretization in space and BDF2 in time. As RBFs we use PHSs,augmented with monomials of up to degree p . The formal order of this spatialdiscretization is p , however for many pricing problems the lack of smoothnessof the initial data limits the actual order obtained in numerical simulations.However, by employing a smoothing technique to the initial data, the formalorder of the discretization is retained.The RBF-FD discretizations have the advantage over standard FD suchthat the nodes do not have to be organized in a Cartesian grid. On the otherhand, the RBF-FD discretizations have the merit that they render sparse dif-ferentiation matrices as opposed to global RBF approximations that lead tofull matrices. Thus, RBF-FD has the possibility to give accurate solutions onnonuniform node layouts, still yielding sparse matrices.As a model problem, we consider pricing of a European type basket optionissued on two underlying assets, resulting in a PDE in two spatial dimensionsand time. By employing a nonuniform node layout that has a denser node dis-tribution where we are most interested in having an accurate solution, togetherwith smoothing of the final condition, the numerical experiments demonstrate13hat our developed method gives a very accurate solution in a short time usingfewer nodes than the methods that we compare with, for this model problem.The fact that we can solve the problem accurately with fewer nodes becomesextremely important when we want to solve problems in higher dimensions, e.g.,for pricing of financial derivatives issued on several underlying assets. Since thenumber of degrees of freedom grows exponentially in the number of dimensions(number of underlying assets), the ability to use fewer nodes per dimension toreach a certain accuracy might lead to the possibility to solve problems thatwould not be possible to solve with traditional techniques. References [1] V. Bayona, N. Flyer, B. Fornberg, and G. A. Barnett. On the role ofpolynomials in RBF–FD approximations: II. Numerical solution of ellipticpdes.
Journal of Computational Physics , 332:257–273, 2017.[2] J. L. Bentley. Multidimensional binary search trees used for associativesearching.
Communications of the ACM , 18(9):509–517, 1975.[3] B. D¨uring and M. Fourni´e. High-order compact finite difference scheme foroption pricing in stochastic volatility models.
Journal of Computationaland Applied Mathematics , 236(17):4462–4473, 2012.[4] B. D¨uring and C. Heuer. High-order compact schemes for parabolic prob-lems with mixed derivatives in multiple space dimensions.
SIAM Journalon Numerical Analysis , 53(5):2113–2134, 2015.[5] N. Flyer, B. Fornberg, V. Bayona, and G. A. Barnett. On the role ofpolynomials in RBF–FD approximations: I. Interpolation and accuracy.
Journal of Computational Physics , 321:21–38, 2016.[6] T. Haentjens and K. J. In ’t Hout. Alternating direction implicit finitedifference schemes for the Heston–Hull–White partial differential equation.
Journal of Computational Finance , 16:83–110, 2012.[7] E. Hairer, S. Nørsett, and G. Wanner.
Solving Ordinary DifferentialEquations I: Nonstiff Problems . Solving Ordinary Differential Equations.Springer, 2008. ISBN 9783540788621.[8] K. In’t Hout and S. Foulon. ADI finite difference schemes for option pricingin the Heston model with correlation.
Int. J. Numer. Anal. Model , 7(2):303–320, 2010. 149] H.-O. Kreiss, V. Thome, and O. Widlund. Smoothing of initial data andrates of convergence for parabolic difference equations.
Communicationson Pure and Applied Mathematics , 23(2):241–259, 1970.[10] E. Larsson, K. ˚Ahlander, and A. Hall. Multi-dimensional option pricingusing radial basis functions and the generalized Fourier transform.
Journalof Computational and Applied Mathematics , 222(1):175–192, 2008.[11] G. Linde, J. Persson, and L. von Sydow. A highly accurate adaptive finitedifference solver for the Black–Scholes equation.
International Journal ofComputer Mathematics , 86(12):2104–2121, 2009.[12] P. L¨otstedt, J. Persson, L. von Sydow, and J. Tysk. Space–time adaptivefinite difference method for European multi-asset options.
Computers &Mathematics with Applications , 53(8):1159–1180, 2007.[13] S. Milovanovi´c. Pricing financial derivatives using radial basis functiongenerated finite differences with polyharmonic splines on smoothly varyingnode layouts. arXiv preprint, arXiv:1808.02365[q-fin.CP] , 2018.[14] S. Milovanovi´c and L. von Sydow. Radial basis function generated finitedifferences for option pricing problems.
Computers & Mathematics withApplications , 75(4):1462–1481, 2018.[15] J. Persson and L. von Sydow. Pricing European multi-asset options using aspace-time adaptive FD-method.
Computing and Visualization in Science ,10(4):173–183, 2007.[16] J. Persson and L. von Sydow. Pricing American options using a space-time adaptive finite difference method.
Mathematics and Computers inSimulation , 80(9):1922–1935, 2010.[17] U. Pettersson, E. Larsson, G. Marcusson, and J. Persson. Improved radialbasis function methods for multi-dimensional option pricing.
Journal ofComputational and Applied Mathematics , 222(1):82–93, 2008.[18] D. M. Pooley, K. R. Vetzal, and P. A. Forsyth. Convergence remedies fornon-smooth payoffs in option pricing.
Journal of Computational Finance ,6(4):25–40, 2003.[19] Y. Saad and M. H. Schultz. GMRES: A generalized minimal residual algo-rithm for solving nonsymmetric linear systems.
SIAM Journal on scientificand statistical computing , 7(3):856–869, 1986.[20] A. Safdari-Vaighani, A. Heryudono, and E. Larsson. A radial basis functionpartition of unity collocation method for convection–diffusion equations15rising in financial applications.
Journal of Scientific Computing , 64(2):341–367, 2015.[21] V. Shcherbakov. Radial basis function partition of unity operator splittingmethod for pricing multi-asset american options.
BIT Numerical Mathe-matics , 56(4):1401–1423, 2016.[22] V. Shcherbakov and E. Larsson. Radial basis function partition of unitymethods for pricing vanilla basket options.
Computers & Mathematics withApplications , 71(1):185–200, 2016.[23] D. Tavella and C. Randall.
Pricing Financial Instruments: The FiniteDifference Method (Wiley Series in Financial Engineering) . Wiley NewYork, 2000.[24] L. von Sydow, L. Josef H¨o¨ok, E. Larsson, E. Lindstr¨om, S. Milovanovi´c,J. Persson, V. Shcherbakov, Y. Shpolyanskiy, S. Sir´en, J. Toivanen, et al.BENCHOP–the BENCHmarking project in Option Pricing.
InternationalJournal of Computer Mathematics , 92(12):2361–2379, 2015.[25] L. von Sydow, J. Toivanen, and C. Zhang. Adaptive finite differences andIMEX time-stepping to price options under Bates model.