key: cord-0597925-7cuocele
authors: Olivares, Pablo; Diaz, Ciro
title: A Finite Elements Strategy for Spread Contract Valuation Via Associated PIDE
date: 2020-09-17
journal: nan
DOI: nan
sha: 531912df7dec894ab5f4f14ae1a8868e42a8bb9c
doc_id: 597925
cord_uid: 7cuocele

We study an efficient strategy based on finite elements to value spread options on commodities whose underlying assets follow a dynamic described by a certain class of two-dimensional Levy models by solving their associated partial integro-differential equation (PIDE). To this end we consider a Galerkin approximation in space along with an implicit scheme for time evolution. Diffusion and drift in the associated operator are discretized using an exact Gaussian quadrature, while the integral part corresponding to jumps is approximated using the symbol method recently introduced in the literature. A system with blocked Toeplitz with Toeplitz blocks (BTTB) matrix is efficiently solved via biconjugate stabilized gradient method (BICSTAB) with a circulant pre-conditioner at each time step. The technique is applied to the pricing of textit{crack} spread options between the prices of futures RBOB gasoline (reformulated blendstock for oxygenate blending) and West Texas Intermediate(WTI) oil in NYMEX.

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 +

A flexible galerkin scheme for option pricing in lévy models

Basket option pricing approximations under jump-diffusion model

Multivariate time changes for lévy asset models: Characterization and calibration

Basket option pricing and implied correlation in a one-factor lévy model

Stepping through fourier space

A fast and accurate fft-based method for pricing early-exercise options under lévy processes

Pricing and hedging in exponential lévy models: review of recent results

Viscosity solutions of nonlinear integro-differential equations

Backward stochastic differential equations and integralpartial differential equations

Optimal stopping of controlled jump diffusi on processes: a viscosity solution approach

On kolmogorov equations for anisotropic multivariate lévy processes

Fast deterministic pricing of options on lévy driven assets

Variational solutions of the pricing pides for european options in lévy models

A feynman-kac-type formula for lévy processes with discontinuous killing rates

Wavelet discretizations of parabolic integrodifferential equations

Wavelet galerkin pricing of american options on lévy driven assets

Computational methods for quantitative finance

Superfast solution of real positive definite toeplitz systems

Asymptotically fast solution of toeplitz and related systems of linear equations

A new algorithm for solving toeplitz systems of equations

A generalization of the levinson algorithm for hermitian toeplitz matrices with any rank profile

A look-ahead bareiss algorithm for general toeplitz matrices

A proposal for toeplitz matrix calculations

Linear and Nonlinear Deconvolution Problems (Optimization)

An introduction to iterative Toeplitz solvers

Bi-cgstab: A fast and smoothly converging variant of bi-cg for the solution of nonsymmetric linear systems

Numerical methods for lévy processes

Classification of lévy processes with parabolic kolmogorov backward equations

A finite element framework for option pricing with the bates model

Linear and non-linear monotone methods for valuing financial options under two-factor, jump-diffusion models

A fast algorithm for solving a toeplitz system of equations

Toeplitz and circulant matrices: A review

An optimal circulant preconditioner for toeplitz systems

The circulant operator in the banach algebra of matrices

Optimal and superoptimal circulant preconditioners

Gmres: A generalized minimal residual algorithm for solving nonsymmetric linear systems