Skorohod and Stratonovich integration in the plane

This article gives an account on various aspects of stochastic calculus in the plane. Specifically, our aim is 3-fold: (i) Derive a pathwise change of variable formula for a path indexed by a square, satisfying some H\"older regularity conditions with a H\"older exponent greater than 1/3. (ii) Get some Skorohod change of variable formulas for a large class of Gaussian processes defined on the suqare. (iii) Compare the bidimensional integrals obtained with those two methods, computing explicit correction terms whenever possible. As a byproduct, we also give explicit forms of corrections in the respective change of variable formulas.


Introduction
Stochastic calculus for processes indexed by the plane (or higher order objects) is notoriously a cumbersome topic. In order to get an idea of this fact, let us start from the simplest situation of a smooth function x indexed by [0, 1] 2 and a regular function ϕ ∈ C 2 (R). Then some elementary computations show that for all 0 ≤ s 1 < s 2 ≤ 1 and 0 ≤ t 1 < t 2 ≤ 1, where we have set [δy] s 1 s 2 ;t 1 t 2 for the planar increment of y in the rectangle [s 1 , s 2 ] × [t 1 , t 2 ], namely [δy] s 1 s 2 ;t 1 t 2 ≡ y s 2 ;t 2 − y s 1 ;t 2 − y s 2 ;t 1 + y s 1 ;t 1 .
This simple formula already exhibits the extra term ϕ (2) (x u;v ) d u x d v x with respect to integration in R, and the mixed differential term d u x d v x is one of the main source of complications when one tries to extend (1) to more complex situations.
Moving to stochastic calculus in the plane, generalizations of (1) to a random process x obviously starts with change of variables formulas involving the Brownian sheet or martingales indexed by the plane. Relevant references include [3,13,19], and some common features of the formulas produced in these articles are the following: • Higher order derivatives of f showing up.
• Mixed differentials involving partial derivatives of x and quadratic variation type elements. • Huge number of terms in the formula due to boundary effects. This non compact form of stochastic calculus in the plane has certainly been an obstacle to its development, and we shall go back to this problem later on.
Some recent advances in generalized stochastic calculus have also paved the way to change of variables formulas in the plane beyond the martingale case. One has to distinguish two type of contributions in this direction: (a) Skorohod type formulas for the fractional Brownian sheet (abbreviated as fBs in the sequel) with Hurst parameters greater than 1/2 have been obtained in [17] thanks to a combination of differential calculus in the plane and stochastic analysis tools inspired by [1]. A subsequent generalization to Hurst parameters smaller than 1/2 is available in [18], invoking the notion of extended divergence introduced in [12]. Notice however that the extended divergence leads to a rather weak notion of integral, and might not be necessary when the Hurst parameters of the fBs are greater than 1/4. (b) The article [4] focuses on pathwise methods for stochastic calculus in the plane, and builds an analog of the rough paths theory for functions indexed by the plane. In particular, generalizations of (1) with Stratonovich type integrals are given for functions with Hölder regularity greater than 1/3. The construction is deterministic and general, and only requires the existence of a stack of iterated integrals of x called rough path, denoted by X. One can show in particular that X exists when x is a fBs.
The current article is a contribution to these recent advances on generalized stochastic calculus in the plane. Namely, we focus on 3 different problems: (i) A complete exposition of the Stratonovich type change of variables formula obtained through rough paths techniques. (ii) Generalization of [17] to a fairly general Gaussian process x. (iii) Comparison of Stratonovich and Skorohod formulas, analogously to the 1 dimensional situation handled in [10]. Before we further comment on these contributions, we now describe our main results more specifically.
1.1. Some general notation. Before we can turn to the description of our main results, we introduce some general notation concerning differential calculus in the plane. Let us mention first that we shall separate as much as possible the first and the second direction of integration, which will be respectively be denoted by direction 1 and direction 2. Thus the evaluation of a function f : [0, 1] 2k → R will be denoted by f s 1 ···s k ;t 1 ···t k . We also set d 12 x for the differential d uv x and d 1 x d 2 x for the differential d u x d v x. In fact, since the differential element d 1 x d 2 x is essential for our purposes, we further shorten it into d12x.
Another notation which will be used extensively throughout the paper is the following: we set y = ϕ(x), and for all j ≥ 1 we write y j for the function ϕ (j) (x). With those first shorthands, equation (1) for a smooth function x : [0, 1] 2 → R can be written as δy = 1 2 This kind of compact notation is of course useful when cumbersome computations come into the picture. Let us anticipate a little on the notation for planar increments which will be introduced at Section 3.1: we denote by P k,l the set of R-valued functions involving k variables in direction 1 and l variables in direction 2, satisfying some vanishing conditions on diagonals. We mostly deal with spaces of the form P 2,2 and introduce some Hölder norms. Namely, if f ∈ P 2,2 (V ), we set f γ 1 ; γ 2 = sup |f s 1 s 2 ; t 1 t 2 | |s 2 − s 1 | γ 1 |t 2 − t 1 | γ 2 ; s 1 , s 2 , t 1 , t 2 ∈ [0, 1] , and we denote by P γ 1 ,γ 2 2,2 (V ) the space of increments in P 2,2 (V ) whose · γ 1 ; γ 2 norm is finite.
1.2. Stratonovich type formula in the Young case. We assume here that x : [0, 1] 2 → R is a path such that the rectangular increments δx of x satisfy δx ∈ P γ 1 ,γ 2 2,2 with γ 1 , γ 2 > 1/2, which corresponds to the case where integration with respect to x can be handled by Young techniques in the plane. Our change of variable formula in this situation relies on the definition of 2 increments x 1;2 , x1 ;2 ∈ P γ 1 ,γ 2 2,2 defined as follows (see also Definition 4.2 for further information): where the integrals can be understood in the Young sense.
With these notations in hand, the change of variables formula can be read as: are well defined in the 2d-Young sense. Moreover: (i) Both z 1 and z 2 can be decomposed as: z 1 = y 1 x 1;2 + ρ 1 , and z 2 = y 2 x1 ;2 + ρ 2 , where ρ 1 , ρ 2 are sums of increments with double regularity (2γ 1 , 2γ 2 ) in at least one direction.
(ii) Provided x is a smooth path, the increments z 1 and z 2 are defined as Riemann-Stieljes integrals.
(iii) If x n is a sequence of smooth functions such that the related increments x n;1;2 , x n;1;2 converge respectively to x 1;2 and x1 ;2 in P γ 1 ,γ 2 2,2 , then z 1,n , z 2,n also converge respectively to z 1 and z 2 .
Observe that this theorem is not new and can be easily recovered from considerations contained in [7,15]. However, we express it here in terms which allow easy generalizations to Skorohod type integrals and to rough situations as well.
1.3. Stratonovich type formula in the rough case. Consider now a function x whose rectangular increments δx only satisfy δx ∈ P γ 1 ,γ 2 2,2 with γ 1 , γ 2 > 1/3. The definition of z 1 , z 2 as in (4) and the equivalent of formula (1) require now a huge additional effort. In particular, the correct definition relies on the introduction of a collection of iterated integrals of x (called rough path above x and denoted by X by analogy with the 1-d case) that we proceed to describe.
The reader will soon observe that the definition of X involves a whole zoology of objects which are somehow tedious to describe. In this article we shall index those objects by the directions of integration, trying to separate as much as possible direction 1 and direction (ii) We manipulate objects which are either iterated integrals or products of iterated integrals. We indicate that one starts a new integral in one specific direction and gets a product of increments by placing a new sign, and this is translated by a · in the indices of x. For instance, modifying our previous example, we define x 1·1;02 ∈ P 3,2 in the following way for a smooth function x: Notice that those breaks in integration can occur at different steps in each direction 1 or 2. The resulting overlapping integrals will be an important source of technical troubles in the remainder of the paper. With these preliminary considerations in mind, our assumptions on the function x are of the following form: The function x is such that δx ∈ P γ 1 ,γ 2 2,2 with γ 1 , γ 2 > 1/3. Moreover, the following rough path X can be constructed out of x: Increment Interpretation Regularity Increment Interpretation Regularity 1 2 d 12 xd 12 x (2γ 1 , 2γ 2 ) x 11;22 1 2 d 12 xd12x (2γ 1 , 2γ 2 ) x1 1;22 1 2 d12xd 12 x (2γ 1 , 2γ 2 ) x11 ;22 In the table above, all increments belong to P 2,2 , so that a regularity (α, β) means that the increment lyes into P α,β 2,2 . Furthermore, the stack X is a geometric rough path, insofar as there exists a regularization x n of x such that lim n→∞ x − x n γ 1 ,γ 2 = 0 and such that all the integrals in X n , constructed out of x n in the Lebesgue-Stieljes sense, converge with respect to their natural respective norms in P γ 1 ,γ 2 2,2 , P 2γ 1 ,γ 2 2,2 , P γ 1 ,2γ 2 2,2 or P 2γ 1 ,2γ 2 2,2 . Note that the natural Hölder norm of a rough path is denoted by N in the sequel. Remark 1.3. As we shall see at Section 5, Hypothesis 1.2 is not completely sufficient in order to settle a satisfying integration theory with respect to x. In fact the rough path X should also include higher order increments like x 11·1;022 or x 1·11;22·2 (and other extra terms). We have only stated Hypothesis 1.2 here in order to keep our exposition into some reasonable bounds. Now we can state our Stratonovich integration theorem in the rough case, which mimics Theorem 1.1: with γ 1 , γ 2 > 1/3 and assume the further rough path Hypothesis 1.2. Consider a function ϕ ∈ C 8 b (R). Then the increments z 1 and z 2 given by (4) are well defined as continuous functions of the rough path X. Moreover: (i) The increment z 1 can be decomposed as: z 1 = y 1 x 1;2 + y 2 x 11;02 + y 2 x 01;22 + y 2 x 11;22 + y 3 x1 1;22 + ρ 1 , and the increment z 2 admits a decomposition of the form z 2 = y 2 x1 ;2 + y 3 x 11;02 + y 3 x 01;22 + y 3 x 11;22 + y 4 x11 ;22 + ρ 2 , where ρ 1 , ρ 2 are sums of increments with triple regularity (3γ 1 , 3γ 2 ) in at least one direction.
(ii) If x n is a sequence of smooth functions such that the related rough path X n converges to X, then z 1,n , z 2,n (defined in the Lebesgue-Stieljes sense) also converge respectively to z 1 and z 2 .
(iii) The change of variables formula (3) still holds true when integrals are understood in the rough path sense.
Obviously, Theorem 1.4 would be of little interest if we could not apply it to processes of interest. To this regard, our guiding example will be the fractional Brownian sheet (fBs in the sequel). Let us recall that this is a centered Gaussian process x defined on [0, 1] 2 , with a covariance function R s 1 s 2 ; where the Hurst parameters γ 1 , γ 2 lye into (0, 1). Many possible representations are available for the fBs, among which we will appeal to the so-called harmonizable representation (see relation (86) below for further details). This allows a natural approximation of x by a sequence of smooth processes x n thanks to a cutoff in frequency, and we recall the following convergence result established in [4]: Let x a fBs with Hurst parameters γ j > 1/3, for j = 1, 2. Define the regularization x n of x given by a frequency cutoff on B(0, n) in the harmonizable representation of x. Then: (i) The family of iterated integral X n defined in (1.2) associated to x n fulfills the relation lim n,m→∞ E[N p (X m − X n )] = 0 for all p ≥ 1, where the norm N is alluded to at Hypothesis 1.2. The limit object X is called rough sheet associated to x.
(ii) Theorem 1.4 applies to the fBs x.
As the reader might imagine, Theorem 1.4 can also be applied to a wide range of Gaussian and non Gaussian processes. We focus here on fBs for sake of simplicity.
1.4. Skorohod integration. One of the main issue alluded to in this article is a comparison between Stratonovich and Skorohod type change of variable formulas when x is a Gaussian process exhibiting some Hölder regularity in the plane. Towards this aim, our global strategy is to use our Theorems 1.1 and 1.4 and compute corrections between Stratonovich and Skorohod type integrals.
We first focus on the Young case, assuming the same regularity conditions as in Section 1.2. We are then able to handle the case of a fairly general centered Gaussian process x whose covariance function R satisfies a factorization property of the form for two covariance functions R 1 , R 2 on [0, 1] and such that R 1 , R 2 ∈ C 1-var ([0, 1] 2 ) (which ensures that x is (γ 1 , γ 2 )-Hölder continuous with γ 1 , γ 2 > 1/2). Notice in particular that the fBs covariance function (8) satisfies condition (9).
The standard growth assumptions on f in order to get a Skorohod formula for f (x) should also be met. They will feature prominently in the sequel, and we proceed to recall them now: Definition 1.6. Let k ∈ N, we will say that a function f ∈ C k (R) satisfies the growth condition (GC) if there exist positive constants c and λ such that , and max l=0,...,k With these notations in hand, and denoting quite informally the Skorohod differentials by d ⋄ (see Section 6.1.1 for further explanations), we can summarize our results in the following: Theorem 1.7. Assume x is a centered Gaussian process on [0, 1] 2 with a covariance function satisfying (9). Consider a function ϕ ∈ C 4 (R) satisfying condition (GC). Then the increments are well defined in the Skorohod sense of Malliavin calculus. Moreover: (i) Some Riemann convergences hold true: if π 1 n and π 2 n are 2 partitions of [s 1 , s 2 ] × [t 1 , t 2 ] whose mesh goes to 0 as n → ∞, then where ⋄ stands for the Wick product in the left hand side of the relations above, and where the convergence holds in both a.s and L 2 (Ω) sense.
Finally, let us move to the Skorohod change of variable in the rough situation. For simplicity of exposition, we have restricted our analysis to the fractional Brownian sheet, mainly because our computations heavily hinges on the explicit regular approximation sequence x n given by the harmonizable representation of fBs (similarly to the construction of the rough path above x). The Skorohod change of variable (consistent with the formulas obtained in [18]) and Skorohod-Stratonovich comparison we obtain in this case are summarized as follows: Theorem 1.8. Assume x is a fractional Brownian sheet on [0, 1] 2 , with γ j > 1/3 for j = 1, 2. Then the increments z 1,⋄ , z 2,⋄ of equation (4) are well defined in the Skorohod sense of Malliavin calculus. Moreover: (i) Both z 1,⋄ and z 2,⋄ can be seen as respective limits of z n,1,⋄ and z n,2,⋄ , computed as in Theorem 1.7 for the regularized process x n .
1.5. Further comments. As the reader might have noticed, our paper gives a rather complete picture of pathwise Stratonovich and Itô-Skorohod integration for processes indexed by the plane. In the case of pathwise integration, we take up the methodology introduced in [4], based on a highly nontrivial extension of the rough path theory, and compute explicitly all the terms involved in our Stratonovich expansion. We wished the exposition to be as clear and self-contained as possible, which explains the length of Sections 4 and 5.
In order to put our strategy for the Itô case into perspective, notice that 2 types of methodologies are usually available for changes of variables in case of a Gaussian process x: (a) Define a divergence type operator δ ⋄ for x and proceed by integration by parts on expressions like E[δf (x) G], where G is a smooth functional of x. This is the strategy invoked e.g. in [1,10].
(b) Base the calculations on the pathwise change of variables formula of type (1). This formula is generally related to some converging Riemann sums like in Theorem 1.1, and one can compute corrections between Wick and ordinary products in relation (6). This is the method implicitly adopted in [17] and we also resort to this second strategy here, which allows to derive our Skorohod formula and its comparison with the Stratonovich formula at the same time.
Unfortunately, the Wick corrections strategy does not work for the rough case, even in the explicit situation of a fractional Brownian sheet. This mainly stems from the fact that convenient Riemann sums related to formula (1) are not available (so far) in the case of Theorem 1.4. This drawback led us to change our strategy again, and proceed by regularization. Indeed, as mentioned before, one can come up with an explicit regular approximation x n of x. For this regularization, we can apply Theorem 1.7 and get some Itô-Stratonovich corrections. Invoking the fact that δ ⋄ is a closable operator, we can then take limits in our operations as n → ∞. This allows to compare the changes of variables formulas (1) and (15), but the interpretation in terms of Riemann-Wick sums is obviously lost in this case. Notice that an approximation procedure (expressed in terms of the extended divergence operator) is also at the heart of [18] for irregular fBs.
Finally, let us say a few words about possible extensions of our work: • Generalizations of Skorohod's change of variable to a Gaussian process without the factorization hypothesis (9) on the covariance function of x are certainly possible. However, at a technical level, one should be aware of the fact that the analysis of mixed terms like x u;v would require tools of Young integration in dimension 4. These techniques have been used e.g in [6], and the elaboration we need would certainly be cumbersome. We have thus sticked to the factorized case for R for sake of readability.
• As mentioned before, our strategy for the Skorohod formula in the rough case relies heavily on a suitable regularization of x. Instead of treating the explicit fBs example, we could have stated some general approximation assumptions satisfied in the fBs case. Once again, we have chosen to specialize our study here for sake of clarity. The general case might be handled in a subsequent paper, and we also hope to design a strategy based on Riemann-Wick sums in the next future.
Here is how our article is structured: We recall some basic notation of algebraic integration in dimension 1 at Section 2, and extend it to integration in the plane at Section 3. The Stratonovich change of variable formula is handled at Section 4 for the Young case and at Section 5 in the rough situation. We then move to Skorohod type formulas at Sections 6 and 7, respectively for the regular and rough cases.

