A general Multidimensional Monte Carlo Approach for Dynamic Hedging under stochastic volatility

In this work, we introduce a Monte Carlo method for the dynamic hedging of general European-type contingent claims in a multidimensional Brownian arbitrage-free market. Based on bounded variation martingale approximations for Galtchouk-Kunita-Watanabe decompositions, we propose a feasible and constructive methodology which allows us to compute pure hedging strategies w.r.t arbitrary square-integrable claims in incomplete markets. In particular, the methodology can be applied to quadratic hedging-type strategies for fully path-dependent options with stochastic volatility and discontinuous payoffs. We illustrate the method with numerical examples based on generalized Follmer-Schweizer decompositions, locally-risk minimizing and mean-variance hedging strategies for vanilla and path-dependent options written on local volatility and stochastic volatility models.

1. Introduction 1.1. Background and Motivation. Let (S, F, P) be a financial market composed by a continuous F-semimartingale S which represents a discounted risky asset price process, F = {F t ; 0 ≤ t ≤ T } is a filtration which encodes the information flow in the market on a finite horizon [0, T ], P is a physical probability measure and M e is the set of equivalent local martingale measures. Let H be an F Tmeasurable contingent claim describing the net payoff whose the trader is faced at time T . In order to hedge this claim, the trader has to choose a dynamic portfolio strategy.
Under the assumption of an arbitrage-free market, the classical Galtchouk-Kunita-Watanabe (henceforth abbreviated by GKW) decomposition yields where L H,Q is a Q-local martingale which is strongly orthogonal to S and θ H,Q is an adapted process. The GKW decomposition plays a crucial role in determining optimal hedging strategies in a general Brownian-based market model subject to stochastic volatility. For instance, if S is a one-dimensional Itô risky asset price process which is adapted to the information generated by a two-dimensional Brownian motion W = (W 1 , W 2 ), then there exists a two-dimensional adapted process φ H,Q := (φ H,1 , φ H,2 ) such that which also realizes In the complete market case, there exists a unique Q ∈ M e and in this case, L H,Q = 0, E Q [H] is the unique fair price and the hedging replicating strategy is fully described by the process θ H,Q . In a general stochastic volatility framework, there are infinitely many GKW orthogonal decompositions parameterized by the set M e and hence one can ask if it is possible to determine the notion of non-selffinancing optimal hedging strategies solely based on the quantities (1.2). This type of question was firstly answered by Föllmer and Sonderman [9] and later on extended by Schweizer [23] and Föllmer and Schweizer [8] through the existence of the so-called Föllmer-Schweizer decomposition which turns out to be equivalent to the existence of locally-risk minimizing hedging strategies. The GKW decomposition under the so-called minimal martingale measure constitutes the starting point to get locally risk minimizing strategies provided one is able to check some square-integrability properties of the components in (1.1) under the physical measure. See e.g [12] and [26] for details and other references therein. Orthogonal decompositions without square-integrability properties can also be defined in terms of the the so-called generalized Föllmer-Schweizer decomposition (see e.g Schweizer [24]).
In contrast to the local-risk minimization approach, one can insist in working with self-financing hedging strategies which give rise to the so-called mean-variance hedging methodology. In this approach, the spirit is to minimize the expectation of the squared hedging error over all initial endowments x and all suitable admissible strategies ϕ ∈ Θ: The nature of the optimization problem (1.3) suggests to work with the subset M e 2 := {Q ∈ M e ; dQ dP ∈ L 2 (P)}. Rheinlander and Schweizer [22], Gourieroux, Laurent and Pham [10] and Schweizer [25] show that if M e 2 = ∅ and H ∈ L 2 (P) then the optimal quadratic hedging strategy exists and it is given by Here θ H,P is computed in terms ofP, the so-called variance optimal martingale measure,ζ realizes (1.5)Z t := EP dP dP F t =Z 0 + t 0ζ dS ; 0 ≤ t ≤ T, and V H,P := EP[H|F · ] is the value option price process underP. See also Cerný and Kallsen [4] for the general semimartingale case and the works [16], [18] and [19] for other utility-based hedging strategies based on GKW decompositions. Concrete representations for the pure hedging strategies {θ H,Q ; Q =P,P} can in principle be obtained by computing cross-quadratic variations d[V H,Q , S] t /d[S, S] t for Q ∈ {P,P}. For instance, in the classical vanilla case, pure hedging strategies can be computed by means of the Feynman-Kac theorem (see e.g Heath, Platen and Schweizer [12]). In the path-dependent case, the obtention of concrete computationally efficient representations for θ H,Q is a rather difficult problem. Feynman-Kac-type arguments for fully path-dependent options mixed with stochastic volatility typically face not-well posed problems on the whole trading period as well as highly degenerate PDEs arise in this context. Generically speaking, one has to work with non-Markovian versions of the Feynman-Kac theorem in order to get robust dynamic hedging strategies for fully path dependent options written on stochastic volatility risky asset price processes.
In the mean variance case, the only quantity in (1.4) not related to GKW decomposition isZ which can in principle be expressed in terms of the so-called fundamental representation equations given by Hobson [14] and Biagini, Guasoni and Pratelli [2] in the stochastic volatility case. For instance, Hobson derives closed form expressions forζ and also for any type of q-optimal measure in the Heston model [13]. Recently, semi-explicit formulas for vanilla options based on general characterizations of the variance-optimal hedge in Cerný and Kallsen [4] have been also proposed in the literature which allow for a feasible numerical implementation in affine models. See Kallsen and Vierthauer [17] and Cerný and Kallsen [5] for some results in this direction. A different approach based on backward stochastic differential equations can also be used in order to get useful characterizations for the optimal mean variance hedging strategies. See e.g Jeanblanc, Mania, Santacrose and Schweizer [15] and other references therein.
1.2. Contribution of the current paper. In spite of deep characterizations of optimal quadratic hedging strategies and concrete numerical schemes available for vanilla-type options, to our best knowledge no feasible approach has been proposed to tackle the problem of obtaining dynamic optimal quadratic hedging strategies for fully path dependent options written on a generic multidimensional Itô risky asset price process. In this work, we attempt to solve this problem with a probabilistic approach. The main difficulty in dealing with fully path dependent and/or discontinuous payoffs is the non-Markovian nature of the option value and a priori lack of path smoothness of the pure hedging strategies. Usual numerical schemes based on PDE and martingale techniques do not trivially apply to this context.
The main contribution of this paper is the obtention of flexible and computationally efficient multidimensional non-Markovian representations for generic option price processes which allow for a concrete computation of the associated GKW decomposition θ H,Q , L H,Q for Q-square integrable payoffs H with Q ∈ M e . We provide a Monte Carlo methodology capable to compute optimal quadratic hedging strategies w.r.t general square-integrable claims in a multidimensional Brownian-based market model.
This article provides a feasible and constructive method to compute generalized Föllmer-Schweizer decompositions under full generality. As far as the mean variance hedging is concerned, we are able to compute pure optimal hedging strategies θ H,P for arbitrary square-integrable payoffs. Hence, our methodology also applies to this case provided one is able to compute the fundamental representation equations in Hobson [14] and Biagini, Guasoni and Pratelli [2] which is the case for the classical Heston model. In mathematical terms, we are able to compute Q-GKW decompositions under full generality so that the results of this article can also be used to other non-quadratic hedging methodologies where orthogonal martingale representations play an important role in determining optimal hedging strategies.
The starting point of this article is based on weak approximations developed by Leão and Ohashi [20] for one-dimensional Brownian functionals. They introduced a one-dimensional space-filtration discretization scheme constructed from suitable waiting times which measure the instants when the Brownian motion hits some a priori levels. In this work, we extend [20] to the multidimensional case as follows: More general and stronger convergence results are obtained in order to recover incomplete markets with stochastic volatility. Hitting times induced by multidimensional noises which drive the stochastic volatility are carefully analyzed in order to obtain Q-GKW decompositions under rather weak integrability conditions for any Q ∈ M e . Moreover, a complete analysis is performed w.r.t weak approximations for gain processes by means of suitable non-antecipative discrete-time hedging strategies for square-integrable payoffs, including path-dependent ones.
It is important to stress that the results of this article can be applied to both complete and incomplete markets written on a generic multidimensional Itô risky asset price process. One important restriction of our methodology is the assumption that the risky asset price process has continuous paths. This is a limitation that we hope to overcome in a future work.
Numerical results based on the standard Black-Scholes, local-volatility and Heston models are performed in order to illustrate the theoretical results and the methodology of this article. In particular, we briefly compare our results with other prominent methodologies based on Malliavin weights (complete market case) and PDE techniques (incomplete market case) employed by Bernis, Gobet and Kohatsu-Higa [1] and Heath, Platen and Schweizer [12], respectively. The numerical experiments suggest that pure hedging strategies based on generalized Föllmer-Schweizer decompositions mitigate very well the cost of hedging of a path-dependent option even if there is no guarantee of the existence of locally-risk minimizing strategies. We also compare hedging errors arising from optimal mean variance hedging strategies for one-touch options written on a Heston model with nonzero correlation.
The remainder of this paper is structured as follows. In Section 2, we fix the notation and we describe the basic underlying market model. In Section 3, we provide the basic elements of the Monte Carlo methodology proposed in this article. In Section 4, we formulate dynamic hedging strategies starting from a given GKW decomposition and we translate our results to well-known quadratic hedging strategies. The Monte Carlo algorithm and the numerical study are described in Sections 5 and 6, respectively. The Appendix presents more refined approximations when the martingale representations admit additional hypotheses.

