{'fraction_non_alphanumeric': 0.04096745297103613, 'fraction_numerical': 0.0075246342191699015, 'mean_word_length': 4.419417475728156, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 4, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 2, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
On a numerical method for solving integro-differential equations with variable coefficients with applications in finance

We propose a new general numerical method aimed to solve integro-differential equations with variable coefficients. The problem under consideration arises in finance where in the context of pricing barrier options in a wide class of stochastic volatility models with jumps. To handle the effect of the correlation between the price and the variance, we use a suitable substitution for processes. Then we construct a Markov-chain approximation for the variation process on small time intervals and apply a maturity randomization technique. The result is a system of boundary problems for integro-differential equations with constant coefficients on the line in each vertex of the chain. We solve the arising problems using a numerical Wiener-Hopf factorization method. The approximate formulae for the factors are efficiently implemented by means of the Fast Fourier Transform. Finally, we use a recurrent procedure that moves backwards in time on the variance tree. We demonstrate the convergence of the method using Monte-Carlo simulations and compare our results with the results obtained by the Wiener-Hopf method with closed-form expressions of the factors.


Introduction
In the years since Dynkin, Feinman, and Katz published their fundamental works, it has been well known, that the solution to a diffusion equation in a given domain can be interpreted as an expected value of a function of Brownian motion at the first exit point from the domain. This idea was later generalized for the processes combining both diffusion and Lévy models, trajectories of which can be discontinuous. Nowadays Lévy processes are widely used for modeling various social end economical processes. Both researchers and financial institutes pay a lot of attention to the processes. From the point of view of financial mathematics, development of means for solving these equations allows to solve an essential applied problem -the problem of estimation of risk-neutral prices of derivative financial instruments in models with jumps, in particularbarrier options.
Thus, calculating functionals from Markov processes with jumps plays an important role in quantitative finance due to the fact that the presence of jumps even on liquid assets became more pronounced in recent decades. This task, nonetheless, is mathematically complex and comes down to solving complex partial integro-differential equations known as forward and backward Kolmogorov equations with certain initial and boundary conditions. As a rule, analytical methods are useless for solving such equations, which brings a necessity to use cutting-edge methods of numerical analysis.
Let's consider a system of stochastic differential equations in Ito form that describes a dynamics of a position S t of a particle on a line, that is connected with variance V t driven by CIR process [1].
where r is a non-negative parameter, Z S t and Z V t are Wiener processes with correlation coefficient ρ. Jumps are defined by Poisson process N t with intensity λ J , that counts a quantity of identically distributed jumps of size J to the moment t. The process N t does not depend on processes Z S t and Z V t , and is independent of J as well. Parameter κ V defines the speed of "regression" of variance process to a "long-run" variance average value θ V , σ V > 0 is referred to as a "volatility" of variance. The value of ζ is fixed by the distribution parameters of the jump size J so that exp(−rt)S t is a martingale.
In financial applications S t in (1) simulates dynamics of an underlying asset, whereas a parameter r ≥ 0 denotes a risk-free interest rate. If jumps J are distributed log-normally, J ∼ LogN(µ J , σ 2 J ), then ζ = e µ J + 1 2 σ 2 J − 1 and we obtain a well-known "Bates model" initially presented in paper [2]. By setting jump intensity to zero while keeping a variance process dynamics, we will obtain a system that represents a Heston model -one of the most widelyused models with stochastic volatility [3].
Replacing the process V t with a positive constant V 0 will make the first equation of the system (1) describe a Lévy process which is referred to as "jump-diffusion" (see, for example, [4]). Particularly, if J is log-normally distributed, we get Merton model [5]. With the intensity of jumps set to zero the model finally comes down to the most deeply explored models in financethe Black-Scholes model [6].
The goal of the research is to develop a numerical method to calculate the following conditional expectation for the class of models described by (1), where H is an absorbing barrier, S T (= inf 0≤t≤T S T ) is infimum process for the process S t . Functionals like (2) arise in financial mathematics when solving a problem of pricing barrier options in models described by the system (1). A "barrier option" is a contract that pays a certain sum G(S T ) on expiration time T , if and only if the underlying asset price S T does not cross a certain barrier H from above (down-and-out barrier option) or from below H (up-andout barrier option). For example, for an option that grant a holder a right to sell a certain asset for a certain price K ("put option"), G(S) = max{0, K − S}, whereas for the option that grants a right to buy a certain asset for a certain price K ("call option"), G(S) = max{0, S − K}.
Suppose τ = T − t, then F (S, V, τ ) (= f (S, V, T − τ )) satisfies the following (see, for example, where f (J) is a probability density function for jump amplitudes J. Because H is an absorbing barrier and S t has an ability to cross the barrier, we need to include condition F (S, V, τ ) = 0, S(τ ) ≤ H.
Most common approaches to solve such equations are the Monte-Carlo method [8], finitedifference schemes (see, for example, [9]) and a hybrid approach that combines tree method and finite-difference schemes [10].
The most significant drawback of Monte-Carlo methods is low computational performance, because for the models with jumps and problems with barriers detailed trajectory simulation becomes inevitable. Finite-difference schemes approach highly depends on exact kind of Lévy model and does not provide much accuracy because a corresponding matrix of a system is dense (unlike diffusion models), even though many authors implicitly use only its tridiagonal part for its inversion. Some authors combine split method and finite-difference schemes while approximating an integral component as a function of a symbol of the corresponding pseudo-differential operator from finite differences [11], but for popular Lévy models it leads to computing matrix logarithms and exponents, which can lead to sufficient computational errors.
For the sake of dimension reduction some authors (see, for example, [12,13]) suggest using Markov chain with continuous time for volatility process approximation. The simplest case of this approach is a binomial model [13]. As a result, we obtain Lévy model with regime-switching volatility. Markov chain is constructed in a way to save corresponding drift and volatility for transitional probabilities from every state.
In a series of articles (see, for example, [14]) a general method of approximate Wiener-Hopf factorization for a wide class of Lévy processes with a constant variation was used to solve problems with barriers for integro-differential equations, to which evaluation of prices for different options comes down. The numerical implementation of Wiener-Hopf operators was performed with the help of Fast Fourier Transform. This made the computational methods comparable to finite-difference schemes in simplicity of implementation, yet significantly faster and more accurate, as numerical experiments show [14].
One of the core ideas of a new numerical method proposed for solving a problem (2) is to use a new enhanced approximate Wiener-Hopf factorization, initially presented in [15]. The new formulae are more general-purposed in comparison to [14] and easier to use, as they do not necessarily require the separation of main part of a characteristic function under factorization. In order to solve a problem after time discretization (for example, Carr randomization) the coefficients of the integro-differential operator obtained are going to be approximated using a tree method which models a volatility process, leading therefore a problem in each state to a model with a constant volatility, and new explicit formulae for Wiener-Hopf factors are going to be used for a respective task. As a result, we obtain a method comparable in speed to a finite-difference scheme, but more general and accurate.

. Process substitution
One important feature of the model considered is an effect of correlation between Wiener processes in the equations for price and variation. To solve these equations one must take it into account while developing a numerical procedure. The method we use is similar to the one in [10], and aimed at constructing a suitable substitution for processes. Let us definê where W t Z t are independent Brownian motions. Let us first also introduce a substitution for the process S t , setting In these terms, Being constructed in such way, the substitution allows to proceed a logarithmic scale and make a normalization with respect to a barrier at the same time, which makes computation more convenient. After making both substitutions, the system (1) will be as follows: Let's introduce a simpler notation for the drifts we use later:

Carr randomization
In order to be able to operate on small timeframes we can use a procedure known as "Carr randomization" ("Canadization"), first introduced in [19] and generalized for a more general case of stochastic control problems in paper [17]. Let's denote a Carr approximation for the function f (S, V, t) at the moment t n , where t i = i∆τ and ∆τ = T N .
-be a sequence of independent exponentially distributed variables with an average value of ∆τ . Let's denote Z n τ = Y tn+τ + ρ σ V V tn+τ . Assuming f N (y, v) = G(H exp(y + ρ σ V v)), we have an opportunity to write down the following: This form of expectation can be iteratively calculated for the case of Lévy model with fixed variance and comes down to a sequence of problems. However, the presence of stochastic variance in our case makes it necessary to use an approximation for the CIR process V , in order to compute it.

Approximation
The next goal we would like to achieve is to reduce the initial problem to a set of problems with lower dimension and obtain a Lévy model with a regime-switching volatility. For this purpose we construct approximation of a CIR process using a Markov chain following the procedure from [13]. Among the advantages of the scheme we can mention fast convergence towards the parent process and the ability to remain adequate in a wide range of parameters. In particular, keeping the Feller condition 2κ V θ V > σ 2 V r is not significant in this case for the procedure to remain correct.
The idea of the approximation is that in every moment of time t n a variation can take one of the states V (n, k). In the moment t n+1 from the vertex (n, k) we can either ho "up" -to a vertex (n + 1, k u ), or "down" -to (n + 1, k d ). The vertices k u and k d are being calculated in accordance with the drift µ(V (n,k) ), using the following rules: Let us define transitional probabilities as follows: .
In order to maintain correct work of the scheme for different parameters values we should introduce additional rules that make impossible for the negative probabilities to appear in some vertices:

Approximate factorization
The tree procedure described above can now be used to obtain problems with a fixed variation. With regard to the Markov chain, the expectation (3) for n = N − 1, .., 0 can be represented in the form: Thus, for each of vertices (n, k), n = N − 1, ...0 two expectations to compute arise: one in the assumption that transition has gone up to a vertex (n + 1, k d ), another one in the assumption that transition has gone down to a vertex (n + 1, k u ).
Let us introduce a shorthand notation f k n (y) = E tn [e −rτn 1 Z n τn >0 f n+1 (Y tn+τn , V (n + 1, k))|Y tn = y]. It is well known (see e.g. [14]) that the expectation f k n (y) is a solution to the following integro-differential equation: where q = ∆t −1 + r and The integro-differential operator L n,k can be interpreted as a pseudo-differential operator with the symbol ψ n,k (ξ) -a characteristic exponent of a process Y t with V t = V (n, k): where φ(ξ) is a characteristic exponent of a compound Poisson process. For example, for a Merton model it looks as follows: To find each of f k n (y), one should solve two problems: where q = ∆t −1 + r; n = 0, 1, ..., N ; k = 0, 1, ..., n. Solutions to each of the problems (4)-(5) mentioned can be written in terms of the Wiener-Hopf operators ε + q and ε − q (see the paper [14]): Further on, consequently calculating f n (y, V (n, k)) = p k ∆t d (n,k) f k d n (y) + p k ∆t u (n,k) f ku n (y) for n = N − 1, .., 0, k = 0, ..., n,, we can get the approximate values for the sought-for functional (2) after returning to the initial notations.
The Wiener-Hopf operators can be represented in the form of pseudo-differential operators [14]: Unlike a much simpler case of Heston model [18], the presence of jumps does not allow to use closed formulae for factors unknown at the moment. Analytical formulae for the factors look as follows: The constants ω + and ω − are so that ω − < 0 < ω + play roles of parameters here and are being set to retain convergence of the respective integrals depending on Lévy process parameters. In the paper [15] a general and convenient for a numerical implementation representation of Wiener-Hopf factors was derived. The function φ + q (ξ) allows analytical continuation to a halfplane ξ > ω − and can be represented as follows: Analogously, φ − q (ξ) allows an analytical continuation to a halfplane ξ < ω + and can be represented as follows: To calculate integrals we use an algorithm of Fast Fourier Transform.

Results and discussion
We provide the results of numerical experiments, that demonstrate the convergence and accuracy of the method proposed. We will refer to our method as M&WHFa, having in mind the use of a Markov Chain approximation and the approximation for Wiener-Hopf factors. To be able to compare the result with the existing methods, we use the European down-and-out put in Heston model with the parameters set as follows: V 0 = 0.01, κ = 2, θ = 0.01, σ = 0.2, ρ = 0.5, r = 0.095. Option parameters are: T = 1, K = 100, H = 90. We also have closed formulae expressions for the Wiener-Hopf factors for the Heston model. So for comparison we will use the algorithm [18] that uses them for computation (which we will refer to as M&FWHFe). As a benchmark for all methods we used a plain Monte-Carlo generator that uses a Mersenne Twister as a pseudo-random generator.
The advanced Alfonsi [8] MC results were obtained using the code implemented into the program platform Premia.
To perform the experiment we used a PC with a configuration: Intel CPU, 2.4 GHz, RAM 8 Gb, Windows 10. The program implementation was made using a scientific Python stack (numpy, scipy and pandas packages). The algorithm parameters for the M&FWHe are as follows: N is the number of time steps, L is a scale coefficient, which controls the localization error of the FFT, and M is a number of points, meant to control its the discretization error. "Stock price" is a price of an underlying asset.
The algorithm parameters for the M&FWHa are as follows: N is the number of time steps, h is a distance between the points on a logarithmic scale, which controls the localization error of the FFT, and M is a number of points, meant to control its the discretization error. "Stock price" is a price of an underlying asset. Table 1 illustrates the convergence of the method, Table 2 and Table 3 shows absolute and relative error values, respectively. It's therefore natural to assume that similar accuracy can be achieved in models where no closed formulae representations for factors is known.

Conclusion
A general method for numerical solving of a special class of problems with barriers for partial integro-differential equations with variable coefficients arising in financial mathematics and other applications has been developed. The numerical method proposed is based on efficient formulae of approximate Wiener-Hopf factorization and allows fast calculation of values of functional of kind (2) for a wide class of models that can be represented as a diffusional process and a non-gaussian Lévy process and does not depend on the exact type of jump model.