Algebraic integration in dimension 1
We recall here the minimal amount of notation concerning algebraic integration theory in R, in order to prepare the ground for further developments in the plane. We refer to [8,9] for a more detailed introduction.
2.1. Increments. The extended pathwise integration we will deal with is based on the notion of increments, together with an elementary operator δ acting on them. The algebraic structure they generate is described in [8,9], but here we present directly the definitions of interest for us, for sake of conciseness. First of all, for a vector space V and an integer k ≥ 1 we denote by C k (V ) the set of functions g : [0, 1] k → V such that g t 1 ···t k = 0 whenever t i = t i+1 for some i ≤ k − 1. Such a function will be called a (k − 1)-increment, and we set C * (V ) = ∪ k≥1 C k (V ). We can now define the announced elementary operator δ on C k (V ): wheret i means that this particular argument is omitted. A fundamental property of δ, which is easily verified, is that δδ = 0, where δδ is considered as an operator from C k (V ) to C k+2 (V ).
Some simple examples of actions of δ, which will be the ones we will really use throughout the paper, are obtained by letting g ∈ C 1 and h ∈ C 2 . Then, for any s, u, t ∈ [0, 1], we have Furthermore, it is easily checked that ZC k+1 (V ) = BC k (V ) for any k ≥ 1. In particular, the following basic property holds: Lemma 2.1. Let k ≥ 1 and h ∈ ZC k+1 (V ). Then there exists a (non unique) f ∈ C k (V ) such that h = δf .
Lemma 2.1 can be rephrased as follows: any element h ∈ C 2 (V ) such that δh = 0 can be written as h = δf for some (non unique) f ∈ C 1 (V ). Thus we get a heuristic interpretation of δ| C 2 (V ) : it measures how much a given 1-increment is far from being an exact increment of a function, i.e., a finite difference.
Notice that our future discussions will mainly rely on k-increments with k ≤ 2, for which we will make some analytical assumptions. Namely, we measure the size of these increments by Hölder norms defined in the following way: for f ∈ C 2 (V ) let Obviously, the usual Hölder spaces C µ 1 (V ) will be determined in the following way: for a continuous function g ∈ C 1 (V ), we simply set and we will say that g ∈ C µ 1 (V ) iff g µ is finite. Notice that · µ is only a semi-norm on where the last infimum is taken over all sequences {h i ∈ C 3 (V )} such that h = i h i and for all choices of the numbers ρ i ∈ (0, z). Then · µ is easily seen to be a norm on C 3 (V ), and we set , and notice that the same kind of norms can be considered on the spaces ZC 3 (V ), leading to the definition of some spaces ZC µ 3 (V ) and ZC 1+ 3 (V ). With these notations in mind the following proposition is a basic result, which belongs to the core of our approach to pathwise integration. Its proof may be found in a simple form in [9].
Proposition 2.2 (The Λ-map). There exists a unique linear map Λ : and Λδ = Id C 1+ 2 (V ) . In other words, for any h ∈ C 1+ 3 (V ) such that δh = 0 there exists a unique g = Λ(h) ∈ C 1+ 2 (V ) such that δg = h. Furthermore, for any µ > 1, the map Λ is continuous from ZC µ 3 (V ) to C µ 2 (V ) and we have Let us mention at this point a first link between the structures we have introduced so far and the problem of integration of irregular functions. Corollary 2.3. For any 1-increment g ∈ C 2 (V ) such that δg ∈ C 1+ 3 , set δf = (Id − Λδ)g. Then where the limit is over any partition Π st = {t 0 = s, . . . , t n = t} of [s, t], whose mesh tends to zero. Thus, the 1-increment δf is the indefinite integral of the 1-increment g.

