Asymptotic Expansion and Weak Approximation for a Stochastic Control Problem on Path Space

The paper provides a precise error estimate for an asymptotic expansion of a certain stochastic control problem related to relative entropy minimization. In particular, it is shown that the expansion error depends on the regularity of functionals on path space. An efficient numerical scheme based on a weak approximation with Monte Carlo simulation is employed to implement the asymptotic expansion in multidimensional settings. Throughout numerical experiments, it is confirmed that the approximation error of the proposed scheme is consistent with the theoretical rate of convergence.


Introduction
Constructing an efficient algorithm for the following stochastic control problem associated with a relative entropy minimization (1) is an interesting topic in various areas, such as probability theory, statistical physics, economics and financial mathematics: Problem (1) appears in the risk-sensitive stochastic control problem, as studied in [1][2][3][4][5][6], where the optimal control is given by minimizing the cost depending on the risk-sensitivity of the policy maker.One of applications related to the problem (1) is the rare event simulation [7,8] in statistical physics, in which accurate approximations of rare event probabilities are studied.In the rare event simulation, importance sampling techniques are proposed by solving (1) through the variational representation based on the large deviation theory (see [9]).Moreover, the relation between the optimal control and data assimilation problems are discussed in [10].
In particular in finance, (1) is closely related to pricing and hedging problems in utility indifference pricing in incomplete market (see [11][12][13][14][15], for example).Since there is no closed-form solution for the stochastic control problem of utility indifference pricing in most cases, various numerical methods for computing indifference prices have been developed.For example, in [11], the mean-variance expansion for utility indifference pricing is proposed by using an expansion approach through Girsanov transformation.In [12], the author provides an alternative approach to the analysis of [11] by using the asymptotic expansion of the corresponding quadratic backward stochastic differential equation.The mean-variance expansion proposed in [11] is generalized in [14] to cover a multidimensional path-dependent payoff in Itô process markets.In [15], the authors extended the results of [11,14] for the case of non-smooth payoffs and apply pricing problems of power derivatives.
In the implementation of the mean-variance expansion of [11,14,15] numerically, a simple approach for computing the mean and the variance terms will be the use of the Euler-Maruyama discretization scheme for stochastic differential equations (SDEs).However, this requires many number of time steps to obtain an accurate result, since it is a first-order time discretization scheme.In other words, for a small number of time steps n, the error term by the Euler-Maruyama discretization may affect the total approximation error, including the mean-variance expansion error and the discretization error.Thus, it is important to improve the convergence rate of approximations for the mean and the variance terms in order to construct an efficient algorithm.
There have been extensive studies on asymptotic expansion methods for small noise diffusions with Malliavin calculus (for instance, [16][17][18]).Moreover, by extending these results, high-order discretization methods for SDEs are developed in various papers (for example, [19][20][21][22][23]).In particular, ref. [24] introduced a new high-order approximation method with respect to a small noise parameter ε and a number of discretization time steps n and implemented the method by deep learning.
In this paper, we show a precise error estimate of the mean-variance expansion of the stochastic control problem under various conditions on functionals on path space based on asymptotic expansion and Malliavin calculus.In particular, we prove the novel fact that the expansion error depends on the regularity of a target functional, which is an extended result of [11,14,15].Then, we implement the mean-variance expansion by using the asymptotic expansion and weak approximation to achieve the high-order approximation error with respect to ε and n based on [24].Numerical experiments confirm the theoretical convergence rate of the proposed method.
The organization of the paper is as follows.After introducing the notations and settings, we provide the main theorem and the approximation method in Section 2. Section 3 shows numerical examples of the proposed method.We conclude the paper in Section 4.