Preliminaries
Throughout this paper, we assume that we are in the usual Brownian market model with finite time horizon 0 ≤ T < ∞ equipped with the stochastic basis (Ω, F, P) generated by a standard p-dimensional is the P-augmentation of the natural filtration generated by B. For a given m-dimensional vector J = (J 1 , . . . , J m ), we denote by diag(J) the m × m diagonal matrix whose -th diagonal term is J . In this paper, for all unexplained terminology concerning general theory of processes, we refer to Dellacherie and Meyer [6].
In view of stochastic volatility models, let us split B into two multidimensional Brownian motions as follows B S := (B (1) , . . . , B (d) ) and B I := (B (d+1) , . . . , B (p) ). In this section, the market consists of d + 1 assets (d ≤ p): one riskless asset given by , and a d-dimensional vector of risky assetsS := (S 1 , . . . ,S d ) which satisfies the following stochastic differential equation Here, the real-valued interest rate process r = {r t ; 0 ≤ t ≤ T }, the vector of mean rates of return b := {b t = (b 1 t , . . . , b d t ); 0 ≤ t ≤ T } and the volatility matrix σ := {σ t = (σ ij t ); 1 ≤ i ≤ d, 1 ≤ j ≤ d, 0 ≤ t ≤ T } are assumed to be predictable and they satisfy the standard assumptions in such way that both S 0 andS are well-defined positive semimartingales. We also assume that the volatility matrix σ is non-singular for almost all (t, ω) ∈ [0, T ] × Ω. The discounted price S := {S i :=S i /S 0 ; i = 1, . . . , d} follows where 1 d is a d-dimensional vector with every component equal to 1. The market price of risk is given by In the sequel, M e denotes the set of P-equivalent probability measures Q such that the respective Radon-Nikodym derivative process is a P−martingale and the discounted price S is a Q-local martingale. Throughout this paper, we assume that M e = ∅. In our setup, it is well known that M e is given by the subset of probability measures with Radon-Nikodym derivatives of the form for some R p−d -valued adapted process ν such that