Products of increments.
For notational sake, let us specialize now to the case V = R, and just write C γ k for C γ k (R). The usual product of two increments considered on (C * , δ) is obtained by gluing one variable in each increment (see e.g [8,9]): Definition 2.4. For g ∈ C n and h ∈ C m , we denote by gh the element of C n+m−1 defined by (gh) t 1 ,...,t m+n−1 = g t 1 ,...,tn h tn,...,t m+n−1 , t 1 , . . . , t m+n−1 ∈ [0, 1].
However, another product (defined without gluing of variables) turns out to be useful for further computations in the plane. This product is called splitting and is defined below: Definition 2.5. For g ∈ C n and h ∈ C m , we denote by S(g, h) the element of C n ⊗ C m defined by [S(g, h)] t 1 ,...,t m+n = g t 1 ,...,tn h t n+1 ,...,t m+n , t 1 , . . . , t m+n ∈ [0, 1].
(23) Notice that S(g, h) can also be considered as an increment in C n+m , except that it is not required to vanish when t n = t n+1 .
The splitting operation will have to be inverted at some point. The inverse operation is called gluing: Definition 2.6. For n, m ≥ 1, we say that f is an element of M n,m if it can be written as a finite linear combination of the following type: For f of the form (24), we then define G(f ) ∈ C n+m−1 as: [G(f )] t 1 ,...,t n+m−1 = j∈J α j g j t 1 ,...,tn h j tn,...,t m+n−1 .
With these definitions in mind, let us remark that if f ∈ M n,m is simply of the form f t 1 ,...,t n+m = g t 1 ,...,tn h t n+1 ,...,t m+n with g ∈ C n and h ∈ C m , then G(f ) = g h.
We now recall some elementary properties concerning products of increments: Proposition 2.7. The following differentiation rules hold true: (1) Let g ∈ C 1 and h ∈ C 1 . Then gh ∈ C 1 and δ(gh) = δg h + g δh.
Proof. We will just prove (25), the other relations being just as simple. If g, h ∈ C 1 , then which proves our claim.

Iterated integrals as increments.
Iterated integrals of smooth functions on [0, 1] are obviously particular cases of elements of C 2 , which will be of interest for us. A typical example of this kind of object is given as follows: consider f j ∈ C ∞ 1 for j = 1, . . . , n and 0 ≤ s 1 < s 2 ≤ 1. For n ≥ 1, we denote by S n (s 1 , s 2 ) the simplex S n (s 1 , s 2 ) = {(σ 1 , . . . , σ n ) ∈ [0, 1] n ; s 1 < σ 1 < · · · < σ n < s 2 } and we set We now introduce some notation for iterated integrals which is much too complicated for integration in dimension 1, but turns out to be useful for integration in the plane. Indeed, we can alternatively denote the increment h 1,...,n defined at (29) by where the integration on the n-dimensional simplex is implicit in both cases. We shall also need a small variant of these conventions: we set where all the products are understood as products of increments as in Definition 2.4. With these conventions in mind, the following relations between multiple integrals and the operator δ will also be useful. The reader is sent to [9] for its elementary proof.

