Multilevel particle filters for Lévy-driven stochastic differential equations

We develop algorithms for computing expectations with respect to the laws of models associated to stochastic differential equations driven by pure Lévy processes. We consider filtering such processes as well as pricing of path dependent options. We propose a multilevel particle filter to address the computational issues involved in solving these continuum problems. We show via numerical simulations and theoretical results that under suitable assumptions regarding the discretization of the underlying driving Lévy proccess, the cost to obtain MSE O(ϵ2)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {O}(\epsilon ^2)$$\end{document} scales like O(ϵ-2)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {O}(\epsilon ^{-2})$$\end{document} for our method, as compared with the standard particle filter O(ϵ-3)\documentclass[12pt]{minimal} \usepackage{amsmath} \usepackage{wasysym} \usepackage{amsfonts} \usepackage{amssymb} \usepackage{amsbsy} \usepackage{mathrsfs} \usepackage{upgreek} \setlength{\oddsidemargin}{-69pt} \begin{document}$$\mathcal {O}(\epsilon ^{-3})$$\end{document}.


Introduction
Lévy processes have become very useful recently in several scientific disciplines.A non-exhaustive list includes physics, in the study of turbulence and quantum field theory; economics, for continuous time-series models; insurance mathematics, for computation of insurance and risk, and mathematical finance, for pricing path dependent options.Earlier application of Lévy processes in modeling financial instruments dates back in [20] where a variance gamma process is used to model market returns.
The typical computational problem in mathematical finance is the computation of the quan- For instance f can be a payoff function which may either depend on the state of Y t at a fixed time or on the whole trajectory Y .Typically one uses the Black-Scholes model, which assumes lognormal distribution for the underyling price process.However, often the asset price exhibits big jumps over the time horizon.The inconsistency of the underlying assumptions of the Black-Scholes model for market data has called for a more realistic models for these data in the literature.
General Lévy processes offer a promising alternative to describe the observed reality of financial market data, as compared to models that are based on standard Brownian motions.
In the application of standard and multilevel particle filter methods to SDEs driven by general Lévy processes, in addition to pricing path dependent options, we will consider filtering of partially observed Lévy process with discrete-time observations.In the latter context, we will assume that the partially observed data are regularly spaced observations z 1 , . . ., z n , where z k ∈ R d is a realization of Z k and Z k |Y kτ has density given by g (z k |Y kτ ), where τ is the time scale.Real S&P 500 stock price data will be used to illustrate with our proposed methods as well as the standard particle filter.We will show how both of these problems can be formulated as general Feynman-Kac type problems [5], with time-dependent potential functions modifying the Lévy path measure.
Recently, [6] provided a detailed analysis of the multilevel Monte Carlo (MLMC) methodology first introduced by [9] for a Lévy-driven stochastic differential equation.The authors in [8] us the MLMC method for general Lévy processes based on Wiener-Hopf decomposition.We extend the methodology described in [6] to particle filtering framework.This is more challenging due to the following reasons.First, one must choose suitable weighting function to prevent the weights in the particle filter being zero (or infinite).Next, one must control the jump part of the underlying Lévy process such that the path of the filter does not blown up as the time parameter increases.
In pricing path dependent options, for example knock out barrier options, we adopt the same strategy described in [14,15] for the computation of the expectation of the functionals of the SDE driven by general Lévy processes.
The rest of the paper is organised as follows.In Section 2, we briefly review the construction of general Lévy processes with application to multilevel settings.The numerical approximation of the Lévy-driven SDE is presented, and the vanilla MLMC method framework explained.We conclude the section by describing the construction of a coupled kernel for Lévy-driven SDEs.Section 3 introduces both the standard and multilevel particle filter methods and its application to the Lévy-driven SDEs.Section 4 features numerical examples of pricing barrier options and filtering of partially observed Lévy processes.The computational savings of the multilevel particle filter over the standard particle filter is illustrated in this section.

Approximating SDE driven by Lévy Processes
In this section, we briefly describe the construction and approximation of a general d -dimensional Lévy process {X t } t∈[0,K] , and the solution Y := {Y t } t∈[0,K] of a d-dimensional SDE driven by X.Consider a stochastic differential equation given by where a : R d → R d×d satisfies the usual growth assumptions on solutions of SDEs, and the initial value is y 0 (assumed known).In particular, we are interested in computing the expectation of bounded functions f :

