Spreads contracts on commodities are the difference between the finished product and the raw material used in its production. Examples of these contracts are cracks, based on the difference between the price of gasoline and oil, sparks, based on the difference between electricity and natural gas prices and crushs based on the difference between soy oil and soy beans. Typically, manufactures look at the difference between a long position on the raw material and a short position on the resulting product.
The goal of this paper is to investigate a pricing approach following the numerical solution of the PIDE associated to option prices driven by two-dimensional Levy models, using finite elements. Its original contribution consists in an efficient strategy that combines a Galerkin discretization in space via the symbol method and BICSTAB along with a circular preconditioner to efficiently solve the associated linear system with non-symmetric and densely populated BTTB matrix involved in the implicit time step algorithm. Additionally, we present a practical application of this strategy to specific crack contracts. The paper provides numerical experiments that validate the good performance of this strategy within the framework of the proposed models.
Levy models offer a more realistic approach to the dynamic of the underlying assets involved in the contract as they allow to capture rapid changes in their movement via random jumps. We propose two types of bivariate stochastic models. The first one is based on a combination of a diffusion and random jumps. In turn, the jumps are divided into idiosyncratic and common ones, each with its specific intensity and length. It corresponds to market information that may impact both or each one separately. The model has been introduced in [2] . The second model belongs to the subclass of time-changed Levy models, where the time index is consider as a non-decreasing process called subordinator. Its extension to a multivariate setting has been consider in [3, 4] . In both models the dependence between assets is captured in a tractable way by the bivariate Brownian motion and the common jump length in the first case and an additional common subordinator in the second.
For exponential Levy models or a jump-diffusion models, the value of an European option can be obtained as the classical solution of the related Kolmogorov equation where the corresponding operator is of partial integro-differential type. For this notion of solution, Fast Fourier Transform (FFT) based methods have been proposed in [5, 6] which are quite easy to implement, flexible and generic although they do not handle boundary conditions efficiently and are only applicable to linear problems with non space-dependent coefficients. Other possible notions of solution are viscosity solutions [7, 8, 9, 10] to be approximated using finite difference methods and variational solutions [1, 11, 12] to be approximated using variational methods. In [13, 14] , important theoretical advances are presented to establish the connection between variational solutions and the corresponding expectations that European option prices represent. The resulting relation is a Feynman-Kac-type representation of the variational solution as a conditional expectation providing a theoretical framework that covers a wide range of Levy models. With these results, the use of the Galerkin method to approximate variational solutions acquires special relevance.
Important features of the Galerkin method are its solid theoretical background an its flexibility to choose approximation spaces and basis to address the particularities of the problem to be solved. For example, specific finite element spaces can be constructed to handle complicated domain geometries, non trivial boundary conditions or to capture different degrees smoothness of the solution. Even so, when applying the Galerkin method to PIDEs some difficulties appear. The matrix associated to the operator is full, which involves two fundamental questions: how to construct it? and (if implicit methods for time evolution are to be considered) how to solve the companion linear system ?. High order algorithms for time evolution are highly desirable since computational effort is substantially lower. By far the favorite approach in pricing has been the use of θ-schemes. They are easy to implement and offer flexibility to obtain implicit algorithms allowing for larger stability regions and thus larger time steps. In particular the choice θ = 1/2 leads to the Crank-Nicholson scheme which is second order accurate and unconditionally stable. Some wavelets-Galerkin methods has been proposed to avoid dealing with full matrices [15, 16, 17] . By using this methods, it is possible to construct a compressed sparse matrix to substitute the original full one. The compressed matrix is also structured and only a few of its entries need to be computed leading to tractable linear systems.
Recently, Gass and Glau presented in [1] a flexible Galerkin scheme for option pricing in multidimensional Levy models, the authors developed an approach to set up a FEM solver for option prices in Levy models that they call it the symbol method. They show a way to implement a Galerkin discretization of the integrodifferential operator whenever its symbol is known and the approximated solution is approximated as the linear combination of the smooth mother basis function's displacements. By taking this approach the stiffness matrix in one dimension is Toeplitz and only its first line and first column are needed to determine it. There exist direct methods to solve a Toeplitz system in O(N log 2 N ) operations [18, 19, 19, 20, 21, 22] and iterative methods to solve it in O(N log(N )) operations [23, 24] . Implementing the symbol method in two dimensions lead us to a BTTB stiffness matrix, an introduction to iterative solvers for Toeplitz and BTTB symmetric systems can be found in [25] .
Inspired by the use of the Preconditioned Conjugated gradient (PCG) method for symmetric Toeplitz and BTTB systems covered in [25] we propose the use of the Preconditioned BICGSATB (Krilov) method introduced in [26] , which works for general (non-symmetric) matrices, along with a circular preconditioner. As mentioned in [1] , the pricing error is sensitive to the accuracy of the stiffness matrix entries and the authors emphasized the necessity of accurately compute them. Pointing in that direction, we discretize the terms in the integro-differential operator not related to the pure jump using exact Gaussian quadrature (as it is suggested in [27] for the case of classical FEM) and suppress their corresponding terms in the symbol. By doing so, some amount of error is avoided since contribution of the diffusion and the drift to the stiffness matrix will be exact.
The paper is organized in the following way: In Section 2 we introduce the main notations and the models to be considered based on two subclasses of exponential Levy models, together with their symbols and their associated PIDEs. Section 3 is devoted to Galerkin method for parabolic problems, the variational formulation and discretization of the PIDE will be described as well as the iterative method to solve the BTTB system and the θ-scheme. In Section 4 we discuss the specific implementation of the method for the two proposed models, estimation problems and the numerical pricing results for crack contracts.
For a Levy process (X t ) t≥0 let the functions ϕ Xt (u) and Ψ X (u) = 1 t log ϕ Xt (u) be their characteristic function and characteristic exponent respectively. Random elements are defined in a filtered space (Ω, A, P, (F t )) t≥0 ) completed in the usual way. For row vectors a and b the product ab refers to its component-wise product while scalar product is denoted a T b, where a T is its transpose. The symbol C 1,2 (R 2 ) denote respectively the set of functions with first and second order continuous derivatives . We define a + := max(a, 0). (Triplet notation) The triplet of a Levy process is denoted by (b, Σ, ν) where b , Σ are the drift and the diffusion coefficients respectively, while ν is the Levy measure
t ) t≥0 be a two dimensional stochastic process whose components represent the prices of two underlying assets, e.g. a commodity and the raw material to produce it at time t, and let
t ) t≥0 be the corresponding log-return process. Both processes are related by
For simplicity we do not include a mean-reversion dynamic which has been reported in commodity movements, but the method above can be adapted to this more realistic situation.
The absence of arbitrage imposes the discounted underlying prices (e −rt S t ) t≥0 to be a Q-martingale in each dimension under a given equivalent martingale measure Q (EMM for short) which, for a Levy process, translates to restrictions on its triplet.
A European type spread call contract with maturity at time T and strike price K is a contingent claim that with a payoff given by:
where the value c > 0 is a conversion factor. Commonly, contracts are written on future prices on both assets. Again, for simplicity we focus on the case of spot prices. The price of a spread contract at time t with initial price S t = x is denoted as C(t, x) and is given by:
Next, we consider specific exponential Levy models that reflect finite and infinite jump activity in the logprices.
We look at a bivariate jump-diffusion process with common and idiosyncratic jumps. Such models have been studied in [2] by considering two Poisson compound processes with a dependence between random sizes in the common jumps. To this end, we define two sequences of independent and identically distributed two-dimensional random vectors (X k ) k∈N and (X 0,k ) k∈N . The components of the vectors in the first sequence are independent among them, with identical cumulative distribution function (c.d.f.) F X := (F (1) (x 1 ), F (2) (x 2 )) and independent also from the components of vectors X 0,k . On the other hand, the vectors X 0,k = (X (1) 0,k , X
0,k ) k∈N have joint cumulative distribution function F X0 (x) for every k ∈ N. Next, we introduce the bivariate homogeneous compound Poisson process (Z t ) t≥0 = (Z (1) t , Z (2) t ) t≥0 with components:
t ) t≥0 is a pair of independent Poisson processes with respective intensities λ j > 0, j = 1, 2, also independent of the Poisson process (N (0) t ) t≥0 , the latter with intensity λ 0 > 0. The processes (N k . We assume the existence of moments up to order two of the later. Furthermore, we set the log-prices as jump-diffusion processes given by a combination of a Brownian motion and a compound Poisson process as:
t ) is a pair of Brownian motions with covariance matrix Σ B with
The process (Y t ) t≥0 has triplet (b, Σ B , ν) with
and a Levy jump measure ν on R 2 given by:
The characteristic function of the log-price process is the product of the characteristic functions of processes Z t and B t . Indeed, for the compound Poisson part is easily seen that
where ϕ X0,j and ϕ X (j) are the respective characteristic functions of common and idiosyncratic jump sizes of the j-th asset. Hence:
Under Q, the characteristic exponent of the log-prices verifies Ψ Q Y (j) (−i) = r. Therefore, the drift changes to 
J ) and X 0,k ∼ N (µ 0,J , Σ 0,J ), such that:
Here δ jl is the usual Kronecker's number. In this case
By equation (6) the characteristic exponent of the log-prices is
where the drift across each dimension in the risk-neutral measure takes the form
We consider now the class of exponential bivariate time-changed Levy models, with the random time described according to a two dimensional subordinator of one-factor type. See [3, 4] for a general formulation of one-factor multidimensional time-changed processes. The dependence between both underlying assets is given through the common time-changed subordinator. In this context, we define the log-price process as:
t ) t≥0 is a vector of subordinator Levy models starting at zero, with independent components and finite second moments. As in the previous subsection the process (B t ) t≥0 is a bivariate standard Brownian motion but with independent components. The pairs µ, σ r = (σ 1,r , σ 2,r ) and d are parameters of the model.
Next, we provide an expression for the characteristic function of the log-price process below. It is a simple extension of well-known results in the one-dimensional case.
Proposition 2.1. Let the process (Y t ) t≥0 be given by equations (1) and (8) . For t ≥ 0 and u ∈ R 2 we have:
Proof. Let the quantity L (0) t represent either the one-dimensional process or the pair (L
But:
After substituting the expression above into equation (11) , from the independence of the subordinator processes, we have
which immediately leads to equations (9) and (10).
Moments of the log-return process (∆Y t ) t≥0 , needed for the parameter estimation in section 4, can be obtained by conditioning on the subordinator process. The result is stated below.
Proposition 2.2. Let the process (Y t ) t≥0 be given by equations (1)- (8) . For any t ≥ 0, k = 1, 2, 3, 4 and
where
Proof. We have the relations ∆Y
By the independence of the subordinators taking conditional expectations on both sides we have that
Hence we have (13) . For the third moment we have
Taking into account that
from where we get (14) . Similarly, the forth moment is:
where:
Gathering all term we get (15) . Finally, for the mixed moment
From the expressions:
we have (16) . 
Notice that :
Hey proposition 1:
The marginal characteristic exponents of (Y t ) t≥0 are obtained also from equat(10) as:
Hence, the martingale conditions Ψ (j) (−i) = r translates into the equations: 
Consider the interval I = [0, T ] ⊂ R, a terminal condition h : R 2 → R and a function C : I × R 2 → R. The Kolmogorov equation is of the form ∂ t C + AC = 0,
where A is the Kolmogorov operator of a Levy model and is given by
for all t ∈ I and x ∈ R 2 , with symbol
If we change variable τ = T − t and define
we have that u satisfies equation.
The Galerkin method is based in the variational formulation of problem (19) and our choice for an implicit θ-scheme requires the resolution of a BTTB system with non-symmetric and densely populated matrix at each time step.
Assuming that (19) has a classical solution u ∈ C 1,2 (R 2 ) (sufficiently regular to apply Ito's formula) while satisfying an appropriate integrability condition on h, there is a Feynman-Kac type representation for u(τ ) so that C(t) has the stochastic representation (2) . But, since we are considering a variational formulation of (19) , such hypothesis turns out to be very restrictive. We are looking for less smooth solutions. Let the Gelfand triplet (V, H, V * ) be composed by separable Hilbert spaces H, V and the dual space V * such that there exist a continuous embedding from V into H. Denote (·, ·) H the inner product of H and ·, · V * ×V the duality pairing. Let L 2 (I, H) be the space of weekly measurable functions u :
Weak solutions of (17) will belong to the Sobolev space
The variational formulation of problem (19) results now in: To find u ∈ W 1 (I, V, H) such that
where h ∈ H, convergence in the second line above is in the norm of H and
We adopt the theoretical framework in [14, 13, 1, 28] and consider the Sobolev-Slobodeckii spaces V = H α η (R 2 ) and H = L 2 η (R 2 ). For concrete definitions of weighted Sobolev-Slobodeckii spaces see [13, 14] . Two matters are to be considered, one related to existence and uniqueness of the solution of (21) and the other related to the equivalence between such a solution and (2) . Both problems has been addressed in [13] and [14] . Well posedness of problem (21) is guarantee by [13, Theorem 5.3 pag. 15] if hypotheses (A1), (A2) and (A3) are to be satisfied by the symbol A whenever the payoff h ∈ L 2 η (R 2 ).
(A3) There exist constants C2 > 0 and C3 ≤ 0 uniformly in time, such that for a certain 0 ≤ β < α 
) for problem (27) . Under the additional condition that h ∈ L 1 η (R 2 ) and as a consequence of [13, Theorem 6.1 pag. 17] we have that the solution u has the stochastic representation 2.
Remark 3.1. It doesn't exist an index α for Example 2.2 so we can not use the results above in this case. Even so, we applied the symbol method to this example and the numerical results in Section 4 show convergence evidence.
We separate the Kolmogorov operator in two parts, one grouping Diffusion, Drift and and the other will be the integral operator A = BS + I
for all t ∈ I and x ∈ R 2 . The bilinear form a can now be written as
The solution of the (21) is defined in the whole R 2 plane, in order to implement a numerical method we need first to localize to a bounded domain. A standard procedure is to find suitable function ψ asymptotically similar to the exact solution when x → ∞ (see [1, 29] ). By subtracting ψ : I × R 2 → R from u we obtain a new equation equivalent to (19) to be solved for v = u − ψ
where h ψ = h − ψ and the right hand side f is obtained by
And we obtain a variational formulation equivalent to (21) by following the steps above: To find v ∈ W 1 (I; V, H) such that
The 
We proceed as usual by truncating the unbounded range R 2 to a bounded computational domain Ω and considering Hilbert spaces H = L 2 (Ω) and V = H 1 0 (Ω) on Ω with the usual inner product (·, ·) L 2 . Since H 1 0 (Ω) is separable it contains a countable basis {ϕ k } k∈N and we restrict our selves to a finite dimensional space V N ⊂ H 1 0 (Ω) with basis {ϕ 1 , ϕ 2 , ..., ϕ N } where we consider approximated solutions of the form u(x) = N k=1 U k (t)ϕ k and an approximation of the initial condition h N = N k=1 h k ϕ k . The Galerkin formulation of (21) now reads :
for all j = 1, ..., N . There might be several possible choices for ψ, in present work we use the terminal condition (payoff). As ψ is just an approximation of the solution outside the domain of interest, a typical practice is to consider solve the problem in the outer domain Ω containing the domain of interest to keep away the error associated to the inaccuracies of the in boundary. The wider the outer domain is the less the propagation of the error inside the domain of interest. The error associated to boundary conditions can be predicted, see [30, Chapter 2.3] . The problem (28)- (29) can be written as a system of ODEs
with U (0) as in (29),
, mass matrix M, stiffness matrix A and J the matrix associated to the pure jump operator.
To solve the ODE system (30) let us consider a uniform partition {t 0 , t 2 , ... , t M } of I with norm's partition ∆t = T /M . We approximate the coefficients U (t) at the partition points, we call U k = U (t k ) to be obtained by solving the following semi-implicit θ-Scheme
Now we specify the approximation space we use along this work. Let Ω be a square domain and a family of meshes T h (Ω) defined by identical squares where h is the norm of the mesh coinciding with the length of the square's side. Let V N be the space of all continuous functions with continuous second derivatives that are piecewise third order polynomials on each element of T h (Ω) and let s 3 be the Irwin-Hall cubic spline on the pattern interval [−2, 2].
Define the mother base ϕ 0 (x) = s 3 (h −1 x 1 )s 3 (h −1 x 2 ) and for each node x (j) associated to the mesh T h we define the nodal function ϕ j as the displacement of ϕ 0
The set of all nodal basis (33) is a basis for V N . So we have set the problem in the classical finite elements framework. The first two integrals in (31) can be solved using Gaussian Legendre quadrature. The Matrix J in (31) will be approximated by the symbol method for stiffness matrices as in [1, Corollary 4.4] .
here J is the symbol of the pure jump process and using [1, Lemma 4.11] we get
for all ξ = (ξ 1 , ξ 2 ) ∈ R 2 . Integrals in (34) are suitable to be computed using Fast Fourier Transform so O (N log(N )) operations are involved in this step. Matrix J is full but is still a BTTB matrix and so are M and A in (31) as well as their linear combinations.
The general m-by-m block Toeplitz matrix with n-by-n Toeplitz blocks is defined
Matrix T mn has a blocked Toeplitz structure where each block T l , for |l| ≤ m − 1, is itself a Toeplitz matrix of order n. There exist direct methods that solve BTTB systems in O(mn log 2 (mn)) (see [31] ). Another approach is to use iterative methods, a brief description of Preconditioned Conjugated Gradient (PCG) based method is presented in [25, Chapter 5] and two different circulat preconditioners are considered to cluster the eigenvalues of T mn around 1. Of course, the use of PCG is suitable when matrix T mn is symmetric which is not our case since the action of the drift breaks-down the operator symmetry. We use the bi-conjugated gradient stabilized method (BICGSTAB) which is suitable for non-symmetric linear systems of equations. Fortunately, the circulant preconditioner c F ⊗F presented in [25] works for general BTTB matrices. Also, BTTB matrices vector product T mn b can be performed in O(mn log(mn)) operations (see [25] ). The evaluation is carried on by first embedding each block T k of matrix T mn in (35) in a circular matrix C k ∈ R n×n as follows, the block T k is defined by its first column v 1 = T k (:, 1) and its first line (with zero in the first position) v 2 = [0, T k (1, 2 : end)], and define C k by defining its first column C k (:, 1) = [v 1 , v 2 ]. In a similar way we define C 2m,2n by defining its first block column as C 2m,2n (:, 1 : 2n) = [C 0 , ..., C 1−m , 0 2n , C 1 , ..., C m−1 ]. Matrix C 2m,2n can be decomposed
where F k is the Fourier matrix of dimension k × k, C is the matrix formed by the first column of the component wise circulant matrix in which each Toeplitz matrix entry is embedded; and P is a permutation matrix that reorders columns of the identity I 4mn in the following way
being i the index of the columns of I 4mn . Matrix P reorders the result of F 2n C in a block diagonal Toeplitz matrix. After that we perform the product C 2m,2n b applying the formula (36) For a detailed explanation see [32, 25] Now we present the definition of the preconditioner c F ⊗F and a formula to obtain it. Consider a general blocked matrix
where A i,j ∈ R n×n . Given 2 unitary matrices V and U , let
where · F is the Forbenius norm. If it is the case that U = F and V = F we have the fallowing formula (see [25, eq. 5.17] )
where Q is an n-by-n circulant matrix given by
It is important to notice that c F ⊗F (T mn ) is a circulant matrix determined by its first column so there is no need to compute all entries in (39) which leads, of course, to a considerable reduction in computational time. Note also that by using iterative methods we just need to perform matrix vector products, so each iteration of the method will be done in O(nm log(nm)) and its overall performance will preserve that order as long as the preconditioner allows for just a few iterations. For a closer approach to this and others circulant preconditioners and its properties see [25] and the precedents works [33, 34, 35] . Other iterative methods can be considered, e.g. Generalized Minimal Residual (GMRES) method introduced in [36] , which is also a Krilov method; or classic non-Krilov iterative methods such as SOR, Gauss-Seidel or Jacobi. In any case, preconditioning will play an essential roll in the method's performance.
By using a Galerkin approximation in space and a θ-scheme (θ = 1 2 , Crank-Nicholson) for time evolution look for an approximate solution of problem (25) as shown in Section 3. By doing so we approximate the price of spread contracts associated to the models introduced in examples 2.1 and 2.2. In particular, we focus on cracks contracts based on the difference between the price of oil and gasoline. Specifically, we consider the prices of futures RBOB gasoline (reformulated blendstock for oxygenate blending) and West Texas Intermediate(WTI) oil in NYMEX.
The contract specifications consist on a series of strike prices around K = 1 US dollar, motivated by typical in-the-money, at-the-money and out-of-the-money contracts. Furthermore, a maturity ranging from 1 month to 1 year, initial underlying prices S (1) 0 = 100 dollars/barrel, S (2) 0 = 2 dollars/gallon with an interest rate of r = 2% are considered.
In Figure 3 the series of daily future prices of WTI in US dollar per barrel, from June 2104-July 2020 is shown, while in Figure 4 the series of future RBOB prices (in US dollar per gallon) over the same period can be observed. A rare negative price for WTI futures is observed on April 20th, 2020. We have eliminated this observation to be able to work with log-returns. It is worth noticing that the volatilities in the prices of both commodities have increased after March 2020 due to the COVID-19 pandemic. Data have been taken from https://ca.investing.com. In crack contracts c = 1 42 barrels/gallon. presence of heavy-tailed distributions. It is three times higher in WIT. As expected both assets are highly correlated. The correlation between prices is about 0.96, while the correlation between log-returns is 0.65.
Parameters in the bivariate jump-diffusion Merton model with common and idiosyncratic jumps introduced in example 2.1. Unfortunately, the number of parameters in the model is exceedingly large for most known estimation methods. Therefore, we rather arbitrarily impose values to some of the parameters, while using a Generalized 
whereȲ and S 2 (Y ) are respectively the sample mean and the sample variance of the vector Y . The time interval, measured in year units, is ∆t = 1 310 corresponding to 310 trading days. From the condition Ψ Y (i) (−i) = r and the drift b j = λ j µ (j) J + λ 0 µ 0,j we add another two constrains to the parameters:
Next, we run the scheme (32) In Table 3 we can see the performance of the algorithm at each grid of norm h = L/2 N . The first column reflects the mesh step size. The second column shows the magnitude of the error in L 2 -norm relative to the reference solution, while the third column provides the average number of iterations for time step of BICGSTAB. Finally, the fourth column shows the running time (in seconds) until a complete set of solutions C(t, x) for values of t ∈ [0, T ] in the mesh. Figure 5 shows the graphic of the reference solution U Nr with N r = 10 and Figure 6 shows the decreasing behavior of the error in norm L 2 of the approximated solution U N with respect to the reference solution U Nr through the mesh refinement.
The c (0) F F pre-conditioner performs very well since the number of iterations for time step does not grow through the mesh refinement, on the contrary, it seems to asymptotically decrease when h → 0 which is a very desirable feature. The convergence rate was 1.8106, very close to the theoretical second order of the Crank-Nicholson scheme. The computational time grows linearly with the mesh refinement. 
We set the interest rate to r = 0.02 and the loading parameters to d 1 = d 2 = 1. Also, we consider the case µ 1 = µ j to simplify calculations. Notice that the remaining parameters to be estimated come from the three subordinators, the drift parameters µ 1 and µ 2 and the volatility parameters σ 1,r and σ 2,r . Once again we implement a GMM approach letting to a least square constrained minimization problem, this time in two steps. First by matching the first four empirical and theoretical moments for the WTI asset prices taking into account risk neutral constraints, then matching the correlation to compute /mu 2 and the first two moments of the RBOB asset to compute the remaining parameters. The first matching is given by the non-linear equations
where m (j) k , j = 1, 2 k = 1, 2, 3, 4 are the empirical moments of log-returns of both commodities and m 12 is its empirical mixed moment. In addition the risk-neutral frame imposes the additional constraints
and the inequality constraint 0 < µ j + σ 2 r,j 2 < min(β 0 , β j ). It leads to a system of six variables with six equations. For convenience we re-parametrize it to new variables Same as in the previous example, we we run the scheme (32) with θ = 1/2 to obtain approximated solutions U N in meshes T h (Ω b ) of norm h = L/2 N for N = 4, 5, 6, 7, 8. Errors in norm · L 2 are computed using a reference solution U Nr where N r = 9. The BICGSTAB method was set run until find a solution with tolerance to the residual relative error of The graphic of the reference solution U Nr with N r = 10 is shown in Figure 7 while the graphic in Figure 8 shows the decreasing behavior of the error in norm L 2 of the approximated solution U N with respect to the reference solution U Nr through the mesh refinement. In Table 5 we can see the performance of the algorithm at each grid of norm h = L/2 N . The first column reflects the mesh step size. The second column shows the magnitude of the error in L 2 -norm relative to the reference solution, while the third column provides the average number of iterations for time step of BICGSTAB. The fourth column shows the running time (in seconds) until a complete set of solutions C(t, x) for values of t ∈ [0, T ] in the mesh. Again, the c
F F pre-conditioner performs very well exhibiting an average value of iterations for time step that not only doesn't grow through the mesh refinement but on the contrary, it seams to asymptotically decreases when h → 0 which is a very desirable feature. The convergence rate was 1.47, a little bit lower than the theoretical second order of the Crank-Nicholson scheme. 
The proposed strategy proved to be efficient for the valuation of spread contracts through the resolution of the associated PIDE. The symbol method allowed the efficient construction of the stiffness matrix, based on the use of the FFT, while the iterative method BICGSTAB for the solution of non-symmetric linear systems accompanied by the circular preconditioner facilitated the implementation of an implicit scheme for the temporal evolution by allowing the efficient resolution of the BTTB system associated to the Galerkin discretization. This strategy can be extended to three dimensions although we must bear in mind that the computational cost grows exponentially with the dimensionality. Other terminal conditions could be considered using this strategy as well as barrier problems.
Unfortunately, the symbol method can not be directly used with most classic finite element approximation spaces since their basis are not smooth enough to efficiently access their symbols through FFT, for example Q r spaces has polynomial interpolation basis are in H 1 0 (Ω). This issue can be addressed by substituting the basis for a mollified version of it as done in [1, Section 4.4.1] . Even if we do that, in most of the cases the basis are obtained through translations of more than one parent function, which leads us to the resolution of a linear system with a blocked matrix where each block is BTTB matrix. The system is also numerically tractable by considering an approach as in the present work but implementation becomes more complicated.
This research has been funded by NSERC and Fields Institute. where a 1 and a 2 are positive constants and the function (r) x is the maximum between zero and r. We want to know what the values of η should be for h ∈ L 1 η (R 2 ) and h ∈ L 2 η (R 2 ).
There is a function y : R → R that separates the plane in two sets, the points at which h is greater than zero and the points at which h is zero. Such a curve is as follows y(x) = log( a 1 e x + K a 2 )
So we can write where A 1 = 1 + η 1 /(1 + η 2 ), A 2 = η 1 /(1 + η 2 ), B 1 = 1 + (1 + η 1 )/η 2 , B 2 = (1 + η 1 )/η 2 C 1 = 1 + η 1 /η 2 , C 2 = η 1 /η 2
In order the integrand to have exponential decay the pairs (A 1 , A 2 ), (B 1 , B 2 ) and (C 1 , C 2 ) need to alternate signs which immediately means that η 1 has to be greater than zero because if not all signs are negative. Now, setting η 1 > 0 we have A 2 < 0, B 2 < 0 and C 2 < 0 and we need to ask for A 1 > 0, B 1 > 0 and C 1 > 0 which lead us to η 1 < −η 2 − 1. Putting all together we need to impose η ∈ (0, a − 1) × (−∞, −a) for any a > 1 in order h ∈ L 1 η (R 2 ). If we repeat the process for · L 2 η (R 2 ) we get the same therms as before (up to a constant) and three new terms from |a 2 e y − a 1 e x − K| 2 ≤ a 2 1 e 2y + a 2 2 e 2x + 2a 1 a 2 e x e y + . . . so we need the following integrals to exist using the same arguments as for L 1 η (R 2 ), we have from the first integral in (43) that η 2 < −2. We obtain another three pairs (D 1 , D 2 ), (E 1 , E 2 ) and (F 1 , F 2 ) that need to alternate signs, where D 1 = 1 + η 1 /(2 + η 2 ), D 2 = η 1 /(2 + η 2 ), E 1 = 1 + (2 + η 1 )/η 2 , E 2 = (2 + η 1 )/η 2 F 1 = 1 + (1 + η 1 )/(1 + η 2 ), F 2 = (1 + η 1 )/(1 + η 2 ).
The choice η ∈ (0, a − 2) × (−∞, −a) for any a > 2 warranties that h ∈ L 2 η (R 2 ) in fact, it also warranties that h ∈ L 1 η (R 2 ) ∩ L 2 η (R 2 ).
Let us denote by m and M the minimum and maximum eigenvalues of Σ B respectively, from (7) we have that |A(u)| = iub Q + 1 2 uΣ B u T +