Approximation of SDE solutions using local asymptotic expansions

We develop an asymptotic expansion for small time of the density of the solution of a non-degenerate system of stochastic differential equations with smooth coefficients, and apply this to the stepwise approximation of solutions. The asymptotic expansion, which takes the form of a multivariate Edgeworth-type expansion, is obtained from the Kolmogorov forward equation using some standard PDE results. To generate one step of the approximation to the solution, we use a Cornish–Fisher type expansion derived from the Edgeworth expansion. To interpret the approximation generated in this way as a strong approximation we use couplings between the (normal) random variables used and the Brownian path driving the SDE. These couplings are constructed using techniques from optimal transport and Vaserstein metrics. The type of approximation so obtained may be regarded as intermediate between a conventional strong approximation and a weak approximation. In principle approximations of any order can be obtained, though for higher orders the algebra becomes very heavy. In order 1/2 the method gives the usual Euler approximation; in order 1 it gives a variant of the Milstein method, but which needs only normal variables to be generated. However the method is somewhat limited by the non-degeneracy requirement.


Introduction
In this paper we describe an approach to the stepwise approximation of solutions of stochastic differential equations using an asymptotic expansion for each step. Consider an SDE dx(t) = a(t, x(t))dt + B(t, x(t))dW (t) (1) where x(t) ∈ R q , a(t, x) ∈ R 2 , B(t, x) is a q×d matrix and W is a d-dimensional standard Brownian motion, with an initial condition x(0) = x (0) . The first step of a stepwise approximation with step h requires the simulation of an approximation to x(h). The standard approach, as described in [9], is based on stochastic Taylor expansions, giving rise to the Euler aproximation x(h) ≈ x (0) + ha(0, x (0) ) + B(0, x (0) )W (h)) and the Milstein approximation where B = (b i j ) and ρ i jk = q m=1 b m j ∂b ik ∂ x m , as the lowest-order approximations. The difficulty with this method is that the stochastic integrals involved are hard to generate, already in the case of the double integrals appearing in (2). In [3] an approach to resolving this difficulty is described, under the assumption that B(0, x (0) ) has rank q, using a perturbation method to generate approximations to the stochastic integrals. This gives an analogue of (2) where X 1 , · · · , X d are independent N (0, h) random variables.
Here we use a different approach, which involves first deriving a asymptotic expansion to the distribution of x(h). For small h a first approximation to this distribution is the normal distribution N (x 0 , h (0) ) where (0) = B(0, x 0 )B(0, x 0 ) t , which suggests an Edgeworth-type asymptotic expansion of the form for the density function f h of h −1/2 (x(h)−x 0 ) where φ is the density of N (0, ) and the S j are polynomials on R q . Provided B(0, x 0 ) has rank q, so that is nonsingular, and under smoothness conditions on a and B, we derive such an expansion from the Kolmogorov forward equation for the density, and make precise in what sense it is an asmptotic expansion of f h . Given the expansion (4), one can then look to approximate x(h) in distribution by a random variable of the form x 0 +h 1/2 (X + k j=1 h j/2 p j (X )) where X is an R q -valued random variable with N (0, (0) ) distribution, and the p j are R q -valued polynomials in q variables. It is shown in [4] that any random variable of this form has a density of the form (4), and the relation between the polynomials p 1 , · · · , p k and the sequence (S j ) in (4) is described there in algebraic terms -it is in effect a multvariate version of the relation between Edgeworth and Cornish-Fisher expansions in one variable.
We show that, in principle at least, this programme can be carried out to give approximations to arbitrary order provided the coefficients a and B are smooth and as long as the essential condition that B has rank q continues to hold. The error in the approximation to x(h) is assessed by a suitable coupling between x(h) and the approximation, using methods from the theory of optimal transport. Sections 2 and 3 give background material on optimal transport and on Edgeworth and Cornish-Fisher-type expansions. Section 4 describes the construction of the expansion (4) and the sense in which it is an effective asymptotic expansion. The application to stepwise approximation schemes is described in Sect. 5 for a single step, and then applied to multistep approximation in Sect. 6, where the main approximation boundss are stated and proved proved. Section 7 concludes with a comparison of this approximation scheme with standard strong and weak approximation methods, indicating that it is stronger than weak approximation but somewhat weaker than standard strong approximation, and mentions some possible directions for further work.
For general background on numerical approximation of SDE solutions see [9]. In a rather different direction we mention [6] which shows how discrete approximations can be used to prove existence of strong solutions of SDE.
I am very pleased to be able to contribute to the Special Issue in honour of Professor Gyöngy, who has contributed so much to the field of Stochastic Analysis.