Lévy Processes
For a general detailed description of the Lévy processes, we shall refer the reader to the monographs of [3,23] and analysis of SDEs driven by Lévy process [1,21].Lévy processes are stochastic processes with stationary and independent increments, which begin almost surely from the ori- There is a strong interplay between Lévy processes and infinitely divisible distributions such that, for any t > 0 the distribution of X t is infinitely divisible.Conversely, if F is an infinitely divisible distribution then there exists a Lévy process X t such that the distribution of X 1 is given by F .This conclusion is the result of Lévy-Khintchine formula for Lévy processes we describe below.Let X be a Lévy process with a triplet (ν, Σ, b), b ∈ R d , 0 < Σ = Σ T ∈ R d×d , and a measure ν, satisfying ν({0}) = 0 and with π the probability law of X t , where The measure ν is called the Lévy measure of X t .The triplet of Lévy characteristics (ν, Σ, b) is simply called Lévy triplet.Note that in general, the Lévy measure ν can be finite or infinite.
If ν(R) < ∞, then almost all paths of the Lévy process have a finite number of jumps on every compact interval and it can be represented as a compensated compound Poisson process.On the other hand, if ν(R) = ∞, then the process has an infinite number of jumps on every compact interval almost surely.Even in this case the third term in the integrand ensures that the integral is finite, and hence so is the characteristic exponent.

Simulation of Lévy Processes
The law of increments of many Lévy processes is not known explicitly.This makes it more difficult to simulate a path of general Lévy processes than for instance standard Brownian motion.
For a few Lévy processes where the distribution of the process is known explicitly, [4,24] provided methods for exact simulation of paths of such processes, which are applicable in financial modelling.For our purposes, the simulation of the path of a general Lévy process will be based on the Lévy-Itô decomposition and we briefly describe the construction below.An alternative construction is based on Wiener-Hopf decomposition.This is used in [8].
The Lévy-Itô decomposition reveals much about the structure of the paths of a Lévy process.
We can split the Lévy exponent, or the characteristic exponent of X t in (2), into three parts where The first term corresponds to a deterministic drift process with parameter b, the second term to a Wiener process with covariance √ Σ, and the last part corresponds to a Lévy process which is a square integrable martingale.This term may either be a compensated compound Poisson process or the limit of such processes, and is the hardest to handle when it arises from such limit.
Thus, any Lévy process can be decomposed into three independent Lévy processes thanks to the Lévy-Itô decomposition theorem.In particular, let {W t } t∈[0,K] denote a Wiener process independent of the process {L t } t∈[0,K] .We define the Lévy process {X t } t∈[0,K] by The Lévy-Itô decomposition guarantees that every square integrable Lévy process has a representation as (3).We will assume that one cannot sample from the law of X t , hence of Y t , without discretization error.It is assume that one will discretize the paths of the Lévy process (e.g.Euler or Milstein [18]) and call the time discretization h l .Analysis of convergence of such numerical methods have been studied extensively in [13,22].
Item (i) provides continuity of the forward map, while (ii) controls the variance of the jumps, and (iii) controls the diffusion and drift components.These assumptions are the same as in the paper [6], and as in that paper we refer to the following general references on Lévy processes for further details [1,3].

Numerical Approximation of a Lévy Process and Lévy-driven SDE
Recall (1) and (3).For pedagogical purposes, assume that K = 1 and consider the evolution of discretized Lévy process and hence the Lévy-driven SDE over the time interval [0, 1].
In order to describe the Euler discretization of the two processes for a given accuracy parameter h l , we need some definitions.Let δ l denote the user specified jump threshold parameter over the time interval [0, 1].Then for δ l > 0, let that is the Lévy measure outside the ball of radius δ l .The parameter δ l will be chosen such that the step-size of the time-stepping method is h l = 1/λ l .The jump time increments are exponentially distributed with parameter λ l so that the number of jumps before time t is a Poisson process N l (t) with intensity λ l .The jump times will be denoted by T l j .The jump heights ∆L l Tj are distributed according to Define The expected number of jumps on an interval of length t is F l 0 t, and the compensated compound Poisson process is an L 2 martingale which converges in L 2 to the Lévy measure L t [1].
The Euler discretization of the Lévy process, and hence the Lévy driven SDE, is given by Algorithm 1.Note that, in the construction of the discretized process, appropriate refinement of the original jump times { T l j } to new jump times {T l j } is necessary to control the discretization error arising from the Brownian motion component, the original drift process, and the drift component of the compound Poisson process.Note that the ∆L l T l j is non-zero only when T l j corresponds with T l m for some m.
The numerical approximation of the Lévy process described in Algorithm 1 gives rise to an approximation of the Lévy-driven SDE as follows.Given where (∆X) l is given by (4).In particular the recursion in (5) gives rise to a transition kernel denoted Q l (u, y) between observation times t ∈ {0, 1, . . ., K}.This kernel Observe that the initial condition for
(B) Generate jump heights: (C) Refinement of original jump times: If T l j = 1, k l = j; Go to (D); Otherwise j = j + 1; Go to start of (C).
(D) Recursion of the process: X is irrelevant for simulation of Y , since only the increments (∆X) l T l m+1 are required, which are simulated independently by adding a realization of Remark 2.1.The numerical approximation of the Lévy process and hence Lévy-driven SDE (1) in Algorithm 1 is the single-level version of a more general coupled discretization [6] which will be described shortly in Section 2.5.This procedure will be used to obtain samples for the plain particle filter algorithm.