Algebraic integration in the plane
This section is devoted to recall the elements of algebraic integration necessary to define an integral of the form [0,1] 2 f (x) dx for a Hölder function x in the plane with Hölder exponent greater than 1/3. This requires a tensorization of the algebraic structures defined in the previous section, plus some extra tools that we proceed to introduce.
3.1. Planar increments. We consider here increments of a variable s (also called direction 1) and a variable t (also called direction 2), with (s, t) ∈ [0, 1] 2 . For a vector space V , we set In the particular case V = R, we simply set P k,l (R) ≡ P k,l . Some partial difference operators δ 1 and δ 2 with respect to the first and second direction can be defined as in the previous section. Namely, for f ∈ P k,l (V ) we set and we define δ 2 similarly. The planar increment δ is then obtained as δ = δ 1 δ 2 . Notice that for f ∈ P 1,1 we have which is the usual rectangular increment of a function f defined on [0, 1] 2 and is consistent with formula (2). Let us label the following notation for further use: Notation 3.1. For j = 1, 2, we set Z j P k,l = P k,l ∩ ker(δ j ) and B j P k,l = P k,l ∩ Im(δ j ). We also write ZP k,l for P k,l ∩ ker(δ) and BP k,l for P k,l ∩ Im(δ).
As in the 1-d case, the Hölder regularity of planar increments is an essential feature of our generalized integration theory. On P 2,2 (V ) and P 3,3 (V ), it is measured by a tensorization of the Hölder norms defined at (18) and (20). Namely, if f ∈ P 2,2 (V ), we set and we denote by P γ 1 ,γ 2 2,2 (V ) the space of increments in P 2,2 (V ) whose · γ 1 ; γ 2 norm is finite. Along the same lines, we say that h ∈ P γ 1 ,γ 2 3 Similar norms, omitted here for sake of conciseness, can be defined on P 2,3 (V ) and P 3,2 (V ).
For Hölder continuous increments with regularity greater than 1, one gets the following inversion properties, which are a direct consequence of the one dimensional Proposition 2.2: We do not include the proof of this proposition for sake of conciseness. Let us just mention that (as the reader might imagine) we have Λ = Λ 1 Λ 2 . It should also be observed that some 2-dimensional Riemann sums are related to the sewing map Λ, echoing Corollary 2.3: Proposition 3.3. Let g ∈ P 2,2 satisfying the following assumptions: where * denotes any kind of Hölder regularity. Then there exists f ∈ P 1,1 such that where π designates a family of rectangular partitions of [s 1 , s 2 ] × [t 1 , t 2 ] whose mesh goes to 0.
Notice that inverse operations of splittings can also be defined for planar increments. For sake of conciseness, we let the patient reader generalize Definition 2.6 in order to get the definition of gluing G 1 , G 2 and G for planar increments.
We close this section by introducing a last product of increments which is labeled for further computations.
Definition 3.6. Let g ∈ P 2,1 and h ∈ P 1,2 . Then g • h is the increment in P 2,2 defined by 3.3. Iterated integrals as increments in the plane. The relationship between iterated integrals and increments in the plane is crucial for us. Generally speaking, an iterated integral is given as follows: consider f j ∈ P ∞ 1,1 for j = 1, . . . , n and (s 1 , s 2 ), (t 1 , t 2 ) ∈ S 2 , where we recall relation (28) defining simplexes. Then we set where we recall from Section 1.1 that d 12 f j σ;τ stands for ∂ 2 στ f j σ;τ . Expression (32) is obviously cumbersome, and it could in particular become clearer by separating the s, σ from the t, τ variables. This is where the conventions introduced in equation (30) turn out to be useful. Namely, one can simply tensorize (30) in order to write the increment h 1,...,n defined at (32) as and notice that we will mainly use the second convention throughout the paper. This notation proves to be particularly convenient when one is faced with partial integrations (as introduced in (31)) in both directions 1 and 2. In order to illustrate this point, let us consider the simple example Let us now describe the algorithm which allows to go from expression (34) to an integral like (32). It can be summarized as follows: • For direction 1, count the number of iterated integrals starting from the left hand side (in our example this number is 2). Then interpret these integrals as integrals on the simplex in direction 1 and write the 1-variables. • Do the same for direction 2. In our example, there is only one integral in this direction, so that variable 2 is frozen in the first differential d 1 f 1 . Applying this algorithm, the reader can easily check that g defined at (34) can be written as For sake of conciseness, we omit generalizations of this simple example.

Planar Young integration
Before going into the computational details of Young integration, let us describe the general strategy we shall follow in order to obtain our Itô-Stratonovich type change of variable formulae in case of non smooth functions x. Indeed, we start from a smooth approximation x n to our path x and we introduce a useful notation for the remainder of the computations: Notation 4.1. We shall drop the index n of approximations in x n , which means that x will stand for a generic smooth path defined on [0, 1] 2 . For a smooth function ϕ : R → R, we also write y for the path ϕ(x) and for all j ≥ 1 we set y j = ϕ (j) (x).
With these notations in hand, for a smooth sheet x and f ∈ C 2 b it is well known that formula (1) holds true. Recall that we have written this relation under the following form, compatible with our convention (33): We shall see that this formula still holds true in the limit for x, except that the integrals involved in the right hand side of (35) have to be interpreted in a sense which goes beyond the Riemann-Stieltjes case. Our main task will thus be to obtain a definition of 1 2 y 1 d 12 x and 1 2 y 2 d12x involving iterated integrals of x and increments of y (or y j for j ≥ 1) only. Though this task might overlap with some aspects of [4], we present it here because it is short enough and allows us to introduce part of our formalism.
Let us introduce what will be later interpreted as the first order elements of the planar rough path above x: with γ 1 , γ 2 > 1/2. We set Notice that for smooth functions we also have We are now ready to express the integrals in (35) in terms of the planar sewing map (Λ i ) i=1,2 and Λ = Λ 1 Λ 2 and the increments introduced in Notation 4.2.

4.1.
Change of variables formula. The following theorem gives the analog of relation (35) in the Young case, and is a way to recast Theorem 1.1. Notice that some extensions of these results are contained in [4].
with γ 1 , γ 2 > 1/2 and ϕ ∈ C 4 (R). With our notations 4.2 and 4.1 in mind, define two increments z 1 , z 2 as Then items (i)-(v) of Theorem 1.1 hold true. Furthermore, for j = 1, 2 the following bound is satisfied: . Proof. We first introduce a formalism which will feature prominently in the rough case treated in the next section.
indices are cumbersome, we shall now show how to translate the above computations with the formalism of Section 3.3: we simply write ≡ y 1 x 1;2 + a 11;02 + a 01;22 + a 11;22 , which is obviously a shorter expression than (37). From now on, we shall carry on our computations with this simplified formalism.

4.2.
Riemann sums decompositions. This section is meant as a preparation for Skorohod type computations. Indeed, change of variables in the Skorohod setting involve some mixed integrals with dx dR terms, for which a suitable representation is required. It will also be convenient for us to express the integral 1 2 y 2 d 1 xd 2 x in different ways, so that we first recall a proposition borrowed from [4]: for ϕ ∈ C 2 (R) and z 2,y ≡ 1 2 y d 1 xd 2 x, understood in the Young sense. Then the following series of identities hold true: where we recall that the notation • has been introduced at Definition 3.6.
The following proposition gives different ways to express the increment z 2,y as limit of Riemann sums. Proposition 4.5. Let 0 < s 1 < s 2 < 1, 0 < t 1 < t 2 < 1 and denote by π 1 = (s i ) i and π 2 = (t j ) some partitions of the intervals [s 1 , s 2 ] and [t 1 , t 2 ] respectively. Then under the assumptions of Proposition 4.4 we have that 12 x st can be written as limit of Riemann sums of the form lim |π 1 |,|π 2 |→0 i,j y s i ;t j δx s i s i+1 ;t j t j+1 , and recalling that z 2,y ≡ 1 2 y d 1 xd 2 x we also have: Finally we shall need an extension of the last three propositions to integrals with mixed driving noises: Then we also have Moreover, taking up the notations of Proposition 4.5 we have the following Reimann-sum representation of our integrals: Then by a simple computation we have that This means that a 1 satisfies the assumptions of Proposition 3.2, and the same is readily checked for a 2 and a 3 . Thus the increments We now identify the increments [Id − Λ 1 δ 1 ][Id − Λ 2 δ 2 ](a j ) by analyzing their Riemann sums. Indeed, a straightforward application of Proposition 3.3 yields the following limits: We now prove that the 3 increments coincide by showing that the differences between Riemann sums vanish when the mesh of the partitions go to zero. Indeed, if we call π 2 the partition in direction 2, observe for instance that The same relation holds true for a 2 − a 3 by symmetry, which ends our proof.

Pathwise Stratonovich formula in the rough case
In this section, we consider a path x ∈ P γ 1 ,γ 2 1,1 with γ 1 , γ 2 > 1/3 and we wish to establish the change of variables formula for y = ϕ(x) (with ϕ ∈ C 8 b (R)) announced in Theorem 1.4. In such a general context, the integration theory with respect to x relies on the existence of a rough path X sitting above x. Now recall that the first elements of X have been introduced at Hypothesis 1.2. However, as mentioned in Remark 1.3, the rough path still has to be completed and we proceed to its description here.
Let us first introduce another indexing convention for the elements of the rough path X, similarly to what is done at Section 1.3: (iii) (Following (ii) at Section 1.3) We shall see that some overlapping integrals in directions regularities related to each increment. As in Hypothesis 1.2, the iterated integrals above are Table 1. Further elements of X