Asymptotic Expansion and Weak Approximation of Stochastic Control Problems
Let C ∞ b (R n ; R m ) be the space of infinitely continuously differentiable functions f : R n → R m with bounded derivatives of all orders.We write be the Borel field over Ω, and let P be the Wiener measure P : B(Ω) → [0, 1].Let F be the completion of B(Ω) with respect to P. Let W = {W t } 0≤t≤T be a d-dimensional Brownian motion on the probability space (Ω, F , P).Let {F t } 0≤t≤T be the filtration generated by W. We assume that {F t } 0≤t≤T contains the P-null sets of F .For a random variable Y : Ω → R N on the probability space (Ω, F , P), let E[Y] denote the expectation of Y and let Var[Y] denote the variance of Y and let ∥X∥ p := E[|X| p ] 1/p , for p ≥ 1.We define the space For more details on Malliavin calculus, see [25,26].We consider an N-dimensional diffusion driven by W: for 0 ≤ t ≤ s ≤ T: where b, . ., d and ε ∈ (0, 1].We assume that σ = [σ 1 , . . ., σ d ] satisfies the uniform elliptic condition.For notational simplicity, we write X x,ε t for X 0,x,ε Let γ > 0. It is known that the free energy of the small noise diffusion has the variational representation: where X x,ε, √ γh is a stochastic system with a control process h ∈ A: The main result is given as follows.
Theorem 1.It holds that: where: Remark 1. Theorem 1 provides a sharp asymptotic expansion for the solution of the stochastic control problem for the small noise diffusion, while the direct estimate of the left-hand side of (6) can cause inefficient computation, which is reported in [7,8].In particular, Theorem 1 provides the theoretical approximation order with respect to both γ and ε for each class of test functions f , which cannot be obtained from the asymptotic analysis in the context of the risk-sensitive control problems in [1][2][3] and the indifference pricing problems in [11,14,15].In the proof of Theorem 1 below, we will take another approach to show the sharp asymptotic expansion bounds (7), and Malliavin calculus plays a crucial role in the error estimate.
Remark 2. In the utility indifference pricing problems, γ is regarded as the risk-aversion parameter of an investor's exponential utility function U(x) = −e −γx and is typically assumed to be small as γ ≈ 0 (see [11,14,15] for more details), which is a natural setting that the investor is not far from risk-neutral (γ = 0 corresponds to the case that the investor is risk-neutral).Thus, the mean-variance expansion is interpreted as the expansion around the sum of the risk-neutral price E[ f (X x,ε T )] and the risk-aversion discount effect − γ 2 Var[ f (X x,ε T )].Theorem 1 tells us that the expansion error depends not only on the risk-aversion parameter γ but also on the smoothness of the payoff function f and the small noise parameter ε, which is a significant information in computing indifference prices in practice.
Proof of Theorem 1.We introduce a perturbed process with δ > 0: in order to expand the minimization problem: which corresponds to (4) if we set δ = √ γ.For notational simplicity, hereafter, we assume N = d = 1 without loss of generality.
The expansion of X x,ε,δh t in D ∞ (Ω) is calculated in the following way: and so on.We introduce the Jacobian of x → X x,ε,δh , i.e., Y x,ε,δh t = ∂ ∂x X x,ε,δh t , whose dynamics are: We will use a notation Y x,ε t = Y x,ε,0 t .The first-order term of the expansion of E[ f (X x,ε,δh T )] with respect to δ is given by: where D t F, 0 ≤ t ≤ T represents the Malliavin derivative process of F ∈ D ∞ (Ω) (for more details, please see [26]).Then, we have the following expansion: where for some C > 0, p, q ≥ 1 independent of f , ε and δ, through the Malliavin integration by parts formula of (2).By the Itô formula, it holds that: where , which corresponds to the Clark-Ocone formula (see [26] for more detail).We should note that if f is a sufficiently smooth function, Remark that in the case N ̸ = d, we have a similar but a more general representation of (18) under the uniform ellipticity.
If we take a sequence of functions which approximates a Schwartz distribution which is regarded as a bounded measurable function, we have that there exists C > 0 such that: and if we take a sequence of functions which approximates f ∈ C Lip (R), by (17), there exists C > 0 such that: Furthermore, it is obvious to apply (17) for the case f is smooth to obtain the desired estimate.To summarize the above discussion, we have the following gradient estimate for the diffusion semigroup which depends on the smoothness condition on f : there exists C > 0 such that: for all s ∈ [0, T].Therefore, we have: and By taking h as h s = −δ(∇P ε T−s f )(X x,ε s )σ ε (X x,ε s ), and combining with ( 13)-( 16) and ( 21), we have: with the error: Finally, setting δ = √ γ, we have: Remark 3. We comment on the advantages of the Malliavin calculus approach (the asymptotic expansion approach [17-19] based on the Watanabe theory [16]) taken in the current paper.While (11) can be obtained from both the Girsanov transform approach [11] and the Malliavin calculus approach, the error estimate (24) and the result of Theorem 1 come from only the latter approach.
Although the Girsanov transform approach is useful to derive the approximation itself, it only shows the error bound of the mean-variance expansion with respect to γ as in [11].On the other hand, the Malliavin calculus approach provides a sharp error bound not only with respect to γ but also ε depending on the smoothness of the test function f .Hence, throughout the proof in Theorem 1, we adopted the unified derivation for the approximation term (11) and the residual term(s) ( 13)-( 15) through Malliavin calculus.Moreover, the Malliavin calculus approach will be a powerful tool to approximate the mean and variance terms of the expansion in Theorem 1.We will see the usefulness of the approach in the following.
In order to implement the asymptotic expansion of Theorem 1 numerically in multidimensional settings, we efficiently approximate the mean and the variance terms by a weak approximation method for the SDE (3).
We expand the N-dimensional diffusion process X as follows: for 0 ≤ t ≤ s ≤ T: Here, we define Xx,ε,n t i+1 as: Xx,ε,n , Xx,ε,n where Xt,x,ε s , 0 ≤ t ≤ s ≤ T is given by: Moreover, we introduce the weight W ε,n T as: where ϑ x,ε t,s satisfies that there exists C > 0 such that: For more details on the derivation and the explicit form of the weight ϑ x,ε t,s , see [24].Using the scheme Xx,ε,n T and the weight W ε,n T , we have the following approximation whose property again depends on the regularity.Corollary 1.It holds that: where: Proof of Corollary 1.By [24], each expectation is discretized with the order O(1/n 2 ) as follows: with: through the applicability of the Malliavin integration by parts in the global error analysis of the weak approximation analysis according to the regularity of f .Then, combining (33) with Theorem 1, the assertion is proved.