Multilevel Monte Carlo Method
Suppose one aims to approximate the expectation of functionals of the solution of the Lévy- the expectation with respect to the density associated with the Euler discretization (4) at level L. The standard Monte Carlo (MC) approximation at time 1 consists in obtaining i.i.d.samples from the density π L 1 and approximating π L 1 (f ) by its empirical average ).
The mean square error of the estimator is Since the MC estimator π L,N L 1 (f ) is an unbiased estimator for π L 1 (f ), this can further be decom- The first term in the right hand side of the decomposition is the variance of MC simulation and the second term is the bias of approximation induced by the approximate Euler discretization with a mixture of random and deterministic grid size.If we want (6) to be O( 2), then it is clearly necessary to choose N ∝ −2 , and then the total cost is ) ∝ −γ for some γ > 0 is the cost to ensure the bias is O( ).Now, in the multilevel Monte Carlo (MLMC) settings, one can observe that the expectation of the finest approximation π L 1 (f ) can be written as a telescopic sum starting from a coarser approximation π 0 1 (f ), and the intermediate ones: The idea of the MLMC method is to approximate the multilevel (ML) identity ( 7) by independently computing each of the expectations of the telescopic sum by a standard MC method.This is possible by obtaining for each l, from a joint measure πl 1 with the appropriate marginals π l 1 and π l−1 1 , for example generated from a coupled simulation the Euler discretization of SDE (1) at successive refinements.The construction of such a coupled kernel is detailed in Section 2.5.Suppose it is possible to obtained such coupled samples at time 1, then in particular, for l = 0, . . ., L, one has independent MC estimates.Let ) , where Analogously to the single level Monte Carlo method, the mean square error for the multilevel estimator ( 8) can be expanded to obtain with the convention that f (Y −1 1 ) ≡ 0. It is observed that the bias term remains the same; that is we have not introduced any additional bias.However, by an optimal choice of N 0:L , one can possibly reduce the computational cost for any pre-selected tolerance of the variance of the estimator, or conversely reduce the variance of the estimator for a given computational effort.
In particular, for a given user specified error tolerance measured in the root mean square error, the highest level L and the replication numbers N 0:L are derived as follows.We make the following assumptions about the bias, variance and computational cost based on the observation that there is an exponential decay of bias and variance as L increases.
Suppose that there exist some constants α, β, γ and an accuracy parameter h l associated with the discretization of SDE (1) at level l such that where α, β, γ are related to the particular choice of the discretization method and cost is the computational effort to obtain one sample Y l , Y l−1 .For example, the Euler-Maruyama discretization method for solution of SDEs driven by Brownian motion gives orders α = β = γ = 1.
The accuracy parameter h l typically takes the form h l = S −l 0 for some integer S 0 ∈ N.
The key observation from the mean-square error of the multilevel estimator (8) − (9) is that the bias is given by the finest level, while the variance is decomposed into a sum of variances of the l th increments.Thus the total variance is of the form and by condition (V l ) above, the variance of the l th increment is of the form V l N −1 l .The total computational cost takes the form C = L l=0 C l N l .In order to minimize the effort to obtain a given mean square error (MSE), one must balance the terms in (9).Based on the condition (B l ) above, a bias error proportional to will require the highest level In order to obtain optimal allocation of resources N 0:L , one needs to solve a constrained optimization problem: minimize the total cost C = L l=0 C l N l for a given fixed total variance or vice versa.Based on the conditions (V l ) and (C l ) above, one obtains via the Lagrange multiplier method the optimal allocation Now targetting an error of size O( ), one sets N l ∝ −2 h (β+γ)/2 l K( ), where K( ) is chosen to control the total error for increasing L. Thus, for the multilevel estimator we obtained: One then sets K( ) = (ii).If β > γ, which correspond to the Milstein scheme, then K( ) ≡ 1, and hence the optimal computational cost is O( −2 ).
(iii).If β < γ, which is the worst case scenario, then it is sufficient to choose K( . In this scenario, one can easily deduce that the total cost is O( −(γ/α+κ) ), where One of the defining features of the multilevel method is that the realizations (Y l 1 , Y l−1

1
) for a given increment must be sufficiently coupled in order to obtain decaying variances (V l ).This means once the increments of the process have been sampled to obtain a draw from π l 1 (f ), a deterministic transformation of the sample points and of the increments should take place to obtain a sample from π l−1 1 (f ).This was particularly clear in the context of stochastic differential equations driven by Brownian motion introduced in [9] (see also [16]), where coarse paths are obtained by summing the fine paths, but it is non-trivial how to proceed in the context of SDEs purely driven by general Lévy processes.Both and [8] suggested A Poisson thinning technique based on Wiener-Hopf factorization has been suggested by [10] for pure-jump diffusion and by [8] for general Lévy processes.In the next section, we explain an alternative construction of a coupled kernel based on the Lévy-Ito decomposition, in the same spirit as in [6].

Coupled Sampling for Levy-driven SDEs
The ML methodology described in Section 2.4 works by obtaining samples from some coupledkernel associated with discretization of (1).We now describe how one can construct such a kernel associated with the discretization of the Lévy-driven SDE.Let u = (y, y ) ∈ R 2d .Define a kernel, , where σ(.) denotes the σ-algebra of measurable subsets, such that The coupled kernel M l can be constructed using the following strategy.Using the same definitions in Section 2.3, let δ l and δ l−1 be user specified jump-thresholds for the fine and coarse approximation, respectively.Define The objective is to generate a coupled pair (Y l,l 1 , Y l,l−1 The parameter δ (h ) will be chosen such that h −1 = ν(B c δ ), and these determine the value of F 0 in ( 13), for ∈ {l, l − 1}.We now describe the construction of the coupled kernel M l and thus obtain the coupled pair in Algorithm 2, which is the same as the one presented in [6].(2) Generate coarse jump times and heights: for j l ∈ {1, . . ., k l l } , and T l,l−1 j l−1 = T l,l j l ; j l−1 = j l−1 + 1; (3) Refine jump times: Set j l−1 = j l = 1 and T l,l−1 If T l,l j l = 1, set k l l = j l , and redefine T l,l i := T l,l i for i = 1, . . ., k l l ; Else j l = j l + 1 and Go to (ii).
(4) Recursion of the process: sample where ∆W T l, The construction of the coupled kernel M l outlined in Algorithm 2 ensures that the paths of finer and coarse processes are correlated enough to ensure that the optimal convergence rates of the multilevel algorithm is achieved.

Multilevel Particle Filter for Lévy-driven SDEs
In this section, the multilevel particle filter will be discussed for sampling from certain types of measures which have a density with respect to a Lévy process.We will begin by briefly reviewing the general framework and standard particle filter, and then we will extend these ideas into the multilevel particle filtering framework.

SDEs
Recall the Lévy-driven SDE (1).It will be assumed that the general probability density of interest is of the form for n ≥ 1, for some given y 0 where Q ∞ (y (i−1) , y) is the transition density of the process (1) as a function of y, i.e. the density of solution Y 1 at observational time point 1 given initial condition Y 0 = y (i−1) .It is assumed that is the conditional density (given y i ) of an observation at discrete time i, so observations (which are omitted from our notations) are regularly observed at times 1, 2, . . . .Note that the formulation discussed here, that is for η∞ n , also allows one to consider general Feynman-Kac models (of the form ( 16)), rather than just the filters that are focussed upon in this section.The following assumptions will be made on the likelihood functions {G i }.Note these assumptions are needed for our later mathematical results and do not preclude the application of the algorithm to be described.
Assumption 3.1.There are c > 1 and C > 0, such that for all n > 0, and In practice, as discussed earlier on Q ∞ is typically analytically intractable (and we further suppose is not currently known up-to a non-negative unbiased estimate).As a result, we will focus upon targets associated to a discretization, i.e. of the type for l < ∞, where Q l is defined by k l iterates of the recursion in (5).Note that we will use ηl n as the notation for measure and density, with the use clear from the context, where l = 0, 1, . . ., ∞.
The objective is to compute the expectation of functionals with respect to this measure, particularly at the last co-ordinate.For any bounded and measurable function f : We will also use the notation ηl Often of interest is the computation of the un-normalized measure.That is, for any bounded and measurable function f : In the context of the model under study, ζl n (1) is the marginal likelihood.

Particle Filtering
We will describe the particle filter that is capable of exactly approximating, that is as the Monte Carlo samples go to infinity, terms of the form ( 18) and ( 19), for any fixed l.The particle filter has been studied and used extensively (see for example [5,7]) in many practical applications of interest.
For a given level l, algorithm 3 gives the standard particle filter.The weights, are defined as with the convention that w l,(i) 0 = 1.Note that the abbreviation ESS stands for effective sample size which measures the variability of weights at time k of the algorithm (other more efficient procedures are also possible, but not considered).In the analysis to follow H = 1 in algorithm 3 (or rather it's extension in the next section), but this is not the case in our numerical implementations.
[5] (along with many other authors) have shown that at step 3 of algorithm 3, the estimate will converge almost surely to (18).In addition, if H = 1 in algorithm 3, will converge almost surely to (19).
Algorithm 3 : Particle filter Iterate 2 and 3 If ESS < H (for some threshold H), resample the particles {Y l,(i) k } N l i=1 and set all weights to w l,(i) k } N l i=1 by using (20).

Multilevel Particle Filter
We now describe the multilevel particle filter of [16] for the context considered here.The basic idea is to run L+1 independent algorithms, the first a particle filter as in the previous section and the remaining, coupled particle filters.The particle filter will sequentially (in time) approximate η0 k and the coupled filters will sequentially approximate the couples (η 0 k , η1 k ), . . ., (η L−1 k , ηL k ).Each (coupled) particle filter will be run with N l particles.
The MLPF is as follows.In the below description, we set H = 1 (as in algorithm 3), but it need not be the case.Recall that the case l = 0 is just a particle filter.For each 1 ≤ l ≤ L the following procedure is run independently.
Hence Theorem D.1 of [16] follows as well with κ = 2, and hence Theorem 4.2 of [16].The additional assumption that γ/α ≤ 2 ensures Corollary D.1 of [16] holds, and the result follows as in the proof of Theorem 4.1 of [16] 4 Numerical Examples In this section, we compare our proposed multilevel particle fitler method with the vanilla particle filter method.A target accuracy parameter will be specified and the cost to achieve an error below this target accuracy will be estimated.The performance of the two algorithms will be compared in two applications of SDEs driven by general Lévy process: filtering of a partially observed Lévy process (S&P 500 stock price data) and pricing of a path dependent option.In each of these two applications, we let X = {X t } t∈[0,K] denote a symmetric stable Lévy process, i.e.X is a (ν, 0, 0)-Lévy process with Lebesgue density of the Lévy measure given by with c > 0, x * > 0 (the truncation threshold) and index φ ∈ (0, 2).The parameters c and x * are both 1 for all the examples considered.The Lévy-driven SDE considered here has the form with y 0 assumed known, and a satisfies Assumption 2.1(i).Notice Assumption 2.1(ii-iii) are also satisfied by the Lévy process defined above.In the examples illustrated below, we take a(Y t ) = Y t , y 0 = 1, and φ = 0.5.
Remark 4.1 (Symmetric Stable Lévy process of index φ ∈ (0, 2)).In approximating the Lévydriven SDE (23), [6] provided in Theorem 2 of the paper, asymptotic error bounds for the strong approximation by the Euler scheme.If the driving Lévy process X t has no Brownian component, that is W = 0, then the L 2 -error rate is bounded by and for W = 0, for a fixed constant C < ∞ (that is the Lipschitz constant), where σ 2 (δ l ) = B c δ l |x| 2 ν(dx) (see also the appendix).Recall that δ l (h l ) is chosen such that h l = 1/ν(B c δ l ).One obtain the analytical expression for some constant C > 0. One can also analytically compute Now, setting h l = 2 −l , one obtains so that the Lévy measure ν(B c δ l ) = 2 l .Then, one can easily bound (25) by ).Using (24)-( 25) and the error bounds for W = 0, one can straightforwardly obtain strong error rates for the approximation of SDE driven by stable Lévy process in terms of the single accuracy parameter h l .This is given by Thus, if b − F l 0 = 0, the strong error rate β associated with a particular discretization is given by Otherwise it is just given by (2 − φ)/φ.

Partially observed data
In this section we consider filtering a partially observed Lévy process.Recall that the Lévydriven SDE takes the form (23) We shall take the test function f (y) = e y for the example considered below.q q q q q q q q q q q q q q 10 2 10 3 10 4

Barrier Option
Here we consider computing the value of a discretley monitored knock out barrier option (see e.g.[11] and the references therein).Let Y 0 ∈ [a, b] for some 0 < a < b < +∞ known and let Q ∞ (y (i−1) , y) be the transition density of the process as in (23).Then the value of the barrier option (up-to a known constant) is for K > 0 given.As seen in [14] the calculation of the barrier option is non-trivial, in the sense that even importance sampling may not work well.We consider the (time) discretized version Define a sequence of probability densities, k ∈ {1, . . ., n} for some non-negative collection of functions Gk (y k ), k ∈ {1, . . ., n} to be specified.Then the value of the time discretized barrier option is exactly where f (y n ) = max{y n −K, 0}.Thus, we can apply the MLPF targetting the sequence {η l k } k∈{1,...,n},l∈{0,...,L} and use our normalizing constant estimator (21).The fitted linear model of log MSE against log cost has a slope of −0.6667 and −0.859 for PF and MLPF respectively.These numerical results are consistent with the expected theoretical asymptotic behaviour of MSE∝Cost −1 for the multilevel method.The single level particle filter achieves the asymptotic behaviour of the standard Monte Carlo method with MSE∝Cost −2/3 .
The following assumption will be made, uniformly over the level l ∈ [0, 1, . ..), which will be omitted for notational simplicity.The following assumptions will be required for the MLPF, for all suitable test functions ϕ ∈ A. where COST[M l ] is the cost to simulate one sample from the kernel M l , which is a random variable since the jump times are random, and E denotes its expectation.
Proposition A.1.Assume 2.1 holds.Let Y l t and Y l t denote the solutions of the discretized process at time t, given initial conditions y and y , respectively.Then there is a C > 0 such that Proof.Denote by T = {T j } the set of all discrete time-steps, which may or may not correspond to a jump time.Letting T j (t) = max{T j ∈ T; T j < t}, and letting ∆ t X = X l t − X l Tj (t) we have tity E [f (Y )], where Y := {Y t } t∈[0,K] solves a stochastic differential equation driven by Lévy 1 arXiv:1804.04444v1[stat.CO] 12 Apr 2018 processes and f ∈ B b (R d × [0, K]), a bounded Borel measurable function on R d × [0, K].
gin and are stochastically continuous.Two important fundamental tools available to study the richness of the class of Lévy processes are the Lévy-Khintchine formula and the Lévy-Itô decomposition.They respectively characterize the distributional properties and the structure of sample paths of the Lévy process.Important examples of Lévy process include Poisson processes, compound Poisson processes and Brownian motions.

Assumption 2 . 1 .
There exists a C > 0 such that (i) |a(y) − a(y )| ≤ C|y − y |, and |a(y)| ≤ C for all y ∈ R ; in order to have variance of O( 2 ).We can identify three distinct cases (i).If β = γ, which corresponds to the Euler-Maruyama scheme, then K( ) = L.One can clearly see from the expression in (10) that L = O(| log( )|).Then the total cost is O( −2 log( ) 2 ) compared with single level O( −3 ).

Figure 2 :Figure 2 .
Figure 2: Mean square error against computational cost for filtering with partially observed data.

For
this example we choose K = 1.25, a = 0, b = 5, y 0 = 1, n = 100.Gk (y k ) = |y k − K| κ k , where κ k is an annealing parameter with κ k → 1 as k → n.The N l are chosen as in the previous example.The error-versus-cost plots for PF and MLPF are shown in Figure 3.Note that the bullets in the graph correspond to different choices of L (for both PF and MLPF, 2 ≤ L ≤ 8).
where f : R d → R is Borel-measurable functions on R d .Typically, one is interested in the expectation w.r.t. the law of exact solution of SDE (1), but this is not always possible in practice.Suppose that the law associated with (1) with no discretization is π 1 , and since we cannot sample from π 1 , we used a biased version π L 1 associated with a given level of discretization of SDE (1) at time 1 . In addition, partial observations {z 1 , . . ., z n } are available with Z k obtained at time k and Z k |Y k has a density function G k (y k ) (the observation is omitted from the notation).The observation density is Gaussian with mean y k and variance 1.We aim to estimate E[f (Y k )|z 1:k ] for some test function f (y).In this application, we consider the real daily S&P 500 log return data (from August to July 24, 2015, normalized to unity variance).