Example:
The typical example studied in the literature is the following one-dimensional stochastic volatility model where Y (1) and Y (2) are correlated Brownian motions with correlation ρ ∈ [−1, 1], µ, a and b are suitable functions such that (S, σ 2 ) is a well-defined two-dimensional Markov process. All continuous stochastic volatility model commonly used in practice fit into the specification (2.1). In this case, p = 2 > d = 1 and we recall that the market is incomplete where the set M e is infinity. The dynamic hedging procedure turns out to be quite challenging due to extrinsic randomness generated by the non-tradeable volatility, specially w.r.t to exotic options.

GKW Decomposition.
In the sequel, we take Q ∈ M e and we set W S := (W (1) , . . . , W (d) ) and W I := (W (d+1) , . . . , W (p) ) where is a standard p-dimensional Brownian motion under the measure Q and filtration F : In what follows, we fix a discounted contingent claim H. Recall that the filtration F is contained in F, but it is not necessarily equal. In the remainder of this article, we assume the following hypothesis.
(M) The contingent claim H is also F T -measurable. For a given Q-square integrable claim H, the Brownian martingale representation (computed in terms of (F, Q)) yields The discounted stock price process has the following Q-dynamics dS t = diag(S t )σ t dW S t , S 0 = x, 0 ≤ t ≤ T, and therefore the Q-GKW decomposition for the pair of locally square integrable local martingales (V , S) is given byV The p-dimensional process φ H,Q which constitutes (2.3) and (2.5) plays a major role in several types of hedging strategies in incomplete markets and it will be our main object of study.
Remark 2.2. If we set ν j = 0 for j = d+1, . . . , p and the correspondent density process is a martingale then the resulting minimal martingale measureP yields a GKW decomposition where L H,P is still a P-local martingale orthogonal to the martingale component of S under P. In this case, it is also natural to implement a pure hedging strategy based on θ H,P regardless the existence of the Föllmer-Schweizer decomposition. If this is the case, this hedging strategy can be based on the generalized Föllmer-Schweizer decomposition (see e.g Th.9 in [24]).