Increment
Interpretation Increment Interpretation assumed to be limits along approximations of x by smooth functions. Furthermore, the rough path X should also contain all the elements x obtained by symmetrizing the increments of Table 1 with respect to 1 ↔ 2, as well as those for which we change the last indices 1; 2 bŷ 1;2. In total, we have to assume the existence of 26 additional increments.
With the additional notation introduced above, we are now ready to give the expression of z 1 in terms of the rough path X.
Proof of Theorem 1.4. Consider a smooth function x and y = ϕ(x). According to our notational conventions of Section 3.3, we wish to express 1 2 y 1 d 12 x in terms of elementary increments of y, plus some elements of the rough path X. As in the Young case, we thus start from expression (38), but now the terms a 11;02 , a 01;22 and a 11;22 need a further decomposition. We proceed to their analysis.
Step 1: First boundary integrals. Consider the term a 11;02 in expression (38). Two building blocks of our methodology can be observed on the analysis of this term in an elementary way: (i) As noticed in the proof of Proposition 4.3, our definition a 11;02 = 1 d 1 y 1 2 d 12 x easily yields the identity δ 1 a 11;02 = δ 1 y 1 δx.
Indeed, a possible way to understand this relation is to write: where we have invoked Proposition 2.8 for the first variable only (ii) If we only consider the integral 1 d 1 y 1 within a 11;02 , we can write where we recall that we have set y 1 = f ′ (x) and y 2 = f ′′ (x) in Notation 4.
Step 2: First decomposition of the double integral. Let us handle now the term a 11;22 , starting by the equivalent of ingredient (ii) above. We thus write an expansion for the term 1 2 d 12 y, which can be handled thanks to relation (35) as We proceed now to the analysis of the term b 11;22 . To this aim, go back to relation (38) in order to write and b 111;222 = 1 2 d 12 y 2 d 12 x d 12 x. We shall treat those 3 terms separately.
Step 3: Analysis of the boundary integrals. The term b 111;022 is a third order iterated integral in direction 1. One can thus try to apply ingredient (i) above, namely compute δ 1 b 111;022 in order to observe if we get an increment in P 3γ 1 , * 3,2 . Unfortunately, this brings some extra complications due to overlapping integrations. This is why we have postponed these computations to Lemma 5.2 below, which asserts that δ 1 b 111;022 ∈ P 3γ 1 ,γ 2
Step 5: Analysis of b1 1;22 . We now go back to the term b1 1;22 introduced at equation (46). Its analysis is very similar to the one we developed for b 11;22 , so that we only sketch the main differences.
First of all, in order to get an equivalent to relation (38), we resort to the differential element d12x. This yields: and we plug this relation into (46) in order to get . Similarly to (50), we thus end up with a relation of the form where ρ y 2 ,1 ∈ P 3γ 1 ,γ 2
Finally, we propagate this identity back into (46), (45), (44) and eventually (38), which yields our identity (7). We leave the details of these last computations to the patient reader.
Step 7: Conclusion. We have just obtained expression (7) for z 1 . As far as z 2 is concerned, we start by replacing relation (38) by: The expansion then proceeds exactly like in the expansion of z 1 , and details are omitted for sake of conciseness. The claims (ii) and (iii) of Theorem 1.4 follow by a limiting procedure.
Now, under Hypothesis 5.1 the term δ 1 y 2 x 11;22 above is easily seen to lye into P 3γ 1 ,2γ 2 3,2 . However, the term h is more troublesome. In order to see this additional problem more clearly, let us specify h a little: for s 1 , s 2 , s 3 ∈ S 3 and t 1 , t 2 ∈ S 2 it can be written as and the regularity of this last increment is far from being obvious to determinate. Furthermore, a closer look at h leads to the following conclusion: the problem in the regularity analysis stems from the overlap in integrations for directions 1 and 2, apparent in the definition of h. In order to fix this new problem (which can obviously only occur when integrals of order 3 or more are showing up), we introduce a third ingredient of our methodology, based on the splitting operation of Definition 3.5: (iii) The splitting operation can be applied to overlapping terms like h 11·1;022 , in order to separate the term 1 d 1 y 2 2 d 12 x from the term 1 d 12 x. This simple trick allows to solve the intricate structure of integration in both directions 1 and 2. Unfortunately, the remaining terms have to be defined separately now, and are not of order 3 anymore. This means that another step of expansion might be necessary for their proper definition. In our running example, this additional step concerns the term 1 d 1 y 2 2 d 12 x in (53). Ingredient (iii) applied to our increment h 11·1;022 yields the following developments: we have and let us add the following comments: (a) In direction 2, the quantity S 1 (h 11·1;022 ) can be simply expressed in terms of iterated integrals of the driving process x and elementary increments of y 1 . (b) The exact expression of S 1 (h 11·1;022 ) can be obtained from (54) by splitting the variables s 2 , s 3 into 3 variables s 2 , s 3 , s 4 : We now set h 11;02 = 1 d 1 y 2 2 d 12 x, which is the increment we shall further expand in (55). This can be performed exactly as for a 11;02 in Step 2, just replacing y by y 2 , and enables the following expression which should be compared to equality (44): h 11;02 = y 3 x 11;02 + Λ 1 r y 2 ,1 δx + δ 1 y 3 x 11;02 , where we recall that x 11;02 = 1 d 1 x 2 d 12 x. We plug this relation into (55) and we obtain In the last expression, h 11·1;022;♯ is defined as follows: we have We can now go back to an expression for h 11·1;022 itself by applying the inverse map G 1 of S 1 . This gives (53), we end up with

Plugging this relation back in equation
which is now clearly an element of P 3γ 1 ,γ 2
Proof. Those two terms contain overlapping integrations in direction 1 and 2 again, and we thus proceed to their dissection through ingredient (iii). In the case of b 11·1;2·22 we will use an additional expansion in direction 1 only, and thus write In order to expand c 11;2·2 one step further in direction 1, observe that Plugging this identity into the definition of c 11;2·2 we get a decomposition as We now apply ingredients (i) and (ii) to analyze those terms. First, one can compute δ 1 c 11;2·2 as Next expand y 3 in direction 1 in both k 1 and k 2 , which yields k 1 = δ 2 y 3 x 11;02 + ρ 1 , and k 2 = y 3 x 11;2·2 + ρ 2 , where ρ 1 and ρ 2 are two remainder terms morally lying in P 3γ 1 ,2γ 2

Skorohod's calculus in the Young case
This section is devoted to relate the Young type integration theory introduced at Section 4 and the Skorohod integral in the plane handled in [17]. Specifically, we shall first generalize the Skorohod change of variables formula given in [17] for a fractional Brownian sheet with Hurst parameter greater than 1/2 to a fairly general Gaussian process. We shall then compare this formula with Theorem 1.1 item (v).
Before going on with our computations, let us label some notations for further use: 6.1. Malliavin calculus framework. We consider in this section a centered Gaussian process {x s;t ; (s, t) ∈ [0, 1] 2 } defined on a complete probability space (Ω, F , P), with covariance function E[x s 1 ;t 1 x s 2 ;t 2 ] = R s 1 s 2 ;t 1 t 2 . We now briefly define the basic elements of Malliavin calculus with respect to x and then specify a little the setting under which we shall work. Starting from the space H, a Malliavin calculus with respect to x can now be developped in the usual way (see [10,14] for further details). Namely, we first define a set of smooth functionals of x by S := {f (I 1 (ψ 1 ), . . . , I 1 (ψ n )); n ∈ N, f ∈ C ∞ b (R n ), ψ 1 , . . . , ψ n ∈ H} and for F = f (I 1 (ψ 1 ), . . . , I 1 (ψ n )) ∈ S we define ψ 1 ), . . . , I 1 (ψ n )) ψ i .
Then D is a closable operator from L p (Ω) into L p (Ω, H). Therefore we can extend D to the closure of smooth functionals under the norm The iteration of the operator D is defined in such a way that for a smooth random variable F ∈ S the iterated derivative D k F is a random variable with values in H ⊗k . The domain D k,p of D k is the completion of the family of smooth random variables F ∈ S with respect to the semi-norm : Similarly, for a given Hilbert space V we can define the space D k,p (V ) of V -valued random variables, and D ∞ (V ) = ∩ k,p≥1 D k,p . Consider now the adjoint δ ⋄ of D. The domain of this operator is defined as the set of u ∈ L 2 (Ω, H) such that E[| DF, u H |] F 1,2 , and for this kind of process δ ⋄ (u) (called Skorohod integral of u) is the unique element of L 2 (Ω) such that The following divergence type property of δ ⋄ will be useful in the sequel: and we also recall the following compatibility of δ ⋄ with limiting procedures: Lemma 6.2. let u n be a sequence of elements in Dom(δ ⋄ ), which converges to u in L 2 (Ω, H). We further assume that δ ⋄ (u n ) converges in L 2 (Ω) to some random variable F ∈ L 2 (Ω). Then u ∈ Dom(δ ⋄ ) and δ ⋄ (u) = F .