Remark 4.
As in the proof of Theorem 1, Malliavin integration by parts formula plays an important role to prove Corollary 1.The detail proof of the small noise expansion error is shown in [17,18] and the global error analysis of weak approximation error is provided in [19][20][21][22][23].These results are essential to show the precise error estimate (34) depending on the regularity of the test function, which is an extension of the error estimate of [24].

Numerical Examples
This section provides numerical experiments to show the validity of the proposed algorithm for indifference pricing problems.
Let (Ω, F , P) be a Wiener space (which is appropriately chosen in each subsection below) on which a Brownian motion is defined.We regard P as the physical probability measure.Let M be the set of equivalent martingale measures.Let U : R → R be the investor's utility function given by U(x) = − exp(−γx), γ > 0.

Indifference Pricing under Black-Scholes Model with a Lipschitz Payoff Function
We consider 2d-dimensional SDE (d-tradable assets S s = (S s,1 , . . ., S s,d ) and d-nontradable assets X x,ε = (X x,1,ε , . . ., X x,d,ε )): for i = 1, . . ., d, µ S , µ X , ρ ∈ R and σ S , σ X > 0, where W = (W 1 , . . ., W 2d ) is a P-dimensional Brownian motion.The model of ( 35) and (36) referred to as the Black-Scholes model is widely used in financial institutions.We define Q ∈ M by: for m = r−µ S σ S and h ∈ A. Under a probability measure Q, we can rewrite the above SDE as: for i = 1, . . ., d, where We consider the case of a basket option of the nontradable assets, i.e., we set the payoff function We assume that the riskfree rate r = 0 for simplicity.Here, the buyer utility indifference price p is given by: We approximate the indifference price by the proposed method.Since f ∈ C Lip (R d ), by Theorem 1, it holds that: In order to estimate the expansion error of (40), we compute the both sides by using the explicit solution of X x,ε,0 obtained by the Itô formula and Monte Carlo simulation with M = 10 8 paths.Moreover, the approximation errors of the mean and variance terms are given by Corollary 1 as where { Xx,ε,n t i } i=0,1,...,n is the approximation process for X x,ε = X x,ε,0 introduced by ( 27) and (28).To check the convergence rate of (41), we employ Monte Carlo simulation with M = 10 8 paths to implement the proposed method, where the reference value (the left hand side of (41)) is obtained by the Itô formula and Monte Carlo simulation with M = 10 8 paths.
First, we estimate the expansion error of (40) with respect to ε and γ.The meanvariance expansion (40) is referred to as "MV-expansion" in the following figures.
Next, we estimate the weak approximation error of (41).The proposed second-order weak approximation method is referred to as "WA 2nd" in the following figures and tables.For comparison, we also compute the approximation error by using the Euler-Maruyama scheme, referred to as "EM".The approximation error of ( 41) is plotted in Figure 3.The figure shows that the error of "WA 2nd" decreases rapidly as the number of time steps increases compared to "EM", which means that the proposed method achieves the second-order accuracy with respect to the number of time-steps n.The results are summarized in Table 1.
Table 1.Numerical error of (41) for each ε with γ = 0.01 under the one-dimensional Black-Scholes model.
By the figures and the tables, we confirm that the proposed method achieves the theoretical rate of convergence with respect to γ, ε and n.
First, we estimate the expansion error of (40) with respect to ε and γ.We perform the same experiment with the parameter γ = 0.01, 0.02, 0.04, 0.08 for the fixed ε = 0.4.The result summarized in Figure 5.In Figures 4 and 5, we can check that the expansion error achieves the theoretical rate of convergence of O(γ 2 ε 3 ).
Next, we estimate the weak approximation error of (41).The discretization error is plotted in Figure 6.The figure shows that the proposed method provides an accurate approximation for (41) with a small number of time steps compared to the Euler-Maruyama scheme.The results are summarized in Table 2.
By the figures and the table, we confirm that the proposed method achieves the theoretical rate of convergence.
In Figures 7 and 8, we can check that the expansion error achieves the theoretical rate of convergence of O(γ 2 ε 3 ).
Next, we estimate the weak approximation error of (41).We set the value we computed in the previous experiment as the reference value of the right hand side of (41).The discretization error is summarized in Table 3.
The table shows that the proposed method approximates (41) more accurately than the Euler-Maruyama scheme with a small number of time steps.Throughout the numerical experiments, we confirm that the proposed method achieves the theoretical rate of convergence consistently with Theorem 1 and Corollary 1 in multidimensional settings with a Lipschitz continuous test function f .