Optimal transport background
We will use some notation and results from optimal transport theory. If P andP are probability measures on R q , then for p ≥ 1 the Vaserstein distance W M (P,P) is defined as inf(E|X −X | M ) 1/M where the inf is over all joint distributions on R q × R q for R q -valued random variables X andX which have marginals P andP. We sometimes abuse notation and write this as W M (X ,X ) or W M (X ,P); we will mainly consider the case when P andP are absolutely continuous w.r.t. Lebesgue measure, with densities f andf , and then we may use the notation W M ( f ,f ). One result we need is that if f andf are probability densities on R q then for any M ≥ 1. This follows from Proposition 7.10 in [15]. We will in fact need a specific coupling (i.e. joint distribution) between two random variables X andX , with densities f andf respectively, which realises (5), and one such is as follows: let Then one readily checks that X andX do have the specified densities, and a straightforward calculation gives from which (5) follows. Details can be found in [15]. Much more on this topic can be found in the books of Rachev and Ruschendorff [12] and Villani [15,16]. And a note on our spelling of 'Vaserstein': his original paper [14] is in Russian, we have used the transliteration 'Vaserstein' (rather then the alternative 'Wasserstein') from the Cyrillic alphabet as that is the one he himself has used in his English-language publications.

Polynomial perturbations of normal distributions
In this section we outline some results from [4] which we will use. Fix q ∈ N and let P denote the space of all real-valued polynomials on R q , and P q the space of R qvalued polynomial functions on R q . We also fix a positive-definite q ×q matrix . Let p 1 , · · · , p k ∈ P q . For ∈ R we define ρ : R q → R q by ρ (x) = x + k j=1 j p j (x). We are interested in the distribution of ρ (X ) where X is an R q -valued random variable with N (0, ) distribution and is close to 0. Pretending that ρ is bijective, then this distribution has a density given by By expanding ρ −1 (y) formally in powers of one can obtain a formal expansion We remark that, if p ∈ P q and det(Dp) is not identically 0, then the zero set of det(Dp)is a closed set of measrue zero, and p is locally invertible on its complement, from which one can deduce that p(X ) has a density. Applying this to ρ , one sees that det(Dρ (x)) = 1 + h( , x) where h is a polynomial, so that if is small enough it cannot vanish identically. Hence ρ (X ) has a density f for small enough.
The construction of the S j from the p j is described in [4]. One sees that for j < n, S j depends only on { p m : m ≤ j}. It follows that, given a sequence p 1 , p 2 , · · · with p j ∈ P q we obtain a well-defined sequence S 1 , S 2 , · · · . Now we introduce the notation P for the set of all sequences (u 1 , u 2 , · · · ) with u j ∈ P, and similarly P q . We write P for the set of sequences satisfying u k (y)φ (y)dy = 0 for all k. Then we define the linear mapping S : P q → P by S ( p 1 , p 2 , · · · ) = (S 1 , S 2 , · · · ) as just described. It is shown in [4] that this mapping is onto. In one dimension it is also (1-1); it gives the connection between Cornish-Fisher expansions and Edgeworth expansions. But when q > 1 it is not (1-1), meaning that there are many Cornish-Fisher expansions for a given Edgeworth expansion.
The calculation of S n from p 1 , · · · , p n is straightforward in principle, though the algebra becomes heavy as n increases. As indicated in [4], we find S 1 (y) = y t −1 p 1 (y) − ∇. p 1 (y) and The following result from [4] demonstrates the sense in which the series (1 + ∞ j=1 j S j )φ can be regarded as an asymptotic expansion for the density f of ρ (X ).