6.1.2.
Wick products. Some of our results below will be expressed in terms of Rieman-Wick sums. We give a brief account on these objects, mainly borrowed from [10,11]. Among functionals F of x such that F ∈ D ∞ , the set of multiple integrals plays a special role. In order to introduce it in the context of a general process x indexed by the plane, consider an orthonormal basis {e n ; n ≥ 1} of H and let⊗ denote the symmetric tensor product. Then is an element of H⊗ n satisfying the relation: Moreover, H⊗ n is the completion of the set of elements like (62) with respect to the norm (63).
For an element f n ∈ H⊗ n , the multiple Itô integral of order n is well-defined. First, any element of the form given by (62) can be rewritten as where the j 1 , . . . , j m are different and k 1 + · · · + k m = n. Then, if f n ∈ H⊗ n is given under the form (64), define its multiple integral as: I n (f n ) = finite f j 1 ,··· ,jm H k 1 (I 1 (e j 1 )) · · · H km (I 1 (e jm )), where H k denotes the k-th normalized Hermite polynomial given by It holds that the multiple integrals of different order are orthogonal and that H⊗ n . This last isometric property allows to extend the multiple integral for a general f n ∈ H⊗ n by L 2 (Ω) convergence. Finally, one can define the integral of f n ∈ H ⊗n by putting I n (f n ) := I n (f n ), wheref n ∈ H⊗ n denotes the symmetrized version of f n . Moreover, the chaos expansion theorem states that any square integrable random variable F ∈ L 2 (Ω, G, P), where G is the σ-field generated by x, can be written as With these notations in mind, one way to introduce Wick products on a Wiener space is to impose the relation for any f n ∈ H⊗ n and g m ∈ H⊗ m , where the multiple integrals I n (f n ) and I m (g m ) are defined by (65). If F = N 1 n=1 I n (f n ) and G = N 2 m=1 I m (g m ), we define F ⋄ G by I n+m (f n⊗ g m ).
By a limit argument, we can then extend the Wick product to more general random variables (see [11] for further details). In this paper, we will take the limits in the L 2 (Ω) topology.
Some corrections between ordinary and Wick products will be computed below. A simple example occurs for products of f (x) by a Gaussian increment. Indeed, for a smooth function f and g 1 , g 2 ∈ H, it is shown in [11] that f (I 1 (g 1 )) ⋄ I 1 (g 2 ) = f (I 1 (g 1 )) I 1 (g 2 ) − f ′ (I 1 (g 1 )) g 1 , g 2 H . (68) We now state a result which is proven in [10, Proposition 4.1].

