Numerical Solution of Nonlinear Backward Stochastic Volterra Integral Equations

: This work uses the collocation approximation method to solve a speciﬁc type of backward stochastic Volterra integral equations (BSVIEs). Using Newton’s method, BSVIEs can be solved using block pulse functions and the corresponding stochastic operational matrix of integration. We present examples to illustrate the estimate analysis and to demonstrate the convergence of the two approximating sequences separately. To measure their accuracy, we compare the solutions with values of exact and approximative solutions at a few selected locations using a speciﬁed absolute error. We also propose an efﬁcient method for solving a triangular linear algebraic problem using a single integral equation. To conﬁrm the effectiveness of our method, we conduct numerical experiments with issues from real-world applications.


Introduction
Backward stochastic differential equations (BSDEs) represent stochastic differential equations with terminal conditions.The existence and uniqueness of solutions for BS-DEs have been proven by Pardoux and Peng, who also developed a general nonlinear BSDE [1,2].Backward stochastic differential equations have numerous applications in finance, stochastic games, and optimal control.The concept of backward stochastic differential equations has been extended to include backward stochastic Volterra integral equations (BSVIEs), which depend on two specific time moments for their drift and diffusion coefficients.Nonlinearities are present in the general of BSDEs as follows.( We will introduce a primary space (Ω, F , {F t } t≥0 , P) corresponding to a total probability space with a filtration, where {F t } t≥0 meets the usual criteria (i.e., it is right continuous and F 0 contains all P-null sets).In constrast, {B(t)} t∈[0,T] defines the Wiener process.The terminal condition ξ is an F T -measurable random variable, and the driver g is a progressively measurable function.The adapted solution of the BSDE (1) is the pair (Y(•), Z(•)) of the adapted processes that satisfy (1).The adapted solution's second component, Z(•), is known as the martingale integrand.
Our investigation is inspired by the method for estimating the BSDEs' adapted solutions [3].We recommend researching backward stochastic Volterra integral equations (BSVIEs) in light of the most recent research of [4][5][6].Pardoux and Peng began this research over a decade ago [7].According to Lin [8], the modified solutions were studied as existence and uniqueness problems under global Lipschitz conditions.The global Lipschitz condition on drift has been eased by Aman and N'zi [9].For a comprehensive explanation of the theory and applications of BSDE (1), including stochastic controls and mathematical finance, the reader might consult El Karoui, Peng, and Quenez's overview paper [10].The emergence of BSVIEs of the form has significantly developed BSDEs.As a natural progression from BSDEs, BSVIEs can be represented as follows. (2) In the literature, (2) is referred to as a Type-II BSVIE.[11][12][13][14][15][16][17][18][19][20] and the references therein).
Approximations for adapted M solutions of Type-II backward stochastic Volterra integral equations were studied, where backward stochastic differential equations converge to the adapted M solution of the original equation [21].In addition, the convolution method has been extended to solve the conditional expectation to solve BSDEs numerically, and a generalized θ scheme has been applied to discretize the backward component [22].
The numerical approximations problem for Type-II BSVIEs has been completely open.There needs to be more quantitative interest in BSVIEs.Hence, with the aid of block pulse functions and their stochastic operational matrix of integration, backward stochastic Volterra integral equations can be effectively solved.These equations can then be reduced to a linear lower triangular system, which can then be solved by forward substitution (See, e.g., [23][24][25][26][27][28][29][30][31][32][33] and the references therein).
The primary characteristic of BSVIEs (2) is that they include memories, which are more accurate to reality.We seek the unknown pair (Y(•), Z(•, •)), where Y(•) and Z(t, •) are adapted for each t ∈ [0, T].In the above, the free term ψ(•), also known as the terminal condition, is allowed to be only a B([0, T]) ⊗ F T -measurable stochastic process (not necessarily F -adapted).Here, B([0, T]) represents the Borel σ field of [0, T].The generator or the driver of the BSVIE is a given map g(•), which can be deterministic or random.The coefficient g(•) is dependent on both t and s, and it g(•) depends not only on Z(t, s), but also on Z(s, t).The drift generally depends on Z(t, s) and Z(s, t).In the case where the driver g is independent of the term Z(s, t), the BSVIE becomes the following: For convenience, we have rewritten the following BSVIE: The structure of this work is as follows.Section 2 covers the basic characteristics of block-pulse functions and an integration operational matrix approximation.In Section 3, the stochastic integration operational matrix is presented.Section 4 solves stochastic Volterra integral equations using the stochastic integration operational matrix via collocation approximation.Section 5 presents an analysis of the solution's general error estimated regularity properties of the solution.In Section 6, we offer numerical results and use numerical examples to demonstrate the accuracy of the suggested approach.

Block-Pulse Functions BPFs
The block-pulse function (BPF) Φ i over the unit interval [0, T) is defined as follows: , and h = T m .The block-pulse functions have the following properties: (1) Disjointness: The BPFs are disjointed with each other in the interval t ∈ [0, T).
(3) The third property is completeness: For every f ∈ L 2 [0, T), when m → ∞, Parseval's identity holds, that is: where The set of function can be described by an m vector.
Therefore, we can write the relationship between BPFs and their integrals in the following matrix form.
The above representation and disjointness property follows and where D F usually denotes a diagonal matrix whose diagonal entries are related to a constant vector

Function Approximations
A real bounded function f (t), where f (t) ∈ L 2 [0, T), can be expanded into a blockpulse series as where f i is the block-pulse coefficient with respect to the ith It can be similarly expanded with respect to BPFs such as where Ψ(s) and Φ(t) are m 1 and m 2 dimensional BPF vectors, respectively, and For convenience, we put m 1 = m 2 = m.

Integration Operational Matrix
From [34], we have: where the operational matrix of integration is given by Therefore, each function f (t) can be expressed as

Stochastic Integration Operational Matrix
The It ô integral of each single BPF φ i (t) can be computed as follows: Now by expressing T 0 φ i (s)dB(s) in terms of the BPFs, we have Therefore, In this case, the stochastic operational matrix of integration can be expressed as follows: where Therefore, the It ô integral of every function f (t) can be manipulated as follows: By approximating the functions Y(t), g(t, s), z 1 (t, s), z 2 (t, s) via BPFs by relations, we have

Implementation in Stochastic Integral Equation
Using the block-pulse operational matrices, we first find the collocation approximation to the functions z 1 (t, s) and z 2 (t, s) for drift and diffusion, respectively, which are defined by From Equations (4) and ( 21), we obtain and We can approximate z 1 (t, s) and z 2 (t, s), and we can assume that g(t, s) is a function of two variables via the block-pulse series as follows: such that m vectors Z 1 , Z 2 , and m × m matrix G are the block-pulse coefficients of z 1 (t, s) and z 2 (t, s) and g(t, s), respectively.By substituting (24) in (22), we obtain In addition, the It ôs integral of ( 22) can be written as where Z1 = diag (Z 1 ), and Z2 = diag(Z 2 ).By substituting ( 25) and ( 26) into ( 23), as well as by replacing with =, we obtain The collocation method with (27), in m nodes t j , s j = j m+1 , j = 1, ..., m, is used for determination of the following: After solving the nonlinear system (28), we obtain Z 1 and Z 2 .Then, the result Y(t) of ( 22) is approximated as follows:

General Error Estimate
In this section, we will provide the general estimate used to determine the convergence.We shall first rely the following presumptions: with respect to y and z, and g possesses first and second partial derivatives that are both continuous and uniformly bounded.
There exists a constant p 0 > 2 and L such that, for continuous process ψ(•), there is an F-adapted value.
The regularity of (Y(•), Z(•, •)) (Theorem 3.7 and 4.1 in [5]) provide the following result of well-posedness of the BSVIE.Theorem 1.A unique solution (Y(•), Z(•, •)) to the BSIVE is admissible under assumptions (I)-(III).In addition, the following estimate holds: Furthermore, the following BSVIE has an adapted solution The convergence speed is calculated using the following result: Lemma 1.For each t, t 0 ∈ [0, T], it holds under assumptions (I)-(III) that where C is a constant.

Numerical Results
In this section, we will provide two numerical examples to illustrate the results obtained in Sections 3 and 4. All computations were carried out in MATLAB R2018a, with a precision of 2.22 × 10 −16 .To compare the values of the approximate and exact solutions at selected points, we used the definition of the absolute error, which is as follows: where X i represents the exact solution, and Xi represents the approximate solution.
The results obtained for k = 0.5, σ = 0.2, and θ(t) = 4 sin(t) in this example are given in Table 1.The approximate and exact solutions' graphs and the absolute error for t = 0.7; k = 0.6, σ = 0.8, m = 100, N 1 = 20, N 2 = 200, and n = 6 are plotted in Figures 1-3, respectively.The accuracy of the generalized absolute error in Table 1 depends on the parameters X i .The error decreases as time steps increase.Example 2 ([24]).Consider the following nonlinear stochastic Volterra integral equation with the exact solution Assume that y(0) = 0.5, r = 0.7, K = 2, n = 10, β = 0.4, T = 0.9, and t, s ∈ [0, 0.9].Table 2 displays the numerical results for various values of m, including the absolute error for the figures-both exact and approximate-for the parameters y(0) = 0.4, r = 1, K = 0.8, n = 100, β = 1, and T = 7. Figures 4-6 show the exact and approximate solutions, as well as the variation trend.

Conclusions
The current study focuses on the so-called Type-II BSVIEs, where the coefficient g(•) is dependent on both t and s, and g(•) depends not only on Z(t, s), but also on Z(s, t).This paper proposed a collocation approximation method to predict an unknown function.By implementing Newton's method, we solve the BSVIEs using block-pulse functions and the corresponding stochastic operational matrix of integration.In addition, we have included examples that demonstrate estimate analysis while highlighting the separate convergence of the two approximating sequences.We also measured the solutions against the values of the exact and approximate solutions at a few selected locations using a specified absolute error.According to the collocation approximation solutions, the issues raised in the work might be applied to Type-I BSVIEs.However, this strategy requires an entirely new methodology, and it is left to future studies to determine the error by computing conditional expectations.It might be the focus of some future research.

Author Contributions:
Conceptualization, E.Y.K.; methodology, M.S. and E.Y.K.; software, M.S.; validation, M.S.; formal analysis and investigation, X.Z.; resources, M.S. and X.Z.; writing-original draft preparation, E.Y.K. and M.S.; supervision, X.Z.All authors have read and agreed to the published version of the manuscript.Funding: This work was supported by the Zhejiang Normal University Postdoctoral Research Fund (Grant No. ZC304022938), the Natural Science Foundation of China (Project No. 61976196) and the Zhejiang Provincial Natural Science Foundation of China under Grant No. LZ22F030003.
Unlike a Type-I BSVIE, solving a Type-II BSVIE requires an extra constraint on the term Z(t, s), where 0 ≤ s ≤ t ≤ T for the equation to demonstrated well-posedness.Researchers have developed the adapted M solution, which was inspired by the duality principle in stochastic control problems of stochastic Volterra integral equations, to solve the Type-II BSVIE.This equation is essential for studying stochastic control and mathematical finance problems.Researchers have used BSVIEs to calculate dynamic risk estimates for position operations and to examine dynamic capital allocations.BSVIE solutions can describe time-inconsistent recursive utility processes of general discounting, and they are strongly linked to time-inconsistent stochastic control problems.Several researchers have proposed differentiability results, investigated stochastic control problems for SVIE and BSVIE systems, and proven numerous comparison theorems for adapted solutions and adapted M solutions to BSVIEs in multidimensional Euclidean spaces.Notably, BSVIE theory is path-dependent, and numerical elements have also been considered (see, e.g.,[5, The Hull-White Model: Hull and White investigated Vasicek model extensions that perfectly fit the basic term structure in 1990.A single-factor interest rate model is the Hull-White model.The short interest rate is assumed to have a normal distribution in this model, and there is no arbitrage assumption.The short interest rate, therefore, satisfies the stochastic

Table 1 .
Mean, standard deviation, and mean confidence interval for error.The graph of absolute error function for Example 1.The trajectory of the approximate solution and exact solution of Example 1. Variation trend of absolute error of Example 1.

Table 2 .
The trajectory of the approximate solution and exact solution of Example 2. Mean, standard deviation, and mean confidence interval for error.
Figure 4.The graph of absolute error function for Example 2.