Proposition 1 With notation as above, for any n and M ≥ 1 there exists C > 0 so that
for all sufficiently small .
Here we will be concerned with the reverse direction, where we have a sequence (S 1 , S 2 , · · · ) giving an asymptotic expansion of a family ( f ) of probability densities, and wish to find a corresponding sequence ( p 1 , p 2 , · · · ) for the purpose of generating an approximation to P . Generally, if ( f ) ∈(0,1) is a family of probability densities on R q and (S n ) n∈N is a sequence in P, then we say that (S n ) is an α -sequence for (P ) if, for any M ≥ 1, there exists C > 0 so that (8) holds for all sufficiently small > 0. Note that an α -sequence is automatically an A -sequence as defined in [4] -the notion of A -sequence is slightly more general and allows for probability measures without densities, but the above notion of α sequence is somewhat simpler and suffices for our purpose here.
Then Theorem 4 of [4] gives the following: ) and (f : ∈ (0, 1)) are families of probability densities on R q having respectively an α -sequence (S 1 , S 2 , · · · ) and an α -sequence (S 1 ,S 2 , · · · ). Suppose also that, for some n ∈ N, we have S j =S j for 1 ≤ j ≤ n. Let M ≥ 1 be given. Then we can find C > 0 such that Note that Proposition 1 says that (S 1 , S 2 , · · · ) is an α -sequence for ρ (X ). The idea now is that, if we have an α -sequence (S 1 , S 2 , · · · ) for some family f , we can find a sequence ( p 1 , p 2 , · · · ) giving the same (S 1 , S 2 , · · · ), and then use the resulting ρ (X ) as an approximation to f ; Proposition 2 then gives a bound for the accuracy of this approximation. In this context we now describe, for later use, a specific coupling between f and ρ (X ) which attains the bound in Proposition 2.
To do this, first fix M ≥ 1 and let r ∈ N with r + 1 ≥ M(n + 1).
. Now if Y has density f then, using (8), the construction described at the end of Section 1 gives an explicit coupling between Y and V with E|Y − V | M ≤ C M(n+1) . Using this coupling to generate Y conditional on V we then get a joint distribution of Y , V , Z for which E|Y − Z | M ≤ C M(n+1) as required. We remark that this deduction of a coupling between Y and Z , given a coupling between V and Z , and one and one between Y and V , is an example of gluing of couplings, as discussed on page 11 of [16].
In the case M = 2 we can give a more precise estimate of the W 2 distance in Proposition 2, using Theorem 11 and equation (22) from [4] (with˜ = ). We find that, under the hypotheses of Proposition 2, we have W 2 (P , General references for Edgeworth-type expansions in several variables are [1] and [11]. For the classical theory of Edgeworth and Cornish-Fisher expansions in one variable see [7] and [8] for example. One method of constructing Cornish-Fisher expansions in several variables is given in [13]. (1), and consider the distribution of the random variable Y t = t −1/2 (x(t) − x (0) ). For the present we make the following assumption: (*) a and B are C ∞ on (0, T ) × R q , for each multi-index α the derivatives D α a and D α B are bounded, and also B has rank q everywhere and (B B t ) −1 is bounded.

Construction of asymptotic expansions
Then Y t will have a density f t and we look for an asymptotic expansion valid for small t > 0, to which we can apply the methods from section 3, with = t 1/2 . To do this we relate f t to the density of the solution of the SDE: let u(t, x) denote the density of x(t). Then u satisfies the forward equation Lu = 0 where ) and then, if we have an expansion (9) we obtain The idea is to substitute the RHS of (11) for u(t, x) in (10), then find polynomials S 1 , S 2 , · · · so that Lu = 0 holds in a formal power series sense, and finally verify that with these (S 1 , S 2 , · · · ) the expansion (9) is a valid asymptotic expansion of ( f t ).
To carry out the first step, we first note that by an orthogonal change of variables we may arrange that (0) is diagonal, with positive entries λ 1 , · · · , λ q . We also impose the requirement that a(t, x) and B(t, x) are infinitely differentiable in a neighbourhood of (0, x (0) ). Then we have Taylor expansions a i (t, σ i jkα t k+|α|/2 y α , where the sums are over all nonnegative integers k and nonnegative multi-indices α of length q. We have σ i jkα = σ jikα and σ i j00 = λ i δ i j . Then the sought expansion (11) can be written as where where Q n and R n are polynomials depending only on S k for k < n. Then substituting (12), (14) and (13) into Lu = 0 using (10), and comparing coefficients of t , and we deduce finally that is the requirement for (10) to hold in a formal power series sense, whereQ n = Q n − 2R n . Now for n ∈ N we define an operator n , acting on polynomials in q variables, by n p(y) = np(y)+ y.∇ p(y)− q i=1 λ i ∂ ii p(y), then we see that for any non-negative multi-inder α we have n (y α ) = (n + |α|)y α + g where g has degree |α| − 2. It then follows by induction on degree that the equation n p = Q has a unique polynomial solution p for any given polynomial Q. Since (15) can be written as n S n =Q n , we conclude by induction on n that there is a unique sequence of polynomials (S n ) satisfying (15).

Recursive calculation of S n
For notational simplicity, we assume here that λ i = 1 for each i (which we can always arrange by taking a non-orthogonal initial coordinate change), and then we find after some calculations that (15) can be written as where A iα and B i jα are operators on polynomials defined by and B i jα S(y) ={(α i α j − δ i j α i )y α−e i −e j − 2α i y α−e i +e j + (y i y j − δ i j )y α }S(y) and where * denotes a sum over all i, k, α with i ∈ {1, · · · , q}, k ≥ 0 and multiindices α such that 2k + |α| < n, while is a similar sum over i, j, k, α such that 0 < 2k + |α| ≤ n; here e i denotes the multi-index with i'th entry 1 and all other entries 0.
In order to apply the recurrence relation (16) for S n it may be helpful to express the RHS in terms of Hermite polynomials and note the easily verified fact that where H α (y) = H α 1 (y 1 ) · · · H α q (y q ). The algebra involved in the this inductive construction is in general quite heavy, as is often the case with Edgeworth-type expansions. The case n = 1 is relatively simple -we find that 1 S 1 (y) = 2 i a i00 y i + where = i, j,s λ i js y s (y i y j − δ is − δ js − δ i j ), with λ i js = σ i j0e s = ∂ s i j (0, x (0) ). We can write = s λ sss H 3 (y s ) + i =s (2λ isi + λ iis )H 2 (y i )y s + * λ i js y i y j y s where * denotes a sum over all cases where i, j, s are distinct. Then we see that is a linear combination of H α terms with |α| = 3, and deduce using (19) that The algebra is simpler in the case q = 1. In this case we have n p(y) = np(y) + yp n (y) − p n (y) and (16) becomes where * is a sum over non-negative integers (k, r ) such that 2k +r < n, while is a sum over pairs such that 0 < 2k+r ≤ n. From this we obtain S 1 (y) = a 00 y+ σ 01 4 H 3 (y) (which also follows from (20)) and We also note two properties of S n which follow easily by induction on n using (16): S n is odd or even according as n is odd or even, and S n has degree at most 3n.
The following theorem makes precise the sense in which the expansion (9) is an asymptotic expansion of the distribution of Y t . (1) on (0, T ) with initial condition x(0) = x (0) , and that E|x(t)| K < ∞ for all K ≥ 1. Suppose also that the coefficients a and B are C ∞ on a neighbourhood of (0, x (0) ), and that B(0, x (0) ) has rank q. Let S n be the unique sequence of polynomials satisfying (15). Then (S n ) is an α (0) -sequence for the distribution of Y t ; in other words, given n ∈ N and M ≥ 1 there exist C > 0 such that the bound

Theorem 1 Suppose x(t) is a solution of an SDE of the form
holds for all t ∈ (0, T ), where P t is the distribution measure of Y t and ν t,n is the measure with density f t,n (y) := φ (0) (y) 1 + n k=1 t k/2 S k (y) Proof We fix n and M. and first assume (*); we shall use the notation C 1 , C 2 , · · · for constants which depend only on n, M and bounds for (B B t ) −1 and for D α a and D α B for finiely many α. Let N = max(q + d, 2n + 6) define u N by the truncated . Then the fact that (10) holds in a formal power series sense implies that where V N is a polynomial in q + 1 variables. We extend ψ N and u N to W := (−∞, T ) × R q by setting then to 0 for t ≤ 0. Then using N > q + 3 one can check that ψ N is C 1 , and bounded with bounded first derivatives on W . And, as a distribution on W , we have Lu N = ψ N + δ 0 , where δ 0 is the unit point mass at (0, 0). Moreover, if 0 < τ < T , then ψ N and its first derivatives are bounded by C 1 τ Then by standard PDE results (for example using Theorem 8.10.1 in [10], with any δ ∈ (0, 1), and noting the comment after equation (8.0.2) that the condition c ≤ 0 can be avoided) we conclude that |u(t, for 0 < t < T , y ∈ R q , and hence Let R = 2M + q + 1, then we have (1 + |y|) R | f t (y)|dy ≤ C 4 , and clearly (1 + |y|) R | f t,N (y)|dy ≤ C 5 , so we have Applying Cauchy-Schwartz to (22)  To treat the general case, given a and B satisfying the hypotheses of the theorem, we findã andB satisfying (*), and such that they agree with a and B respectively for |x − x 0 | < δ, for some δ > 0. LetP t be the distribution ofỸ t . Then, for given M, we have C such that for some constant K . Then Cauchy-Schwartz gives (1 + |y|) M d|P t −P t |(y) ≤ K (n+1)/2 , which together with (24) gives (21).
We remark that the hypotheses of this theorem allow cases where the distribution P t does not have a density. For example the one-dimensional SDE dx = |x| 1/2 dW with initial condition x(0) = −1 has a solution which stays at 0 after first reaching 0, and for this solution P(x(t) = 0) > 0 for all t > 0 so P t has a point mass at 0.
In fact the application of theorem 1 in the following sections only uses the case (*), but we have included the general case in case the asymptotic expansion may have some independent interest.

One-step approximation
Now we apply the results described in the previous sections to the construction of an approximation scheme. Given an SDE of the form (1), we find the polynomials S n in the expansion (9) of the density of Y t , then use the methods of [4] as outlined in Sect. 3 above to construct a Cornish-Fisher type expansion X + t 1/2 p 1 (X ) + t p 2 (X ) + · · · , where X is an N (0, (0) ) random variable, which expansion can be truncated to give a random variable whose distribution approximates that of Y t . To illustrate the procedure we consider the relatively simple case of p 1 . We asume for simplicity that d = q, and then we assume as usual that B(0, x (0) is invertible. As before, by suitable choice of coordinates we can assume that B(0, x (0) ) = I , and the (0) = I . We recall the expression (20) for S 1 , and note that the requirement for p 1 is that y t p 1 (y)−∇. p 1 (y) = S 1 (y). One choice for p 1 given by the approximation scheme (3), as follows: the scheme (3) which corresponds to p 1 (y) i = a i + 1 2 d j,k=1 ρ i jk (y j y k − δ jk ). Then we obtain Recalling that B = I we see that λ i js = ρ is j + ρ jsi and then the above agrees with (20) as claimed.
We note in passing that the scheme (3) also occurs in a different setting, in the antithetic multilevel method for weak approximation in [5].
As mentioned above, the choice of p 1 for a given S 1 is far from unique if q > 1; one can add to p 1 any p ∈ P q such that y t p(y) − ∇, p(y) = 0, without changing the corresponding S 1 . An example is the scheme derived by geometric methods in [2], equation (3), which corresponds to the choice p 1 (y) i = a i + 1 2 d j,k=1 (ρ i jk + ρ jki − ρ jik )(y j y k − δ jk ): it is easily checked that this gives the same S 1 as above.
In general, given the Edgeworth-type expansion with S 1 , S 2 , · · · ) as in Sect. 4, for any choice of ( p 1 , p 2 , · · · ) such that S (0) ( p 1 , p 2 , · · · ) = (S 1 , S 2 , · · · ), we can construct an approximation for any given n by where X is an N (0, (0) ) random variable. Then, writing Z (t) for the RHS of (27), from Propositions 1 and 2, together with Theorem 1, we have W M (x(t), Z (t)) = O(t (n+2)/2 ) for any M ≥ 1. And the discussion following Proposition 2 gives (for any given M ≥ 1) a specific coupling between x(t) and Z (t) which attains this bound. Note that the last term in (27) does not give any improvement in this W M bound, but it does give the bound |Ex(t) − EZ (t)| = O(t (n+3)/2 ) which is needed to get an optimal bound for the multi-step approximation in the next section.
We remark that if n is odd then E p n+1 (X ) = S n+1 (y)ydy = 0 since S n+1 is even, so the last term on the right of (27) appears only if n is even. (kh) from which, using the bounds (29), it is straightforward to deduce that and from this that the bound (30) holds for indidual k. To get the maximal bound, note that x(kh) − u (k) can be written as is a martingale and (30) follows from the maximal martingale inequality and the uniform bound for m j .
The global condition (*) is rather restrictive and we now give an application of Theorem 2 which covers a rather more general situation. We suppose a and B are C ∞ and B has rank q on [0, T ] × U where U is an open set in R q , and consider (1) with x (0) ∈ U . We fix compact sets K and K * such that x (0) ∈ int(K ) and K ⊂ int(K * ) ⊂ U . Then for a given time-step h > 0 and n ∈ N we define a modified version of the recursion (28) as follows: as long as the RHS of (28)is in K * , we define x (k+1) by (28), but as soon as the RHS is outside K * we set x ( j) = x (k) for all j > k.
Then we assert the following: Theorem 3 With the above assumptions and notation, given M ≥ 1 and n ∈ N there exists C > 0 such that, for any h > 0, there is a probability space, equipped with a filtration F t , on which the solution x(t) to (1) is defined, along with a sequence of random vectorsx (1) ,x (2) , · · · , having the same joint distribution as (x )1) , x (2) , · · · ) as defined above, and such thatx (k) is F kh -measurable, and where τ is the escape time of x(t) from K .
Proof We can find modified coefficientsã(t, x) andB(t, x) which satisfy (*) and agree with a and B on K * . Then we apply Theorem 2 to the modified system. The vectorsx (1) ,x (2) , · · · from the modified system will have the required distribution, except that anyx (k) which is outside K * has to be modified. Since x(kh) ∈ K for k ≤ τ/h, it follows from (30) that the probability of such a modification being needed is ≤ Cδ −M h M(n+1)/2 where δ is the distance of K from the complement of K * . When no modification is needed the bound (30) applies, and (31) follows.
The idea is that one can generate random variables x (1) , x (2) , · · · using (28), and identify them withx (1) ,x (2) , · · · since the distributions are the same. In this sense the bounds (30) and (31) may be regarded as strong error bounds. In the concluding section we compare this method with standard strong and weak approximation methods and discuss the extent to which it qualifies as a strong approximation.