The Random Skeleton and Weak Approximations for GKW Decompositions
In this section, we provide the fundamentals of the numerical algorithm of this article for the obtention of hedging strategies in complete and incomplete markets.
3.1. The Multidimensional Random Skeleton. At first, we fix once and for all Q ∈ M e and a Q-square-integrable contingent claim H satisfying (M). In the remainder of this section, we are going to fix a Q-Brownian motion W and with a slight abuse of notation all Q-expectations will be denoted by E. The choice of Q ∈ M e is dictated by the pricing and hedging method used by the trader.
In the sequel, [·, ·] denotes the usual quadratic variation between semimartingales and the usual jump of a process is denoted by ∆Y t = Y t − Y t− where Y t− is the left-hand limit of a cadlag process Y . For a pair (a, b) ∈ R 2 , we denote a ∨ b := max{a, b} and a ∧ b := min{a, b}. Moreover, for any two stopping times S and J, we denote the stochastic intervals For a fixed positive integer k and for each j = 1, 2, . . . , p we define T k,j 0 := 0 a.s. and In the sequel, we define A k := (A k,1 , . . . , A k,p ) as the p-dimensional step process given componentwise by for k, n ≥ 1 and j = 1, . . . , p. We split Here, denotes the smallest sigma-algebra generated by the union. One can easily check that F k,j With a slight abuse of notation we write F k,j t to denote its Q-augmentation satisfying the usual conditions.
Let us now introduce the multidimensional filtration generated by A k . Let us consider for n ≥ 1. In this case, T k is the partition generated by all stopping times defined in (3.1). The finite-dimensional distribution of W (j) is absolutely continuous for each j = 1, . . . , p and therefore the elements of T k are almost surely distinct for every k ≥ 1. Moreover, the following result holds true.
Lemma 3.1. For every k ≥ 1, the set T k is an exhaustive sequence of F k -stopping times such that sup n≥1 |T k n − T k n−1 | → 0 in probability as k → ∞.
Proof. The following obvious estimate holds in probability as k → ∞ and T k n → ∞ a.s as n → ∞ for each k ≥ 1. Let us now prove that T k = {T k n ; n ≥ 0} is a sequence of F k -stopping times. In order to show this, we write (T k n ) n≥0 in a different way. This sequence can be defined recursively as follows Next, let us define a family of F k T k 1 -random variables related to the index j which realizes the hitting time T k 1 as follows k,j for any j = 1, . . . , p. Then, we shift A k as follows for t ≥ 0. In this case, we conclude thatÃ k 1 is adapted to the filtration {F k t+T k 1 ; t ≥ 0}, the hitting time In the sequel, we define a family of F k T k 2 -random variables related to the index j which realizes the hitting time T k 2 as follows for every t ≥ 0. In this case, we conclude thatÃ k 2 is adapted to the filtration {F k t+T k 2 ; t ≥ 0}, the hitting time By induction, we conclude that (T k n ) n≥0 is a sequence of F k -stopping times. With Lemma 3.1 at hand, we notice that the filtration F k is a discrete-type filtration in the sense that We now embed the process X into the quasi left-continuous filtration F k by means of the following operator Since X is a F-martingale, then the usual optional stopping theorem yields the representation Therefore, δ k X is indeed a Q-square-integrable F k -martingale and we shall write it as and the integral in (3.4) is computed in the Lebesgue-Stieltjes sense.
Remark 3.1. Similar to the univariate case, one can easily check that F k → F weakly and since X has continuous paths then δ k X → X uniformly in probability as k → ∞. See Remark 2.1 in [20].
Based on the Dirac process D j δ k X, we denote In order to work with non-antecipative hedging strategies, let us now define a suitable F k -predictable version of D k,j X as follows One can check that D k,j X is F k -predictable. See e.g [11], Ch.5 for details.
Example: Let H be a contingent claim satisfying (M). Then for a given j = 1, . . . , p, we have One should notice that (3.5) is reminiscent from the usual delta-hedging strategy but the price is shifted on the level of the sigma-algebras jointly with the increments of the driving Brownian motion instead of the pure spot price. For instance, in the one-dimensional case (p = d = 1), we have and hence a natural procedure to approximate pure hedging strategies is to look at D k,1 X T k,1 1 /S 0 σ 0 at time zero. In the incomplete market case, additional randomness from e.g stochastic volatilities are encoded by E[H|F k 1 is determined not only by the hitting times coming from the risky asset prices but also by possibly Brownian motion hitting times coming from stochastic volatility.
In the next sections, we will construct feasible approximations for the gain and cost processes based on the ratios (3.5). We will see that hedging ratios of the form (3.5) will be the key ingredient to recover the gain process in full generality.

3.2.
Weak approximation for the hedging process. Based on (2.3), (2.4) and (2.5), let us denote In order to shorten notation, we do not write (φ H,Q,S , φ H,Q,I ) in (3.6). The main goal of this section is the obtention of bounded variation martingale weak approximations for both the gain and cost processes, given respectively, by We assume the trader has some knowledge of the underlying volatility so that the obtention of φ H,S will be sufficient to recover θ H . The typical example we have in mind are generalized Föllmer-Schweizer decompositions, locally-risk minimizing and mean variance strategies as explained in the Introduction. The scheme will be very constructive in such way that all the elements of our approximation will be amenable to a feasible numerical analysis. Under very mild integrability conditions, the weak approximations for the gain process will be translated into the physical measure.
The weak topology. In order to obtain approximation results under full generality, it is important to consider a topology which is flexible to deal with nonsmooth hedging strategies θ H for possibly non-Markovian payoffs H and at the same time justifies Monte Carlo procedures. In the sequel, we make use of the weak topology σ( The subspace of the squareintegrable F-martingales will be denoted by H 2 (F). It will be also useful to work with σ(B 1 , Λ ∞ )topology given in [20]. For more details about these topologies, we refer to the works [6,7,20]. It turns out that σ(B 2 , M 2 ) and σ(B 1 , Λ ∞ ) are very natural notions to deal with generic square-integrable random variables as described in [20].
In the sequel, we recall the following notion of covariation introduced in [20].
. . , p be a sequence of stochastic integrals and Proof. The proof follows easily from the arguments given in the proof of Prop. 3.2 in [20] by using the fact that {W (j) ; 1 ≤ j ≤ p} is an independent family of Brownian motions, so we omit the details.
In the sequel, we present a key asymptotic result for the numerical algorithm of this article.
Proof. We divide the proof into three steps. Throughout this proof C is a generic constant which may defer from line to line.

STEP1.
We claim that for each j = 1, . . . , p. In order to prove (3.10), we begin by noticing that Lemma 3.1 states that the elements of T k are F-stopping times. By assumption, X is Q-square integrable martingale and hence one may use similar arguments given in the proof of Lemma 3.1 in [20] to safely state that the following estimate holds Now, we notice that the sequence F k converges weakly to F, X is continuous and therefore δ k X → X uniformly in probability (see Remark 3.1). Since X ∈ B 2 (F), then a simple application of Burkhölder inequality allows us to state that δ k X converges strongly in B 1 (F) and a routine argument based on the definition of the B 2 -weak topology yields holds weakly in L 1 for each t ∈ [0, T ] and j = 1, . . . , p due to the pairwise independence of {W (j) ; 1 ≤ j ≤ p}. Summing up (3.11) and (3.13), we shall apply Lemma 3.2 to get (3.10).

STEP 2.
In the sequel, let (·) o,k and (·) p,k be the optional and predictable projections w.r.t F k , respectively. Let us consider the F k -martingales given by a.s for each n, k ≥ 1 and j = 1 . . . , p (see e.g chap.5, section 5 in [11]). Moreover, by the very definition (3.14) Therefore, Jensen inequality yields where in (3.15) we have used (3.14) and the fact that (D k,j X) 2 a.s for each n, k ≥ 1 and j = 1 . . . , p. We shall write J k in a slightly different manner as follows The above identities, the estimates (3.11), (3.15) and Remark 3.1 yield STEP 3. We claim that for a given g ∈ L ∞ , t ∈ [0, T ] and j = 1 . . . , p we have By using the fact that D k,j X is F k -optional and D k,j X is F k -predictable, we shall use duality of the F k -optional projection to write In order to prove (3.18), let us check that The same trick we did in (3.16) together with (3.14) yield as k → ∞ because X has continuous paths (see Remark 3.1). This proves (3.19). Now, in order to shorten notation let us denote I k,j by the expectation in (3.20). Cauchy-Schwartz and Burkholder-Davis-Gundy inequalities jointly with (3.17) and (3.11) yield We shall proceed similar to Lemma 4.1 in [20] to safely state that Stronger convergence results can be obtained under rather weak integrability and path smoothness assumptions for representations (φ 1 , . . . , φ p ). We refer the reader to the Appendix for further details.

Weak dynamic hedging
In this section, we apply Theorem 3.1 for the formulation of a dynamic hedging strategy starting with a given GKW decomposition where H is a Q-square integrable European-type option satisfying (M) for a given Q ∈ M e . The typical examples we have in mind are quadratic hedging strategies w.r.t a fully path-dependent option. We recall that when Q is the minimal martingale measure then (4.1) is the generalized Föllmer-Schweizer decomposition so that under some P-square integrability conditions on the components of (4.1), θ H is the locally risk minimizing hedging strategy (see e.g [12], [24]). In fact, GKW and Föllmer-Schweizer decompositions are essentially equivalent for the market model assumed in Section 2. We recall that decomposition (4.1) is not sufficient to fully describe mean variance hedging strategies but the additional component rests on the fundamental representation equations as described in Introduction. See also expression (6.4) in Section 6. For simplicity of exposition, we consider a financial market (Ω, F, P) driven by a two-dimensional Brownian motion B and a one-dimensional risky asset price process S as described in Section 2. We stress that all results in this section hold for a general multidimensional setting with the obvious modifications.
In the sequel, we denote for k, n ≥ 1.  Proof. We have E| dP dQ | 2 = E P | dP dQ | 2 dQ dP = E P dP dQ < ∞. To shorten notation, let Y k t := t 0 D k,1 s XdA k,1 s and Y t := t 0 θ H dS for 0 ≤ t ≤ T. Let G be an arbitrary F-stopping time bounded by T and let g ∈ L ∞ (P) be an essentially P-bounded random variable and F G -measurable. Let J ∈ M 2 be a continuous linear functional given by the purely discontinuous F-optional bounded variation process where the duality action · , · is given by J, N = E T 0 N s dJ s ; N ∈ B 2 (F). See [20] for more details. Then Theorem 3.1 and the fact dP dQ ∈ L 2 (Q) yield Then from the definition of the σ(B 1 , Λ ∞ )-topology based on the physical measure P, we shall conclude the proof.

Remark 4.1. Corollary 4.1 provides a non-antecipative Riemman-sum approximation for the gain process
· 0 θ H t dS t in a multi-dimensional filtration setting where none path regularity of the pure hedging strategy θ H is imposed. The price we pay is a weak-type convergence instead of uniform convergence in probability. However, from the financial point of view this type of convergence is sufficient for the implementation of Monte Carlo methods in hedging. More importantly, we will see that θ k,H can be fairly simulated and hence the resulting Monte Carlo hedging strategy can be calibrated from market data.

Remark 4.2.
If one is interested only at convergence at the terminal time 0 < T < ∞, then assumption (4.2) can be weakened to E P | T 0 θ H t dS t | < ∞. Assumption E P dP dQ < ∞ is essential to change the Q-convergence into the physical measure P. One should notice that the associated density process is no longer a P-local-martingale and in general such integrability assumption must be checked case by case. Such assumption holds locally for every underlying Itô risky asset price process. Our numerical results suggest that this property behaves well for a variety of spot price models.
Of course, in practice both the spot prices and trading dates are not observable at the stopping times so we need to translate our results to a given deterministic set of rebalancing hedging dates.

Hedging Strategies.
In this section, we provide a dynamic hedging strategy based on a refined set of hedging dates Π := 0 = s 0 < . . . < s p−1 < s p = T . For this, we need to introduce some objects. For a given s i ∈ Π, we set W si,· is an (F j si,t ) 0≤t≤T −si -Brownian motion for each j = 1, 2 and independent from F j si , where F j si,t := F j si+t for 0 ≤ t ≤ T − s i . Similar to Section 3.1, we set T k,1 si,0 := 0 and T k,j si,n := inf{t > T k,j si,n−1 ; |W For a given k ≥ 1 and j = 1, 2, we define H k,j si,n as the sigma-algebra generated by {T k,j si, ; 1 ≤ ≤ n} and W ; 1 ≤ ≤ n. We then define the following discrete jumping filtration F k,j si,t := H k,j si,n a.s on {T k,j si,n ≤ t < T k,j si,n+1 }. In order to deal with fully path dependent options, it is convenient to introduce the following augmented filtration G k,j si,t := F j si ∨ F k,j si,t ; 0 ≤ t ≤ T − s i , for j = 1, 2. The bidimensional information flows are defined by F si,t := F 1 si,t ⊗ F 2 si,t and G k si,t := G k,1 si,t ⊗ G k,2 si,t for 0 ≤ t ≤ T − s i . We set G k si := {G k si,t ; 0 ≤ t ≤ T − s i }. We shall assume that they satisfy the usual conditions. The piecewise constant martingale projection A k,j si based on W (j) si is given by si,T −si |G k,j si,t ]; 0 ≤ t ≤ T − s i . We set {T k si,n ; n ≥ 0} as the order statistic generated by the stopping times {T k,j si,n ; j = 1, 2, n ≥ 0} similar to (3.3).
If H ∈ L 2 (Q) and X t = E[H|F t ]; 0 ≤ t ≤ T , then we define δ k si X t := E[H|G k si,t ]; 0 ≤ t ≤ T − s i , so that the related derivative operators are given by An G k si -predictable version of D k,j si X is given by In the sequel, we denote where σ si,· is the volatility process driven by the shifted filtration {F si,t ; 0 ≤ t ≤ T − s i } and S si,· is the risky asset price process driven by the shifted Brownian W   for k ≥ 1and Π. At first, we recall that {T k,1 si,n − T k,1 si,n−1 ; n ≥ 1, s i ∈ Π} is an i.i.d sequence with absolutely continuous distribution. In this case, the probability of the set {T k,1 si,n ≤ s i+1 − s i } is always strictly positive for every Π and k, n ≥ 1. Hence, R(θ k,H , Π, k) is a non-degenerate subset of random variables. By making a change of variable on the Itô integral, we shall write Let us fix Q ∈ M e . By the very definition, si , for each s i ∈ Π with the discretization of the Brownian motion given by A k,1 si . Moreover, using the fact that E| dP dQ | 2 < ∞ and repeating the argument given by  In practice, one may approximate the gain process by a non-antecipative strategy as follows: Let Π be a given set of trading dates on the interval [0, T ] so that |Π| = max 0≤i≤p |s i − s i−1 | is small. We take a large k and we perform a non-antecipative buy and hold-type strategy among the trading dates [s i , s i+1 ); s i ∈ Π in the full approximation (4.6) which results Convergence (4.5) implies that the approximation (4.9) results in unavoidable hedging errors w.r.t the gain process due to the discretization of the dynamic hedging, but we do not expect large hedging errors provided k is large and |Π| small. Hedging errors arising from discrete hedging in complete markets are widely studied in literature. We do not know optimal rebalancing dates in this incomplete market setting, but simulation results presented in Section 6 suggest that homogeneous hedging dates work very well for a variety of models with and without stochastic volatility. A more detailed study is needed in order to get more precise relations between Π and the stopping times, a topic which will be further explored in a forthcoming paper.
Let us now briefly explain how the results of this section can be applied to well-known quadratic hedging methodologies.
Generalized Föllmer-Schweizer: If one takes the minimal martingale measureP, then L H in (4.1) is a P-local martingale and orthogonal to the martingale component of S. Due this orthogonality and the zero mean behavior of the cost L H , it is still reasonable to work with generalized Föllmer-Schweizer decompositions under P without knowing a priori the existence of locally-risk minimizing hedging strategies.
Local Risk Minimization: One should notice that if θ H dS ∈ B 2 (F), L H ∈ B 2 (F) under P and dP dP ∈ L 2 (P), then θ H is the locally risk minimizing trading strategy and (4.1) is the Föllmer-Schweizer decomposition under P.
Mean Variance hedging: If one takesP, then the mean variance hedging strategy is not completely determined by the GKW decomposition underP. Nevertheless, Corollary 4.2 still can be used to approximate the optimal hedging strategy by computing the density processZ based on the so-called fundamental equations derived by Hobson [14]. See

The algorithm
In this section we present the basic algorithm to evaluate the hedging strategy for a given Europeantype contingent claim H ∈ L 2 (Q) satisfying assumption (M) for a fixed Q ∈ M e at a terminal time 0 < T < ∞. The structure of the algorithm is based on the space-filtration discretization scheme induced by the stopping times {T k,j m ; k ≥ 1, m ≥ 1, j = 1, . . . , p}. From the Markov property, the key point is the simulation of the first passage time T k,j 1 for each j = 1 . . . , p for which we refer the work of Burq and Jones [3] for details.
(1) One chooses k ≥ 1 which represents the level of discretization of the Brownian motion.
(3) One simulates the family {η k,j ; ≥ 1} independently from {T k,j − T k,j −1 ; ≥ 1}. This i.i.d family {η k,j ; ≥ 1} must be simulated according to the Bernoulli random variable η k,j 1 with parameter 1/2 for i = −1, 1. This simulates the jump process A k,j for j = 1, . . . , p. The next step is the simulation of D k,j X T k,j 1 where the conditional expectations in (3.5) play a key role. For this, we need to simulate H based on {S t ; 0 ≤ t ≤ T } as follows. ( Step 2) Simulation of the risky asset price process {S i ; i = 1, . . . , d}.
(1) Generate a sample of A k,i according to Step 1 for a fixed k ≥ 1.
(2) With the partition T k at hand, we can apply some appropriate approximation method to evaluate the discounted price. Generally speaking, we work with some Itô-Taylor expansion method. The multidimensional setup requires an additional notation as follows. In the sequel, t k,j denotes the realization of the T k,j by means of Step 1, t k denotes the realization of T k based on the finest random partition T k . Moreover, any sequence (t k 1 < t k 2 < . . . < t k,j 1 ) encodes the information generated by the realization of T k until the first hitting time of the j-th partition. In addition, we denote t k,j 1− as the last time in the finest partition previous to t k,j 1 . Let ν k = (ν 1,k , ν 2,k ) be the pair which realizes Based on this quantities, we define η k t k as the realization of the random variable η k,ν 1,k ν 2,k . Recall expression (3.2). ( Step 3) Simulation of the stochastic derivative D k,j X T k,j 1 .
Based on Steps 1 and 2, for each j = 1, . . . , p one simulates D k,j X T k,j The correspondent simulated pure hedging strategy is given by Remark 5.1. In order to compute the hedging strategy θ k,H hedg over a trading period {s i ; i = 0, . . . , q}, one perform the algorithm described above but based on the shifted filtration and the Brownian motions W (j) si for j = 1, . . . , p as described in Section 4.1.
Remark 5.2. In practice, one has to calibrate the parameters of a given stochastic volatility model based on liquid instruments such as vanilla options and volatility surfaces. With those parameters at hand, the trader must follow the steps (5.1) and (5.4). The hedging strategy is then given by calibration and the computation of the quantity (5.4) over a trading period.

Numerical Analysis and Discussion of the Methods
In this section, we provide a detailed analysis of the numerical scheme proposed in this work. 6.1. Multidimensional Black-Scholes model. At first, we consider the classical multidimensional Black-Scholes model with as many risky stocks as underlying independent random factors to be hedged (d = p). In this case, there is only one equivalent local martingale measure, the hedging strategy θ H is given by (3.6) and the cost is just the option price. To illustrate our method, we study a very special type of exotic option: a BLAC (Basket Lock Active Coupon) down and out barrier option whose payoff is given by It is well-known that for this type of option, there exists a closed formula for the hedging strategy. Moreover, it satisfies the assumptions of Theorem 7.2. See e.g Bernis, Gobet and Kohatsu-Higa [1] for some formulas.
For comparison purposes with Bernis, Gobet and Kohatsu-Higa [1], we consider d = 5 underlying assets, r = 0% for the interest rate and T = 1 year for the maturity time. For each asset, we set initial values S i 0 = 100; 1 ≤ i ≤ 5 and we compute the hedging strategy with respect to the first asset S 1 with discretization level k = 3, 4, 5, 6 and 20000 simulations.
Following the work [1], we consider the volatilities of the assets given by σ 1 = 35%, σ 2 = 35%, σ 3 = 38%, σ 4 = 35% and σ 5 = 40%, the correlation matrix defined by ρ ij = 0, 4 for i = j, where σ i = (σ i1 , · · · , σ i5 ) and we use the barrier level L = 76. Table 1 provides the numerical results based on the algorithm described in Section 5 for the pointwise hedging strategy θ H . Due to Theorem 7.2, we expect that when the discretization level k increases, we obtain results closer to the true value and this is what we find in our Monte Carlo experiments. The standard deviation and percentage % error in Table 1 are related to the average of the hedging strategies calculated by Monte Carlo and the difference between the true and the estimated hedging value, respectively. In Figure 1, we plot the average hedging estimates with respect to the number of simulations. One should notice that when k increases, the standard error also increases, which suggests more simulations for higher values of k.

Hedging Errors.
Next, we present some hedging error results for two well-known non-constant volatility models: The constant elasticity of variance (CEV) model and the classical Heston stochastic volatility model [13]. The typical examples we have in mind are the generalized Föllmer-Schweizer, local risk minimization and mean variance hedging strategies, where the optimal hedging strategies are computed by means of the minimal martingale measure and the variance optimal martingale measure, respectively. We analyze digital and one-touch one-dimensional European-type contingent claims as follows Remark 6.1. When no locally-risk minimizing strategy is available, we also expect to obtain low hedging errors when dealing with generalized Föllmer-Schweizer decompositions due to the orthogonal martingale decomposition. In the mean variance hedging case, two terms appear in the optimal hedging strategy: the pure hedging component θ H,P of the GKW decomposition under the optimal variance martingale measureP andζ as described by (1.4) and (1.5). For the Heston model,ζ was explicitly calculated by Hobson [14]. We have used his formula in our numerical simulations jointly withθ k,H underP in the calculation of the mean variance hedging errors. See expression (6.4) for details.
6.2.1. Constant Elasticity of Variance (CEV) model. The discounted risky asset price process described by the CEV model under the physical measure is given by where B is a P-Brownian motion. The instantaneous sharpe ratio is ψ t = bt−rt σS (β−2)/2 t such that the model can be rewritten as where W is a Q-Brownian motion and Q is the equivalent local martingale measure. For both the digital and one-touch options, we consider the parameters r = 0 for the interest rate, T = 1 (month) for the maturity time, σ = 0.2, S 0 = 100 and β = 1.   Table 3. Hedging error of one-touch option for the CEV model.

Heston's Stochastic Volatility Model.
Here we consider two types of hedging methodologies: Local-risk minimization and mean variance hedging strategies as described in the Introduction and Remark 6.1. The Heston dynamics of the discounted price under the physical measure is given by t ,ρ = 1 − ρ 2 , with (B (1) B (2) ) two independent P-Brownian motions and κ, m, β 0 , µ are suitable constants in order to have a well-defined Markov process (see e.g Heston [13]). Alternatively, we can rewrite the dynamics as Local-Risk Minimization. For comparison purposes with Heath, Platen and Schweizer [12], we consider the hedging of a European put option H written on a Heston model with correlation parameter ρ = 0. We set S 0 = 100, strike price K = 100, T = 1 (month) and we use discretization levels k = 3, 4 and 5. We set the parameters κ = 2.5, θ = 0.04, ρ = 0, σ = 0.3, r = 0 and Y 0 = 0.02. In this case, the hedging strategy θ H,P based on the local-risk-minimization methodology is bounded with continuous paths so that Theorem 7.2 applies to this case. Moreover, as described by Heath, Platen and Schweizer [12], θ H,P can be obtained by a PDE numerical analysis. Table 4 presents the results of the hedging strategyθ k,H 0,0 by using the algorithm described in Section 5. Figure 2 provides the Monte Carlo hedging strategy with respect to the number of simulations of order 10000. We notice that our results agree with the results obtained by Heath, Platen and Schweizer [12] by PDE methods. In this case, the true value of the hedging at time t = 0 is approximately −0.44. The standard errors in Table 4 are related to the hedging and prices computed, respectively, from the Monte Carlo method described in Section 5.  Hedging with generalized Föllmer-Schweizer decomposition for one-touch option. Based on Corollary 4.2, we also present the hedging error associated to one-touch options for a Heston model with non-zero correlation. We simulate the hedging error along the interval [0, 1] using k = 3, 4 as discretization levels and 1 and 2 hedging strategies per day with parameters κ = 3.63, θ = 0.04, ρ = −0.53, σ = 0.3, r = 0, b = 0.01, Y 0 = 0.3 and S 0 = 100 where the barrier is 105. The hedging error result for the one-touch option is summarized in Table 5. The standard deviations in Table 5 are related to the hedging error. To our best knowledge, there is no result about the existence of locally-risk minimizing hedging strategies for one-touch options written on a Heston model with nonzero correlation. As pointed out in Remark 6.1, it is expected that pure hedging strategies based on the generalized Föllmer-Schweizer decomposition mitigate very-well the hedging error. This is what we get in the simulation results.    Table 5. Hedging error with generalized Föllmer-Schweizer decomposition: Onetouch option with Heston model. follows from Remark 6.1. The quantityζ is not related to the GKW decomposition but it is described by Theorem 1.1 in Hobson [14] as follows. The processζ appearing in (1.4) and (1.5) is given by where F is given by (see case 2 of Prop. 5.1 in Hobson [14]) . The initial conditionZ 0 is given byZ The hedging error results are summarized in Table 6 where the standard deviations are related to the hedging error. In comparison with the local-risk minimization methodology, the results show smaller percentual errors when k increases. Also, in all the cases, we had smaller values of the standard deviation which suggests the mean variance methodology provides more accurate values of the hedging strategy.  Table 6. Hedging error in the mean variance hedging methodology for one-touch option with Heston model.

Appendix
This appendix provides a deeper understanding of the Monte Carlo algorithm proposed in this work when the representation (φ H,S , φ H,I ) in (3.6) admits additional integrability and path smoothness assumptions. We present stronger approximations which complement the asymptotic result given in Theorem 3.1. Uniform-type weak and strong pointwise approximations for θ H are presented and they validate the numerical experiments in Tables 1 and 4 in Section 6. At first, we need of some technical lemmas.
Proof. It is sufficient to prove for p = 2 since the argument for p > 2 easily follows from this case. Let H be the linear space constituted by the bounded R 2 -valued F-progressive processes φ = (φ 1 , φ 2 ) such that (7.1) holds with X = X 0 + · 0 φ 1 s dW In order to check (7.1) for such φ, we only need to show for j = 1 since the argument for j = 2 is the same. With a slight abuse of notation, any sub-sigma algebra of F T of the form Ω * 1 ⊗ G will be denoted by G where Ω * 1 is the trivial sigma-algebra on the first copy Ω 1 .
At first, we split Ω = = 0 a.s on the set {T k n−1 ≤ J < T k n = T k,1 1 }. We also have, for some constant C which does not depend on m, k ≥ 1.
Lemma 7.3. Assume that φ H,j ∈ B 2 (F) for some j = 1, . . . , p. Then there exists a constant C such that Proof. By repeating the argument employed in Lemma 7.1 for k ≥ 1, n > 1 and j ∈ {1, . . . , p}, we shall write The following result allows us to get a uniform-type weak convergence of D k,j X under very mild integrability assumption.
Theorem 7.1. Let H be a Q-square integrable contingent claim satisfying assumption (M) and assume that H admits a representation φ H such that φ H,j ∈ B 2 (F) for some j ∈ {1, . . . , p}. Then lim k→∞ D k,j X = φ H,j weakly in B 2 (F).
Proof. Let us fix j = 1, . . . , p. From Lemma 7.3, we know that {D k,j X; k ≥ 1} is bounded in B 2 (F) and therefore this set is weakly relatively compact in B 2 (F). By Eberlein Theorem, we also know that it is B 2 (F)-weakly relatively sequentially compact. From Theorem 3.1, lim k→∞ D k,j X = φ H,j weakly in L 2 (Leb × Q) and since · B 2 is stronger than · L 2 (Leb×Q) , we necessarily have the full convergence lim k→∞ D k,j X = φ H,j in σ(B 2 , M 2 ).
Next, we analyze the pointwise strong convergence for our approximation scheme. Proof. In the sequel, C will be a constant which may differ from line to line and let us fix j = 1, . . . , p. For a given k ≥ 1, it follows from Lemma 7.1 that We recall that T k,j φ H,j uτ j − φ H,j Therefore, the right-hand side of (7.7) vanishes if, and only if, t = 0 is a Lebesgue point of u → ψ H,j (u), i.e., The estimate (7.7), the limit (7.8) and the weak convergence of F k T k,j 1 to the initial sigma-algebra F 0 yield lim k→∞ D k,j X T k,j ; k ≥ 1 then we conclude the proof.
Remark 7.1. At first glance, the limit (7.5) stated in Theorem 7.2 seems to be rather weak since it is not defined in terms of convergence of processes. However, from the purely computational point of view, we shall construct a pointwise Monte Carlo simulation method of the GKW projectors in terms of D k,j X T k,j 1 given by (3.5). This substantially simplifies the algorithm introduced by Leão and Ohashi [20] for the unidimensional case under rather weak path regularity.
Remark 7.2. For each j = 1, . . . , p, let us define ψ H,j (t 0 , u) := E|φ H,j t0+τ j u − φ H,j t0 | 2 , for t 0 ∈ [0, T ], u ≥ 0. One can show by a standard shifting argument based on the Brownian motion strong Markov property that if there exists a representation φ H such that u → ψ H,j (t 0 , u) is cadlag for a given t 0 then one can recover pointwise in L 1 -strong sense the j-th GKW projector for that t 0 . We notice that if φ H,j belongs to B 2 (F) and it has cadlag paths then u → ψ H,j (t 0 , u) is cadlag for each t 0 , but the converse does not hold. Hence the assumption in Theorem 7.2 is rather weak in the sense that it does not imply the existence of a cadlag version of φ H,j .