Indifference Pricing under Constant Elasticity Model (CEV Model) with a Bounded Measurable Payoff Function
We consider the following two-dimensional SDE: where W = (W 1 , W 2 ) is a P-dimensional Brownian motion.The dynamics of ( 43) is called the constant elasticity of variance model (CEV model), which is a generalized model of the Black-Scholes model.We define Q ∈ M by: for m = r−µ S σ S and h ∈ A. Under a probability measure Q, we can rewrite the above SDE as: where As the previous section, we assume that the riskfree rate r = 0 for simplicity.We consider the case of a digital option of the nontradable asset, i.e., we set the payoff function f : R → R as f (x) = 1 {x>K} .
Here, the buyer utility indifference price p is given by: We approximate the indifference price by the proposed method.Since f ∈ B b (R), it holds that: To estimate the expansion error of (49), we compute both sides by using the Euler-Maruyama discretization scheme with n = 2 10 time steps and Monte Carlo simulation with M = 2 × 10 9 paths.The mean-variance expansion (49) is referred to as "MV-expansion (EM n = 2 10 )" in the following figures.
Moreover, the approximation errors of the mean and variance terms are given by Corollary 1 as: where { Xx,ε,n t i } i=0,1,...,n is the approximation process for X x,ε,0 introduced by ( 27) and (28).To check the convergence rate of (50), we employ Monte Carlo simulation with M = 2 × 10 9 paths to implement the proposed method, where the reference value (the left hand side of (50)) is obtained by using the Euler-Maruyama discretization scheme with n = 2 10 time steps and Monte Carlo simulation with M = 2 × 10 9 paths.
We set the parameters as T = 2.0, x = 100.0,K = 100.0,µ S = µ X = 0.0, σ S = σ X = 0.3, ρ = 0.0, β = 0.5.We estimate the approximation error of (49) with respect to γ. Figure 9 plots the results for ε = 0.1, 0.2, 0.4, 0.6, 0.8, 1.0 for the fixed γ = 0.025.In Figure 9, a linear regression line (which is referred to as "Regression" in Figure 9) is added to confirm the empirical convergence rate.We can check that the coefficient of ε is quite small and the regression line seems to be flat.By the experiment, we confirm that the expansion error is consistent with Theorem 1.
In Figure 10, we can check that the expansion error achieves the theoretical rate of convergence of O(γ 2 ).Since f ∈ C Lip (R d ), by Theorem 1, we have: We check the expansion error of (59) by computing the both sides by using the Euler-Maruyama discretization scheme with n = 2 10 time steps and Monte Carlo simulation with M = 10 8 paths.The mean-variance expansion (59) is referred to as "MV-expansion (EM n = 2 10 )" in the following figures.
Moreover, the approximation errors of the mean and variance terms are given by Corollary 1 as: where { Xx,ε,n t i } i=0,1,...,n is the approximation process for X x,ε,0,0 introduced by ( 27) and (28).To check the convergence rate of (60), we implement the proposed method by Monte Carlo simulation with M = 10 8 paths, where the reference value (the left hand side of (60)) is obtained by using the Euler-Maruyama discretization scheme with n = 2 10 time steps and Monte Carlo simulation with M = 10 8 paths.
Next, we estimate the weak approximation error of (60).The approximation error of (60) is plotted in Figure 13.The figure shows that the error of "WA 2nd" decreases rapidly as the number of time steps increases compared to "EM", which means that the proposed method achieves the second-order accuracy with respect to the number of time-steps n.The results are summarized in Table 5.Throughout the numerical experiments, we confirm that the proposed method provides a consistent result with Theorem 1 and Corollary 1 even when f is a non-smooth function.

Conclusions
In the paper, we showed a precise error estimate of the mean-variance expansion of the stochastic control problem under general conditions on the terminal test function based on asymptotic expansion and Malliavin calculus.In particular, we proved that the expansion error depends on the smoothness of the test function, which is an extension of [11,14,15].Moreover, an efficient algorithm has been introduced based on the asymptotic expansion method and weak approximation for the small noise diffusion.Numerical experiments confirmed the theoretical convergence rate of γ, the small noise parameter ε and the number of time-steps.For future work, it is worth considering to apply the proposed method to important problems in statistical physics, such as the rare event simulation studied by [7,8].