Concluding remarks
Here we compare the approximation scheme described above to standard strong and weak convergence schemes, and also mention some open questions.
Standard strong approximation schemes, such as described in [9], typically involve approximating the solution x(t) on [0, T ] of an SDE by approximations x (k) for k = 1, · · · , N , where the x (k) are calculated from simulated increments of the Brownian path, and possibly simulated iterated integrals of the path. One then aeeks a bound of the form E max N k=1 |x (k) − x(kh)| M ≤ Ch Mγ where γ is the order of the method. For weak approximation of order γ one is concerhed with the estimation of E f (x(T )) for suitable f and only requires bounds of the from |E f (z (N ) − E f (x(T ))| ≤ Ch γ . The bounds given in Theorems 2 and 3 then look more like standard strong bounds than weak bounds, but differ in that the random variablesx (k) are not measurable w.r.t. the probability space generated by the Brownian path, and their connection to the Brownian path is not explicit. One might desribe our method as strong approximation to a weak solution, the 'weak' referring to the fact that the solution is on a larger probability space than the Brownian path. It is natural to ask if one can get the same bounds but with x (k) measurable w.r.t. the σ -field generated by the path W . This seems quite possible but would require a different type of coupling.
To illustrate that our method is stronger than standard weak approximation we note that if f is a Lipschitz function on R q then (30) with M = 1 implies that |E f (x (N ) ) − E f (x(T ))| = O(h (n+1)/2 ), so that one can obtain weak approximation of arbitrary order for Lipschitz f , whereas standard weak approximation requires increasing differentiability of f for higher order.
Some general discussion of the interpretation of error bounds based on couplings and Vaserstein distances as a substitute for strong bounds can be found in Section 12 of [3].
Another question on our coupling construction is whether the samex (k) could be chosen for all M -the construction given above is essentially dependent on M.
A different aspect of our approximation scheme which may be worth further investigation is the choice of ( p 1 , p 2 , · · · ) corresponding to given (S 2 , S 2 , · · · ), i.e. the choice of a Cornish-Fisher expansion for a given Edgeworth expanson, which in dimension > 1 is not unique. As mentioned at the end of Section 3, W 2 bounds are easier to estimate precisely than W M in general. For example, in the context of (27), one could ask for a choice of p 1 , · · · , p n which minimises the leading order term in the expansion of W 2 (x(t), Z (t)) in powers of t. This is an essentially algebraic problem which merits exploration.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.