Further assumptions and preliminary results.
In order to simplify our computations, let us introduce some additional assumptions on the covariance R: Hypothesis 6.4. The covariance R of our centered Gaussian process x belongs to the space C 1-var ([0, 1] 4 ), and satisfies a factorization property of the form 1]. In addition, setting R i a = R i aa for a ∈ [0, 1] and i = 1, 2, we assume that a → R i a is differentiable and we suppose that 1]). The first consequence of our Hypothesis 6.4 is that the regularity of x corresponds to the Young type regularity of Section 4. Indeed, it is readily checked that relation (69) yields Since x is Gaussian, an easy application of Kolmogorov's criterion ensures that for arbitrarily small ǫ 1 , ǫ 2 > 0. This enables us to appeal to Young's integration theory in order to define integrals of the form 1 2 ϕ(x) d 12 x.
Next we also need an integral semi-norm dominating Hölder's norms in the plane. This is given by the following Garsia type result: Lemma 6.6. Let p > 1, θ 1 , θ 2 > 0 and y ∈ P 1,1 . The following relation holds true: We now turn to a consequence of our additional Hypothesis 6.4 on embedding properties of the Hilbert space H defined above: Proof. Consider a step function f = ij a ij 1 ∆ ij related to a partition (∆ ij ) ij of [0, 1] 2 . We have The general case now easily follows by density of the step functions in H.
Let us now recall that we work under the usual assumptions for Skorohod type change of variables formulae given at Definition 1.6 and referred to as (GC) condition in the sequel. |f (x s;t )| r < ∞, for all r ≥ 1.
We now state an approximation result in H which proves to be useful in order to get our Itô type formula. Proposition 6.8. Let x be a centered Gaussian process on [0, 1] satisfying Hypothesis 6.4 and f ∈ C 1 (R) such that the growth condition (GC) is fulfilled for ϕ and ϕ (1) . Consider a rectangle ∆ = [s 1 , s 2 ] × [t 1 , t 2 ] and π 1 = (s i ) i , π 2 = (t j ) j two respective dissections of the intervals [s 1 , s 2 ] and [t 1 , t 2 ]. Then where we have used Notation 6.1 for the rectangles ∆ i,j .
Proof. Observe first that from which the following estimation is easily obtained: Hence if we take expectations in this last estimation and resort to Hölder's inequality, we obtain that Now the r.h.s of this inequality goes to zero when the mesh of the partitions π 1 , π 2 goes to zero by continuity properties of x (see (70)). Our claim thus easily stems from the embedding (72).
6.2. Itô-Skorohod type formula. We now turn to one of the main aim of this article, namely the proof of a Skorohod type change of variable formula for a general Gaussian process x defined on [0, 1] 2 , under our assumptions 6.4. Our starting point is the relation between z 1 and its Skorohod equivalent.
Proof. Consider a sequence of partitions π n = (π 1 n , π 2 n ) whose mesh go to 0 as n → ∞. The generic elements of π n will be denoted by (s i , t j ). Owing to formula (61), we have that We now treat those two terms separately.
The L 2 (Ω) convergence of the Riemann sums defining A n 1 is more cumbersome, and we have to go back to the definition of the Young integral given by Theorem 4.3. Indeed one can write s 2 s 1 t 2 t 1 y 1 s;t d 12 x s;t = πn ∆ ij y 1 s;t d 12 x s;t , and thanks to (36) we get that s 2 Furthermore, some partial summations can be performed on the terms Id − Λ i δ i for i = 1, 2 and setting D n 1 ≡ s 2 s 1 t 2 t 1 y 1 s;t d 12 x s;t − A n 1 we deduce that Recalling the Hölder regularity (70) of our process x, we thus obtain Now taking expectations in the last relation and using Hölder's inequality we end up with for n large enough.
Step 2: Estimation of A n 2 . Recall that A n 2 is defined by In order to treat this term, first remark that for k = 1, 2 we have Injecting this relation in the definition of the term A n 2 and recalling that we have set R k a ≡ R k aa , we obtain We will now show that where the limits are understood in the almost sure and L 2 sense. Indeed, it is easily understood that the terms A n 22 , A n 23 , A n 24 are remainder terms: according to Hypothesis 6.4 we have that |ρ i ab | |a − b| γ i , and we get the following inequality for A n 22 : This relation, plus the condition (GC) on f , obviously entails that lim n→∞ A n 22 = 0 in the almost sure and L 2 (Ω) sense. The case of A n 23 , A n 24 follow exactly along the same lines. We now focus on the term A n 21 : observe that Invoking the same estimates as before for the Hölder norm of x and condition (GC) on f , the proof of our assertion (77) is now completed.
Step 3: Conclusion. Let us summarize the results obtained in the last two steps: plugging relation (77) into the definition of A n 2 and recalling the limiting behavior of A n 1 established at Step 1, we have obtained that where the convergence is understood in both a.s and L 2 (Ω) sense. Furthermore, Proposition 6.8 asserts that πn y 1 s i ;t j 1 ∆ ij converges in L 2 (Ω, H) to y 1 · 1 ∆ . This finishes our proof of relation (74) thanks to a direct application of Lemma 6.2.
As far as expression (12) with Wick-Riemann sums is concerned, recall that we have proved that δ ⋄ (y 1 · 1 ∆ ) = lim Now invoke Proposition 6.3 for k = 1 in order to state that , which ends the proof. Proposition 6.9 gives a meaning to the increment z 1,⋄ and compares them to the corresponding Stratonovich increment z 1 . In order to compare change of variables formulae, we still have to define Skorohod integrals of the form z 2,⋄ , which is what we proceed to do now.
To this aim, let us start by some formal considerations: it is easily conceived that where, similarly to [17], we set where we integrate firstly in (u ′ , v) and then in (u, v ′ ), and where the notation δ ⋄,2 specifies that we perform double integrals in the Skorohod sense. Our objective in what follows is to give a rigorous meaning to equation (78).
Lemma 6.10. Take up the notation of Proposition 6.9, and consider f ∈ C 3 (R) satisfying condition (GC). For a sequence of partitions (π n ) n≥1 whose mesh goes to 0 define Then a πn converges to N(y) in L 2 (Ω, H ⊗2 ) as n goes to infinity.
Proof. First notice that the tensor norm of an element K ∈ H ⊗2 can be bounded as: and thus, Our claims are now easily derived: on the one hand the right hand side of (81) converges to zero when n → ∞ if u ′ = s i and v ′ = t j for all i, j. Then using inequality (80) and dominated convergence we obtain that a πn converges a.s to N(y) in H ⊗2 . On the other hand, in order to obtain the convergence in L 2 (Ω, H ⊗2 ) it suffices to use the fact that f satisfies condition (GC) and apply once again dominated convergence.
Now we are able to define our mixed integral in the Skorohod sense and connect it to the equivalent integral in the Young theory: Proposition 6.11. Assume x is a centered Gaussian process on [0, 1] 2 with a covariance function satisfying (9). Consider a function ϕ ∈ C 4 (R) satisfying condition (GC) and a rectangle ∆ = [s 1 , s 2 ] × [t 1 , t 2 ]. Then we have that N(y) ∈ Dom(δ ⋄,2 ), and if we define z 2,⋄ = δ ⋄,2 (N(y)) the following relation holds: ;v are defined according to Proposition 4.6. Moreover, relation (13) holds true in the L 2 (Ω) and almost sure sense.
Proof. Like for Proposition 6.9, our strategy is as follows: consider a sequence π n = (π 1 n , π 2 n ) whose mesh go to 0 as n → ∞ and set a n ≡ a πn defined by (79). We have seen at Lemma 6.10 that lim n→∞ a n = N(y) in L 2 (Ω, H ⊗2 ). We shall now study the convergence of δ ⋄ (a n ) by means of Wick-Stratonovich corrections. Then we will conclude by invoking Proposition 6.2.
Step 1: Wick-Stratonovich corrections. According to relation (67) and Proposition 6.3 for k = 2 we obtain δ ⋄,2 (a n ) = πn δ ⋄,2 (y 2 We now use Theorem 4.10 in [11] in order to get that δ ⋄,2 (a n ) can be decomposed as: Like in the proof of Proposition 6.9, we treat those 4 terms separately.
Step 2: Estimation of B n 1 , . . . , B n 4 . The term B n 1 can be decomposed as Moreover, the second term in the r.h.s is the same as A n 2 in the proof of Proposition 6.9, while the convergence for πn y 2 s i ;t j δ 1 x s i s i+1 ;t j δ 2 x s i ;t j t j+1 follows exactly along the same lines as A n 1 in the same proof. We thus leave to the patient reader the task of showing that and we concentrate now on the other terms in (84). The term B n 2 = πn y 3 where we recall that we have set ρ k It is now easily seen that the almost sure and L 2 convergence of B n 2 are obtained with the same kind of considerations as for A n 2 in the proof of Proposition 6.9. We get that and B n 3 , B n 4 are also handled in the same way.
Step 3: Conclusion. Thanks to Step 1 and Step 2, we have obtained that δ ⋄ (a n ) converges to the right hand side of relation (82) as n → ∞, in both almost sure and L 2 senses. As mentioned before, this limiting behavior plus the convergence of a n to N(y) established at Lemma 6.10 yield relation (82) by a direct application of Proposition 6.2. Furthermore, relation (13) is also a direct consequence of relation (83).
Notice that our formula (82) involves some mixed integrals of the form which are defined as Young type integrals. The following proposition, whose proof is similar to Propositions 6.9 and 6.11 and is left to the reader for sake of conciseness, gives a meaning to the analogue integrals in the Skorohod setting.
Proposition 6.12. Let f ∈ C 4 (R) be a function satisfying condition (GC). Then for every where δ ⋄,u is the divergence operator associated to the process (x u;v ) v∈[0,t] . We can thus define where the convergence holds in both L 2 (Ω) and almost sure senses. In addition, we have the following identity: Finally, the integral Remark 6.13. We have defined all the integrals we needed in order to prove our Skorohod change of variable formula (14). Indeed, the proof of formula (14) is now easily deduced by injecting the identities of Propositions 6.9, 6.11 and 6.12 in the Stratonovich type formula (3).

Skorohod's calculus in the rough case
Our goal in this section is to extend the formulae given in Propositions 6.9 and 6.11 to rougher situations, namely for Gaussian processes in the plane with Hölder regularities smaller than 1/2. This is however a harder task than in the Young case, and this is why we introduce 2 simplifications in our considerations: (1) Instead of dealing with a general centered Gaussian process whose covariance admits the factorization property of Hypothesis 6.4, we handle here the case of a fractional Brownian sheet (x s;t ) (s,t)∈[0,1] 2 with Hurst parameters γ 1 , γ 2 ∈ (1/3, 1/2]. (2) The definition of our Skorohod integrals with respect to x is obtained in the following way: we first regularize x as a smooth process x n . For this process we can still use the formulae of Propositions 6.9 and 6.11 like in the Young case. We shall then perform a limiting procedure on these formulae (this is where the specification of a concrete approximation is important), which will give our Stratonovich-Skorohod corrections. Notice however that the interpretation in terms of Riemann-Wick sums will be lost with this strategy.
As in the previous section, we start our considerations by specifying the Malliavin framework in which we are working.

Further Malliavin calculus tools.
Recall that the covariance function of our fractional Brownian sheet x is given by (8). We can thus consider a Hilbert space H x related to x exactly as in Section 6.1, where we now stress the dependence in x of H x in order to differentiate it from the Hilbert space related to white noise. In particular we denote by I x 1 the isometry between H x and the first chaos generated by x.
However, the Malliavin structure related to the harmonizable representation of x will also play a prominent role in the sequel. Namely, it is well known (see e.g. [16]) that for s, t ∈ [0, 1], x can be represented as where c γ 1 ,γ 2 is a normalization constant whose exact value is irrelevant for our computations, W is the Fourier transform of the white noise on R 2 , and Q s;t is a kernel defined by This induces us to consider the canonical Hilbert space related toŴ , that is HŴ = L 2 (R 2 ). The relations between Malliavin calculus with respect toŴ and x are then summarized in the next lemma: Lemma 7.1. Denote by D x,k,p (resp. DŴ ,k,p ) the Sobolev spaces related to x (resp.Ŵ ), and recall the notation L 1,2 = DŴ ,1,2 (L 2 (R)) borrowed from [14]. For φ : [0, 1] 2 → R, set where we recall that Q is defined by (87). Then the following holds true: (i) We can represent the space H x as the closure of the set of step functions under the norm 2 ). In addition, for any smooth function F and any H x -valued square integrable random variable u the following identity holds: (iii) As far as divergence operators are concerned, the relation is Dom(δ x,⋄ ) = K −1 Dom(δŴ ,⋄ ), and δ x,⋄ (u) = δŴ ,⋄ (Ku).
] be a step function. We have that: which easily yields our first claim (i). Let now F be a smooth functional of x of the form F = f (x s 1 ;t 1 , . . . , x sn;tn ). Then and since K1 [0,s l ]×[0,t l ] = Q s l ,t l we end up with which gives our assertion (ii) by density of smooth functionals. Relation (iii) is easily derived from (ii) by duality.
Notice that the preceding result can be extended to second order derivatives thanks to a simple tensorization trick. We label here the result for further use: Then for any smooth functional F and any (H x ) ⊗2 -valued square integrable random variable u we have:

7.2.
Embedding results. Similarly to [2], we now give an embedding result for the space H x which proves to be useful for further computations.
Proof. In this proof we only consider the case γ 1 , γ 2 < 1/2. Indeed, if γ 1 = 1/2 or γ 2 = 1/2 then our process x is simply a Brownian motion in the first or in the second direction, and this situation is handled by L 2 norms.
Let us now introduce two more semi-norms: With these notations in hand, the following embedding result is easily deduced from Lemma 7.3.
Corollary 7.4. Let γ 1 , γ 2 ∈ (0, 1/2) and u ∈ P α 1 ,α 2 1,1 such that 0 < 1 2 − α i < γ i . Then we have the following embedding: 7.3. Strategy and preliminary results. The strategy we shall develop in order to extend Proposition 6.9 (and also Proposition 6.11) to the rough case is based on a regularization of x. Specifically, for a strictly positive integer n, set x n s;t = c γ 1 ,γ 2 |ξ|,|η|≤n where we recall that x and Q are respectively defined by (86) and (87). For fixed n, it is readily checked that x n is a regular Gaussian process. Its covariance function is given by R n s 1 s 2 ;t 1 t 2 = R 1,n s 1 s 2 R 2,n t 1 t 2 , where and hence R n is a regular function which satisfies Hypothesis 6.4. One can thus apply Proposition 6.9 and obtains the following Skorohod-Stratonovich comparison: for ϕ ∈ C 6 (R) satisfying condition (GC), and where we set again ∆ = [s 1 , Our goal is now to take limits in equation (96). A first observation in this direction is that equation (96) involves Skorohod integrals with respect to x n . The fact that a different integral has to be defined for each n is somehow clumsy, and this is why we have decided to express all integrals with respect toŴ in the remainder of our computations. Namely, the same computations as for equations (89) and (90) entail that δ x n ,⋄ (y n ) = δŴ ,⋄ (K n y n ), where K n is the operator defined by With this representation in hand, our limiting procedure can be decomposed as follows: • Take L 2 limits in the right hand side of equation (96) by means of rough paths techniques.
• Show that K n y n converges in L 2 (Ω, L 2 (R)) to Ky.
Thanks to the closability of δŴ ,⋄ , this will show the convergence of δ x n ,⋄ (y n 1 [0,1] 2 ) to δ x,⋄ (y1 [0,1] 2 ) and our Skorohod-Stratonovich correction formula will be obtained in this way. We now state and prove 3 useful lemmas for our future computations. The first one deals with convergence of covariance functions: where R i,n is defined by (95).
Proof. The definitions (86) of x, (94) of x n plus formula (87) for Q allow to write, for all n ≥ 1: where we have set I n = |x|,|y|≥n dx dy |x| 1+ǫ |y| 1+ǫ (this quantity is obviously finite). Hence by Gaussian hypercontractivity and Lemma 6.6 we obtain Since lim n→∞ I n = 0, we thus get lim n→∞ E[[ x n − x p γ 1 −ǫ,γ 2 −ǫ ] = 0, which is our first claim. We now focus on the exponential integrability of sup x n . Notice that for a fixed n one can easily get those exponential estimates thanks to Fernique's lemma. However, we claim some uniformity in n here, and we thus come back to uniform estimates of moments in order to prove (98). Let then r = max(⌊ 1 γ 1 −ǫ ⌋, ⌊ 1 γ 2 −ǫ ⌋) + 1 and remark that x n ∞ ≤ x n ǫ,ǫ . We thus use a decomposition of the form E[e λ sup (s,t)∈[0,1] 2 |x n s;t | 2 ] = I 1,n (λ) + I 2,n (λ), where We now bound those 2 terms separately: one the one hand, it is readily checked that for ǫ < min(γ 1 , γ 2 ). On the other hand, the bound on I 2,n (λ) is obtained invoking Lemma 6.6 again. Indeed, starting from expression (71) and introducing a standard Gaussian random variable N , it is easily seen that Furthermore, it can be shown that C 2l,ǫ,ǫ is of the form M l for a given M > 1. Thus from which the relation I 2,n (λ) < ∞ is easily obtained. This finishes the proof of (98).
With classical considerations concerning compositions of Hölder functions with non linearities, we finally get the following result which is labelled for further use. Its proof is omitted for sake of conciseness.

7.4.
Itô-Skorohod type formula. We now turn to the limiting procedure in equation (96), beginning with the term involving covariances only: Proposition 7.8. Let f ∈ C 6 (R) be a function satisfying condition (GC) with a small parameter λ > 0, and x n be the regularized version of x defined by (94). Then the following convergence: holds in L 2 (Ω).
Proof. The integrals involved in (101) are all of Young type. Owing to Proposition 4.6, we thus have: By continuity of the sewing map, the desired convergence will thus stem from the relations lim n→0 A 1,n = 0 and lim n→0 A 2,n = 0, where for ǫ > 0 we set: Now the relation lim n→0 A 1,n = 0 is obviously a direct consequence of Lemma 7.5. As far as A 2,n is concerned, we start from relation (100) and apply Hölder's inequality. This yields 16 ]. Then according to Proposition 7.6 we see that the r.h.s of this last equation vanishes when n goes to infinity, which proves our claim.
Proof. Let us start from the corrections for the regularized process x n , for which we can appeal to Proposition 6.9. We obtain relation (96), written here again for convenience: z 1,n,⋄ s 1 s 2 ;t 1 t 2 ≡ δ x n ,⋄ (y n,1 · 1 ∆ ) = ∆ y n,1 s;t d 12 x n s;t − 1 4 ∆ y n,2 s;t d 1 R 1,n s d 2 R 2,n t .
Now putting together Proposition 7.8 and the continuity of the rough path integral, we get convergence of the r.hs of (103) in L 2 (Ω). Thus one can write, in the a.s and L 2 (Ω) sense: where the integral with respect to x is interpreted in the sense of Theorem 1.4. Let us further analyze the convergence of δ x n ,⋄ (y n ): recall that this quantity can be written as δŴ ,⋄ (K n y n ), where K n is defined by (97) or specifically as (K n y n )(ξ, η) = ıξ ıη |ξ| γ 1 +1/2 |η| γ 2 +1/2 ∆ y n uv e ıξu+ıηv dudv 1 (|ξ|,|η|≤n) Hence, owing to closability of the operator δŴ ,⋄ , the proof of (102) is reduced to show that K n y n,1 converges in L 2 (Ω; L 2 (∆)) to Ky 1 . Now expression (104) easily entails that On the other hand, Corollary 7.4 asserts that y − y n H N γ 1 −ǫ,γ 2 −ǫ (y n − y), and the r.h.s of this relation vanishes as n → ∞ thanks to Proposition 7.6. This concludes our proof.
In order to complete our comparison between Itô and Stratonovich formulae, we still have to compare the Skorohod type increment z 2,⋄ and the rough integral z 2 . As a previous step, let us give an intermediate result concerning some mixed integrals in R, x: Proposition 7.10. Let ϕ ∈ C 6 (R) and recall that for the fractional Brownian sheet x we have R i u = u 2γ i for i = 1, 2. Then the integral y 1 R 1 d 1 x d 2 R 2 = [(Id − Λ 1 δ 1 )(Id − Λ 2 δ 2 )] y 1 R 1 δ 1 xδ 2 R 2 + 1/2y 2 R 1 (δ 1 x) 2 δ 2 R 2 (105) is well defined a.s, in the sense of Proposition 3.3. Moreover the following convergence takes place in L 2 (Ω): Finally, the same kind of result is still verified when one interchanges directions 1 and 2 in relation (105).
To begin with, note that the convergence (108) easily stems from Lemma 7.5 for the terms R, Proposition 7.6 for the terms x and Lemma 7.7 for the terms y. In order to prove (109), we invoke again the integral representation (107) for both h n and h. Then some elementary considerations (omitted here for sake of conciseness) allow to reduce the problem to the following relation: This last relation is a direct consequence of Proposition 7.6 and composition with non linearities, whenever f satisfies the growth condition (GC) with a small parameter λ > 0. The proof is now finished.
The last step in order to go from Theorem 7.11 to Theorem 1.8 is to convert Stratonovich into Skorohod type integrals in the right hand side of relation (110). To this aim, let us first recall the following one-parameter result: Proposition 7.12. Let B a fractional brownian motion with hurst parameters 1/2 ≥ γ > 1/3 then we have that u → ϕ(B u )u 2γ ∈ Dom(δ ⋄,B ) and we have that Proof. Use exactly the same arguments of the Proposition (7.9) for the one parameters setting. Now the Corollary below is the key to the conversion of Theorem 7.11 into Theorem 1.8: Corollary 7.13. For γ i > 1/3 and ϕ ∈ C 6 (R) then for every v ∈ [0, 1] u → y 3 u;v u 2γ 1 ∈ Dom(δ ⋄,x.;v ) and the following formula hold true Proof. we recall that x u,v = law v γ 2 B u with B is a fBm with hurst parameter γ 1 and then it suffice to use the proposition (7.4)