{'fraction_non_alphanumeric': 0.031765814511961574, 'fraction_numerical': 0.012878580620319989, 'mean_word_length': 4.491431115067884, 'pattern_counts': {'":': 0, '<': 2, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 20, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Improving the Robustness of Trading Strategy Backtesting with Boltzmann Machines and Generative Adversarial Networks

This article explores the use of machine learning models to build a market generator. The underlying idea is to simulate artificial multi-dimensional financial time series, whose statistical properties are the same as those observed in the financial markets. In particular, these synthetic data must preserve the probability distribution of asset returns, the stochastic dependence between the different assets and the autocorrelation across time. The article proposes then a new approach for estimating the probability distribution of backtest statistics. The final objective is to develop a framework for improving the risk management of quantitative investment strategies, in particular in the space of smart beta, factor investing and alternative risk premia.


Introduction
In machine learning, we generally distinguish two types of statistical modeling (Jebara, 2004): • the generative approach models the unconditional probability distribution P (X) given a set of observable variables X; • the discriminative approach models the conditional probability distribution P (Y | X) given a set of observable variables X and a target variable Y .
In finance, we generally observe only one sample path of market prices. For instance, if we would like to build a trading strategy on the S&P 500 index, we can backtest the strategy using the historical values of the S&P 500 index. We can then measure the performance and the risk of this investment strategy by computing annualized return, volatility, Sharpe ratio, maximum drawdown, etc. In this case, it is extremely difficult to assess the robustness of the strategy, since we can consider that the historical sample of the S&P 500 index is one realization of the unknown stochastic process. Therefore, portfolio managers generally split the study period into two subperiods: the 'in-sample' period and the 'out-of-sample' period. The objective is to calibrate the parameters of the trading strategy with one subperiod and measure the financial performance with the other period in order to reduce the overfitting bias. However, if the out-of-sample approach is appealing, it is limited for two main reasons. First, by splitting the study period into two subperiods, the calibration procedure is performed with fewer observations, and does not generally take into account the most recent period. Second, the validation step is done using only one sample path. Again, we observe only one realization of the risk/return statistics. Of course, we could use different splitting methods, but we know that these bootstrap techniques are not well-adapted to times series and must be reserved for modeling random variables. For stochastic processes, statisticians prefer to consider Monte Carlo methods. Nevertheless, financial times series are difficult to model, because they exhibit non-linear autocorrelations, fat tails, heteroscedasticity, regime switching and non-stationary properties (Cont, 2001).
In this article, we are interested in generative models in order to obtain several training/validation sets. The underlying idea is then to generate artificial but realistic financial market prices. The choice of generative model and market data leads naturally to the concept of market generator introduced by Kondratyev and Schwarz (2019). If the market generator is able to replicate the probability distribution of the original market data, we can then backtest quantitative investment strategies on several financial markets. The backtesting procedure is then largely improved, since we obtain a probability distribution of performance and risk statistics, and not only one value. Therefore, we can reduce the in-sample property of the backtest procedure and the overfitting bias of the parameters that define the trading strategy. However, the challenge is not simple, since the market generator must be sufficiently robust and flexible in order to preserve the uni-dimensional statistical properties of the original financial time series, but also the multi-dimensional dependence structure.
This paper is organized as follows. Section Two reviews the more promising generative models that may be useful in finance. In particular, we focus on restricted Boltzmann machines, generative adversarial networks and Wasserstein distance models. In Section Three, we apply these models in the context of trading strategies. We first consider the joint simulation of S&P 500 and VIX indices. We then build an equity/bond risk parity strategy with an investment universe of six futures contracts. Finally, Section Four offers some concluding remarks.

Generative models
In this section, we consider two main approaches. The first one is based on a restricted Boltzmann machine (RBM), while the second one uses a conditional generative adversarial network (GAN). Both are stochastic artificial neural networks that learn the probability distribution of a real data sample. However, the objective function strongly differs. Indeed, RBMs consider a log-likelihood maximization problem, while the framework of GANs corresponds to a minimax two-player game. In the first case, the difficulty lies in the gradient approximation of the log-likelihood function. In the second case, the hard task is to find a learning process algorithm that solves the two-player game. This section also presents an extension of generative adversarial networks, which is called a convolutional Wasserstein model.

Restricted Boltzmann machines
Restricted Boltzmann machines were initially invented under the name 'Harmonium' by Smolensky (1986). Under the framework of undirected graph models 1 , an RBM is a Markov random field (MRF) associated with a bipartite undirected graph 2 . RBMs are made of only two layers as shown in Figure 1. We distinguish visible units belonging to the visible layer from hidden units belonging to the hidden layer. Each unit of visible and hidden layers is respectively associated with visible and hidden random variables. The term 'restricted ' comes from the facts that the connections are only between the visible units and the hidden units and that there is no connection between two different units in the same layer.
We would like to model the distribution of m visible variables V = (V 1 , V 2 , . . . , V m ) representing the observable data whose elements V i are highly dependent. A first way to directly model these dependencies is to introduce a Markov chain or a Bayesian network. In this case, networks are no longer restricted and those methods are computationally expensive particularly when V is a high-dimensional vector. The RBM approach consists in introducing hidden variables H = (H 1 , H 2 , . . . , H n ) as latent variables which will indirectly capture dependencies. Therefore, the hidden layer can be considered as an alternative representation of the visible layer.
1 Fundamental concepts of undirected graph model are explained in Appendix A.1 on page 48. 2 A bipartite graph is a graph whose nodes can be divided into two disjoint and independent sets U and V such that every edge connects a node in U to one in V.

Bernoulli RBMs
Definition A Bernoulli RBM is the standard type of RBMs and has binary-valued visible and hidden variables. Let us denote by v = (v 1 , v 2 , . . . , v m ) and h = (h 1 , h 2 , . . . , h n ) the configurations of visible variables V = (V 1 , V 2 , . . . , V m ) and hidden variables H = (H 1 , H 2 , . . . , H n ), where v i and h j are the binary states of the i th visible variable V i and the j th hidden variable H j such that (v, h) ∈ {0, 1} m+n . The joint probability distribution of a Bernoulli RBM is given by the Boltzmann distribution: (v,h) where the energy function E (v, h) is defined as: where a = (a i ) and b = (b j ) are the two vectors of bias terms associated with V and H and W = (w i,j ) is the matrix of weights associated with the edges between V and H. The normalizing constant Z is the partition function and ensures the overall distribution sums to one 3 . It follows that the marginal probability distributions of the visible and hidden unit states are: (v,h) and: The underlying idea of a Bernoulli RBM is to learn the unconditional probability distributions P (v) of the observable data.
Conditional Distributions According to Long and Servedio (2010), the partition function Z is intractable in the case of Bernoulli RBMs, since its calculus requires summing up 2 m+n elements. Therefore, the probability distribution P (v) is also intractable when m + n increases. However, we can take advantage of the property of the bipartite graph structure of the RBM. Indeed, there are no connections between two different units in the same layer. The probabilities P (v i | h) and P (h j | v) are then independent for all i and j. It follows that: and: With these properties, we can find some useful results that help when computing the gradient of the log-likelihood function on page 5. For instance, we can show that 4 : h P (h | v) h j = P (h j = 1 | v) 3 We have Z = v,h e −E (v,h) . 4 See Appendix A.2.1 on page 51.
A neural network perspective of RBMs In Appendix A.2.2 on page 51, we show that: and: where σ (x) is the sigmoid function: σ (x) = 1 1 + e −x Thus, a Bernoulli RBM can be considered as a stochastic artificial neural network, meaning that the nodes and edges correspond to neurons and synaptic connections. For a given vector v, h is obtained as follows: is a binarizer function and σ is the sigmoid activation function. In addition, we may also go backward in the neural network as follows: According to Fischer and Igel (2014), "an RBM can [then] be reinterpreted as a standard feed-forward neural network with one layer of nonlinear processing units".
Training process Let θ = (a, b, W ) be the set of parameters to estimate. The objective is to find a value of θ such that P θ (v) ≈ P data (v). Since the log-likelihood function (θ | v) of the input vector v is defined as (θ | v) = log P θ (v), the Bernoulli RBM model is trained in order to maximize the log-likelihood function of a training set of N samples v (1) , . . . , v (N ) : where: Hinton (2002) proposed to use gradient ascent method with the following update rule between iteration steps t and t + 1: where η (t) is the learning rate parameter, ∆θ (t) = N s=1 ∇ θ (t) v (s) and ∇ θ (v) is the gradient vector given in Appendix A.2.3 on page 52. Ackley et al. (1985) and Hinton and Sejnowski (1986) showed that the expectation over P (v) can be approximated by Gibbs sampling, which belongs to the family of MCMC algorithms. The goal of Gibbs sampling is to simulate correlated random variables by using a Markov chain. Usually, we initialize the Gibbs sampling with a random vector and the algorithm updates one variable iteratively, based on its conditional distribution given the state of the remaining variables. After a sufficiently large number of sampling steps, we get the unbiased samples from the joint probability distribution. Formally, Gibbs sampling of the joint probability distribution of n random variables X = (X 1 , X 2 , . . . , X n ) consists in sampling

Gibbs sampling
Algorithm 1 Gibbs sampling initialization: Sample x (k) as follows: Let us consider a Bernoulli RBM as a Markov random field defined by the random variables V = (V 1 , V 2 , . . . , V m ) and H = (H 1 , H 2 , . . . , H n ). Since an RBM is a bipartite undirected graph, we can take advantage of conditional independence properties between the variables in the same layer. At each step, we jointly sample the states of all variables in one layer, as shown in Figure 2. Thus, Gibbs sampling for an RBM consists in alternating 5 between sampling a new state h for all hidden units based on P (h | v) and sampling a new state v for all visible units based on denote the states of the visible layer and the hidden layer at time step k. For each unit, we rely on the fact that the conditional probabilities P v | v (k−1) are easily tractable according to Equations (2) and (3). We start by initializing the state of the visible units and we can choose a binary random vector for the first time step. At time step k, here are the steps of the algorithm: 1. We go forward in the network by computing simultaneously for each hidden unit j the following probability p h Contrastive divergence algorithm As shown in the previous paragraph, it is possible to use Gibbs sampling to get samples from the joint distribution P (v, h). However, the computational effort is still too large since the sampling chain needs many sampling steps to obtain unbiased samples. To address this issue, Hinton (2002) initialized the Gibbs sampling with a sample from the real data, instead of using a random vector and suggested the use of only a few sampling steps to get a sample that could produce a sufficiently good approximation of log-likelihood gradient. This faster method is called contrastive divergence algorithm. Given a training set of N samples v (1) , . . . , v (N ) , we obtain: We note P model (v) = P θ (v). Since P data is independent of θ, maximizing the log-likelihood function (5) is equivalent to minimizing the Kullback-Leibler divergence between P data and P model : In contrastive divergence algorithms, we note the distribution of starting values as P (0) , which is also the distribution of real data in training set: P (0) = P data . Let P (k) and P (∞) be the distribution after running k steps of Gibbs sampling and the equilibrium distribution. Compared to P (0) , P (k) is k steps closer to the equilibrium distribution P (∞) , so the divergence measure KL P (0) P (∞) should be greater than or equal to KL P (k) P (∞) . In particular, when the model is well-trained, P (0) , P (k) , and P (∞) should have the same distribution as P data . In this case, KL P (0) P (∞) is equal to 0 and the difference KL P (0) P (∞) − KL P (k) P (∞) is also equal to 0. Thus, according to Hinton (2002), we can find the optimal parameter θ by minimizing the difference KL P (0) P (∞) −KL P (k) P (∞) instead of minimizing directly KL P (0) P (∞) in Equation (7). This quantity is called a contrastive divergence since we compare two KL divergence measures. Therefore, we can rewrite the objective function in Equation (7) as follows: Again, we can use the gradient descent method to find the minimum of the objective function. In this case, the update rule between iteration steps t and t + 1 is: where the gradient vector ∂ θ CD (k) θ (t) is given in Appendix A.2.4 on page 54. Thus, Algorithm (2) summarizes the k-step contrastive divergence algorithm for calibrating θ .

Gaussian-Bernoulli RBMs
A Bernoulli RBM is limited to modelling the probability distribution P (v) where v is a binary vector. To address this issue, given an RBM with m visible units and n hidden units, we can associate a normally distributed variable to each visible unit and a binary variable to each hidden unit. This type of RBMs is called Gaussian-Bernoulli RBM. We are free to choose the expression of the energy function as long as it satisfies the Hammersley-Clifford theorem 6 and its partition function is well defined. For instance, Cho et al. (2011) defined the energy function of this RBM as follows: where a i and b j are bias terms associated with visible variables V i and hidden variables H j , w i,j is the weight associated with the edge between V i and H j , and σ i is a new parameter associated with V i . Following the same calculus in Equation 3, we can show that the conditional probability is equal to: Moreover, according to Krizhevsky (2009), the conditional distribution of V i given h is Gaussian and we obtain: As we have previously seen, we can maximize the log-likelihood function in order to train a Gaussian-Bernoulli RBM with m visible units and n hidden units. As all visible units are Algorithm 2 k-step contrastive divergence algorithm input: v (1) , . . . , v (N ) initialization: We note ∆a i = ∂ ai CD (k) , ∆b j = ∂ bi CD (k) and ∆w i,j = ∂ wi,j CD (k) , and we set ∆a i = ∆b j = ∆w i,j = 0 end for end for end for Output: gradient estimation ∆a i , ∆b j and ∆w i,j associated with continuous probability distribution, the log-likelihood function (θ | v) of an . . , V m ). In Appendix A.2.5 on page 55, we give the expression of the gradient vector. Therefore, there is no difficulty to use the gradient ascent method or the k-step contrastive divergence algorithm to train a Gaussian-Bernoulli RBM.
Remark 1. If we normalize the data of the training set using a z-score function, we can set the standard deviation σ i to 1 during the training process. This reduces the number of parameter and accelerates the convergence of the training process.

Conditional RBM structure
Traditional Bernoulli and Gaussian-Bernoulli RBMs can only model the static dependence of variables. However, in the practice of financial data modeling, we also want to capture the temporal dependencies between variables. To address this issue, Taylor et al. (2011) introduced a conditional RBM structure by adding a new layer to the Gaussian-Bernoulli RBM. Thus, this conditional RBM is made of three parts as shown in Figure 3: 2. a visible layer with m units; 3. a conditional layer with d conditional meta-units {C 1 , . . . , C d }. Since we model the temporal structure, we note the observation at time t as v t = (v t,1 , v t,2 , . . . , v t,m ). The conditional layer is fed with d past values (v t−1 , . . . , v t−d ) that are concatenated into a d × m-dimensional vector 7 and we note it as c t . The conditional layer is fully linked to visible and hidden layers with directed connections. Let us denote by Q = (q d×m,m ) the weight matrix connecting the conditional layer to the visible layer and P = (p d×m,n ) the weight matrix connecting the conditional layer to the hidden layer.
Hidden layer (n units)

Visible layer (m units)
A conditional RBM contains both undirected and directed connections in the graph. Thus, it can't be defined as an MRF or Bayesian network. However, conditionally on c t , we can consider both visible and hidden layers as an undirected graph and we can compute P (v t , h t | c t ) instead of computing P (v t , h t , c t ) in order to take still advantage of undirected graph properties. Thus, according to the Hammersley-Clifford theorem, the conditional probability distribution has the following form: vt,i,ht,j ,ct) 7 A conditional meta-unit is then composed of m values: where Z is the partition function. Taylor et al. (2011) proposed a form of energy function that is an extension of the Gaussian-Bernoulli RBM energy: are dynamic bias terms with respect to a and b. Thus, the energy functionẼ g corresponds to a Gaussian-Bernoulli RBM energy function E g by replacing constant biases a and b by dynamic biasã t andb t . By updating these terms in Equations (9) and (10), we obtain: and: Again, the partition function Z is still intractable: ht vt e −Ẽg(vt,ht,ct ) dv t Remark 2. As for Gaussian-Bernoulli RBMs, we can use the k-step contrastive divergence algorithm for the training process with some modifications of the gradient vector 8 .

Generative adversarial networks
Generative models are an important part of machine learning algorithms that learn the underlying probability distributions of the real data sample. In other words, given a finite sample data with a distribution P data (x), can we build a model such that P model (x; θ) ≈ P data (x)? The goal is to learn to sample a complex distribution given a real sample. Generative adversarial networks (GANs) belong to the class of generative models and move away from the classical likelihood maximization approach, whose objective is to estimate the parameter θ. GAN models enable the estimation of the high dimensional underlying statistical structure of real data and simulate synthetic data, whose probability distribution is close to the real one. To assess the difference between real and simulated data, GANs are trained using a discrepancy measure. Two widely classes of discrepancy measures are information-theoretic divergences and integral probability metrics. The choice of the GAN objective function associated with the selected discrepancy measure explains the multitude of GAN models appearing in the machine learning literature. However, the different GAN models share a common framework. Indeed, Appendix A.3 page 58 shows how the original formulation of GAN models is a particular case of the theory of φ-divergence and probability functional descent. Moreover, the link between influence function used in robust statistics and the formulation of the discriminator and generator is derived.
According to Goodfellow et al. (2014), "a generative adversarial process trains two models: a generative model G that captures the data distribution, and a discriminative model D that estimates the probability that a sample came from the training data rather than G. The training procedure for G is to maximize the probability of D making a mistake. This framework corresponds to a minimax two-player game. In the space of arbitrary functions G and D, a unique solution exists, with G recovering the training data distribution and D equal to 1 /2 everywhere". Therefore, a generative adversarial model consists of two neural networks. The first neural network simulates a new data sample, and is interpreted as a data generation process. The second neural network is a classifier. The input data are both the real and simulated samples. If the generative model has done a good job, the discriminative model is unable to know whether an observation comes from the real dataset or the simulated dataset.

Adversarial training problem
We use the framework of Goodfellow et al. (2014) and Wiese et al. (2020). Let Z be a random noise space. z ∈ Z is sampled from the prior distribution P noize (z). The generative model G is then specified as follows: where X denotes the data space and Θ G defines the parameter space including generator weights and bias that will be optimized. The generator G helps to simulate data x 0 . The discriminative model is defined as follows: where x corresponds to the set of simulated data x 0 and training (or real) data x 1 . In this approach, the statistical model of X corresponds to: The probability that the observation comes from the real (or true) data is equal to p 1 = D (x 1 ; θ d ), whereas the probability that the observation is simulated is given by p 0 = D (x 0 ; θ d ). If the model (15) is wrong and does not reproduce the statistical properties of the real data, the classifier has no difficulty in separating the simulated data from the real data, and we obtain p 0 ≈ 0 and p 1 ≈ 1. Otherwise, if the model is valid, we must verify that: The main issue of GANs is the specification of the two functions G and D and the estimation of the parameters θ g and θ d associated to G and D. For the first step, Goodfellow et al. (2014) proposed to use two multi-layer neural networks G and D, whereas they consider the following cost function for the second step: The optimization problem becomes: In other words, the discriminator is trained in order to maximize the probability to correctly classify historical samples from simulated samples. The objective is to obtain a good classification model since the maximum value C (θ g , θ d ) with respect to θ d is reached when: In the meantime, the generator is trained in order to minimize the probability that the discriminator is able to perform a correct classification or equivalently to maximize the probability to fool the discriminator, since the minimum value C (θ g , θ d ) with respect to θ g is reached when: Remark 3. In Appendix A.5 on page 67, we show that the cost function is related to the binary cross-entropy measure or the opposite of the log-likelihood function of the logit model. Moreover, the cost function can be interpreted as a φ-divergence measure as explained in Appendix A.3.1 on page 58.

Solving the optimization problem
The minimax optimization problem is difficult to solve directly, because the gradient vector ∇ θ C (θ g , θ d ) is not well informative if the discriminative model is poor. Therefore, the traditional way to solve this problem is to use a two-stage approach: 1. In a first stage, the vector of parameters θ g is considered to be constant whereas the vector of parameters θ d is unknown. This implies that the minimax problem reduces to a maximization step:θ where C (max) (θ d | θ g ) corresponds to the cost function by assuming that θ g is given.
2. In a second stage, the vector of parameters θ d is considered to be constant whereas the vector of parameters θ g is unknown. This implies that the minimax problem reduces to a minimization step:θ where C (min) (θ g | θ d ) corresponds to the cost function by assuming that θ d is given.
The two-stage approach is repeated until convergence by setting θ g to the valueθ (min) g calculated at the minimization step and θ d to the valueθ (max) d calculated at the maximization step. The cycle sequence is given in Figure 4. A drawback of this approach is the computational time. Indeed, this implies to solve two optimization problems at each iteration.
Therefore, Goodfellow et al. (2014) proposed another two-stage approach, which converges more rapidly. The underlying idea is not to estimate the optimal generator at each iteration. The objective is, rather, to improve the generative model at each iteration, such that the new estimated model is better than the previous. The convergence is only needed for the discriminator step (Goodfellow et al., 2014, Proposition 2). Moreover, the authors applied a mini-batch sampling in order to reduce the computational time. It follows that the cost functions become: and: where m is the size of the mini-batch sample and c is a constant that does not depend on the generator parameters θ g : In Algorithm (3), we describe the stochastic gradient optimization method proposed by Goodfellow et al. (2014) when the learning rule corresponds to the steepest descent method. But other learning rules can be used such as the momentum method or the adaptive learning method. In practice, these algorithms may not converge 9 . Another big issue known as 'mode collapse ' concerns the diversity of generated samples. In this case, the simulated data drawn by the generator exhibit small and limited differences, meaning all generated samples are almost identical. This is why researchers have proposed alternatives to Algorithm (3) in order to solve these two problems (Mao et al., 2017;Karras et al., 2018;Metz et al., 2017).

Time series modeling with GANs
It is possible to generate fake time series from a single random noise vector z without any labels. However, this implies the use of complex structures for both generator and

Algorithm 3 Stochastic gradient optimization for GAN training
The goal is to compute the optimal parametersθ g andθ d We note n max the number of iterations to apply to the maximization problem (discriminator step) and n min the number of iterations to apply to the minimization problem (generator step) m is the size of the mini-batch sample The starting values are denoted by θ discriminator. These structures can cause the models to be computationally expensive, particularly when using convolutional networks. Moreover, such structures don't address the issue concerning the lack of data. Let us remember that we only have one single historical scenario. A less expensive alternative is to use labels combined with simple multi-layer perceptrons for the generator and the discriminator. This approach has been introduced by Mirza and Osindero (2014) and is called conditional generative adversarial networks (cGANs).
Before training a GAN, the data need to be labelled. An additional vector that encodes structural information about real historical data must then be defined. It will help the GAN to generate specific scenarios depending on labels defined before. Therefore, the learning process will be supervised. The previous framework remains valid, but the generative and discriminative models are written as G z | v; θ g and D (x| v; θ d ), where v ∈ V is the label vector. Therefore, the cost function becomes: The label vector may encode various types of information, and can be categorical or continuous. For example, if we consider the S&P 500 index, we can specify V = {−1, 0, +1} depending on its short trend. Let y t be the value of the S&P 500 index and m t the corresponding 10-day moving average. v t is equal to −1 if the S&P 500 index exhibits a negative trend, 0 if it has no trend, and +1 otherwise 10 . An example of continuous labels is a vector composed of the last p values (Koshiyama et al., 2019). This type of labels acts as a time memory and helps to reproduce auto-correlation patterns of the stochastic process.

Wasserstein GAN models
In Section 3, we will see that the basic GAN model using the cross-entropy as the loss function for the discriminator suffers from three main problems. First, the visualization of the training process is not obvious. Traditionally, we focus on the loss error curve to decide whether or not the network is trained correctly. In the case of financial applications, looking at binary cross entropy is not relevant. Computer vision algorithms can trust generated images to evaluate the GAN, financial time series are extremely noisy. Therefore, a visual evaluation is much more questionable. Second, the basic GAN model is not suitable for generating multi-dimensional time series. A possible alternative is to modify the GAN's structure. For instance, we can replace the binary cross entropy function by the Wasserstein distance for the training error (Arjovsky et al., 2017a,b), or temporal information of multidimensional time series can be encoded using a more complex structure such as convolutional neural networks  or recurrent neural networks (Hyland et al., 2017). Third, the mode collapse phenomenon must be addressed in order to manage the poor diversity of generated samples, because we would like to have several simulated time series, as we have when performing Monte Carlo methods.

Optimization problem
Let P and Q be two univariate probability distributions. The Wasserstein (or Kantorovich) distance between P and Q is defined as: where F (P, Q) denotes the Fréchet class of all joint distributions F (x, y) whose marginals are equal to P and Q. In the case p = 1, the Wasserstein distance represents the cost of the optimal transport problem 11 (Villani, 2008): because F (x, y) describes how many masses are needed in order to transport one distribution to another (Arjovsky et al., 2017a,b). The Kantorovich-Rubinstein duality introduced by Villani (2008) allows us to rewrite Equation (17) as follows: In the case of a generative model, we would like to check whether sample and generated data follows the same distribution. Therefore, we obtain P = P data and Q = P model . In the special case of a GAN model, P model (x; θ g ) is given by G (P noise (z) ; θ g ). Arjovsky et al. (2017a,b) demonstrated that if G (z; θ g ) is continuous with respect to θ, then W (P data , P model ) is continuous everywhere and differentiable almost everywhere. This result implies that GAN can be trained until it has converged to its optimal solution contrary to basic GANs, where the training loss were bounded and not continuous 12 . Moreover, Gulrajani et al. (2017) adapted Equation (18) in order to recover the two player min-max game: The Wasserstein GAN (WGAN) plays the same min-max optimization problem as previously:θ This implies that the discriminator needs to be a 1-Lipschitz function.
Remark 4. The discriminator function is not necessarily a classifier, since the output of the function D can take a value in [0, 1]. It can then be a general scoring method.
Remark 5. The Kantorovich-Rubinstein duality makes the link between Wasserstein distance and integral probability metrics presented in Appendix A.3.2 on page 63.

Properties of the optimal solution
In Wasserstein GAN models, the cost function becomes: 11 See Appendix A.6 on page 68 for an introduction of Monge-Kantorovich problems, and the relationship between optimal transport and the Wasserstein distance.
12 There is also no problem in computing the gradient. and we have: We recall that the discriminator ϕ needs to satisfy the Lipschitz property otherwise the loss gradient will explode 13 . A first alternative proposed by Arjovsky et al. (2017a,b) is to clip the weights into a closed space [−c, c]. However, this method is limited because the gradient can vanish and weights can saturate. A second alternative proposed by Gulrajani et al. (2017) is to focus on the properties of the optimal discriminator gradient. They showed that the optimal solution ϕ has a gradient norm 1 almost everywhere under P and Q. Therefore, they proposed to add a regularization term to the cost function in order to constraint the gradient norm to converge to 1. This gradient penalty leads to define a new cost function: where P mixture is a mixture distribution of P and Q. More precisely, Gulrajani et al. (2017) proposed to sample x 2 as follows: Remark 6. Gulrajani et al. (2017) found that a good value of the coefficient λ is 10.

Convolutional neural networks
A convolutional neural network (CNN) is a class of deep neural networks, where the inputs are transformed using convolution and filtering operators. For instance, this type of neural networks is extensively used in computer vision. Recently, Radford et al. (2016) implemented a new version of GANs using convolutional networks as generator and discriminator in order to improved image generation 14 . In this approach, the underlying idea is to extract pertinent features. In finance, inputs are different and correspond to n x -dimensional time series representing asset prices over n t days. Therefore, inputs are represented by a matrix belonging to M nt,nx (R), and are characterized by the time axis, where order matters and the asset axis, where order doesn't matter. Therefore, the input transformation is referring to 1-dimensional convolution (1D-CNN).

Extracting features using convolution
Convolution is no more than a linear transformation of a given input. But, contrary to simple multilayer perceptrons, convolutional operations preserve the notion of ordering according to a specific axis. To achieve this, a weight matrix of a given length n k called kernel will slide across the input data. At each location, a continuous part of the input data is selected, and the kernel receives this vector as input in order to produce a single output by matrix product. This process is repeated with n f different filters of similar dimension. Consequently, we obtained n f down-sampled versions of the original input. The way the kernel can slide across the input data is controlled by the stride n s . It is defined as the distance between two consecutive positions of the kernel (Dumoulin and Visin, 2016). The higher the stride value, the more important the sub-sampling. If s is equal to one, all the data are considered and there is no sub-sampling. Finally, the notion of padding allows us to address the situation when the kernel arrive at the end or the beginning of an axis. So, padding is defined as the number of zeros concatenated at the beginning or at the end of a given axis. In finance, we want to down-sample a given multi-dimensional time series belonging to M nt,nx (R) into a subspace M n t ,nx (R) with the condition that n t < n t . With a 1dimensional convolution, output time axis is automatically defined with respect to the dimension of the kernel. The first dimension of the kernel is free, whereas the second dimension is set to the number of time series n x . Thus, the kernel is defined by a weight matrix belonging to M n k ,nx (R). In order to sample the input data, the stride is chosen, such that s ≥ 1. The padding is set, such that n p rows of zeros is padded at the (bottom and top) limits of the input data. Finally, the output will belong to M n t ,nx (R), such that: This type of convolution will be used to build the discriminator that should output a single scalar. In this case, down-sampling real or fake time series become essential.

Up-sampling a feature using transpose convolution
This transformation also called 'deconvolution' is useful to build the generator when we would like to up-sample a random noise in order to produce a fake time series belonging to M nt,nx (R). More generally, it could be used as the decoding layer in all forms of autoencoders. Transpose convolution is the exact inverse transformation of the convolution previously defined. Considering an input data belonging to M nt,nx (R), we obtain an output data belonging to M n t ,nx (R) such that: The term transpose comes from the fact that convolution is in fact a matrix operation. When we perform a convolution, input matrix flattens into a n t n x -dimensional vector. All the parameters (stride, kernel settings or padding) are encoded into the weight convolution matrix belonging to M ntnx,n t (R). Therefore, the output data is obtained by computing the product between the input vector and the convolution matrix. Performing a transpose convolution is equivalent of taking the transpose of this convolution matrix.

Financial applications
3.1 Application of RBMs to generate synthetic financial time series A typical financial time series of length T may be described by a real-valued vector y = (y 1 , ..., y T ). As we have introduced in Sections 2.1.1 and 2.1.2 on page 4, Bernoulli RBMs take binary vectors as input for training and Gaussian RBMs take the data with unit variance as input in order to ignore the parameter σ. Therefore, data preprocessing is necessary before any training process of RBMs. For a Bernoulli RBM, each value of y t needs to be converted to an N -digit binary vector using the algorithm proposed by Kondratyev and Schwarz (2019). The underlying idea consists in discretizing the distribution of the training dataset and representing them with binary numbers. For instance, N may be set to 16 and a Bernoulli RBM, which takes a one-dimensional time series as training dataset, will have 16 visible layer units. In the same way, the visible layer will have 16 × n neurons in the case of an n-dimensional time series. Moreover, samples generated by a Bernoulli RBM after Gibbs sampling are also binary vectors and we need another algorithm to transform binary vectors into real values. These transformation algorithms between real values and binary vectors are detailed in Appendix A.7 on page 72. For a Gaussian RBM, we need only to normalize data to unit variance and scale generated samples after Gibbs sampling.
For the training process of RBMs, we use the contrastive divergence algorithm CD (1) to estimate the log-likelihood gradient and the learning rate η (t) is set to 0.01. All models are trained using mini-batch gradient descent with batch size 500 and we apply 100 000 epochs to ensure that models are well-trained. Using a larger mini-batch size will get a more stable gradient estimate at each step, but it will also use more memory and take a longer time to train the model.
After having trained the RBMs, we may use these models to generate synthetic (or simulated) samples to match training samples. Theoretically, a well-trained RBM should be able to transform random noise into samples that have the same distribution as the training dataset after doing Gibbs sampling for a sufficiently long time. Therefore, the trained RBM is fed by samples drawn from the Gaussian distribution N (0, 1) as input data. We then perform a large number of forward and backward passes through the network to ensure convergence between the probability distribution of generated samples and the probability distribution of the training sample. In the particular case of Bernoulli RBMs, after the last backward pass from the hidden layer to the visible layer, we need to transform generated binary vectors into real-valued vectors.

Simulating multi-dimensional datasets
In this first study, we test Bernoulli and Gaussian RBMs on a simulated multi-dimensional dataset to check whether RBMs can learn the joint distribution of training samples that involves marginal distributions and a correlation structure. We simulate 10 000 observations of 4-dimensional data with different marginal distributions, as shown in Figure 6. The first dimension consists of samples drawn from a Gaussian mixture model 15 . The samples in the second dimension are drawn from a Student's t distribution with 4 degrees of freedom. For the third and fourth dimensions, we draw respectively samples from the Gaussian distribution N (0, 2). The empirical probability distributions of these 4-dimensional data are reported in Figure 6. In addition, we use a Gaussian copula to construct a simple correlation structure of the simulated samples with the correlation matrix of the Gaussian copula

100%
-60% 0% 30% given in Figure 7. We set a strong negative correlation −60% between the Gaussian mixture model and the Student's t distribution and a strong positive correlation 50% between the two Gaussian distributions. Moreover, the Gaussian distribution in the fourth dimension has a more complicated correlation structure 16 than the Gaussian distribution in the third dimension, which is independent from the Gaussian mixture model and the Student's t distribution or the first and second dimensions. Here, the challenge is to learn both the marginal distribution and the copula function.
Before implementing the training process, we need to update the parameters for the RBMs to adjust multi-dimensional input. In the case of the Bernoulli RBM, each value should be transformed into a 16-digit binary vector and we need to concatenate them to form a 64-digit binary vector. So, the visible layer of the Bernoulli RBM has 64 neurons. For the Gaussian RBM, we have simply 4 visible layer units. According to our experience, we need to set a large number of hidden layer units in order to learn marginal distributions and the correlation structure at the same time. Therefore, we choose 256 hidden layer units for the Bernoulli RBM and 128 hidden layer units for the Gaussian RBM. We also recall that the number of epochs is set to 100 000 in order to ensure that models are well-trained. After the training process, we perform 1 000 steps of the Gibbs sampling on a 4-dimensional random noise to generate 10 000 samples with the same size as the training dataset. These simulated samples are expected to have not only the same marginal distributions but also the same correlation structure as the training dataset. Figures 8 and 9 compare the histograms and QQ-plots between training samples and generated samples after 1 000 Gibbs sampling steps in the case of the Bernoulli RBM. According to these figures, we observe that the Bernoulli RBM can learn pretty well each marginal distribution of training samples. However, we also find that the Bernoulli RBM focuses on extreme values in the training dataset instead of the whole tail of the distribution and this phenomenon is more evident for heavy-tailed distributions. For example, in the case of the Student's t distribution (Panel (b) in Figure 9), the Bernoulli RBM tries to learn the behavior of extreme values, but ignores the other part of the distribution tails. Comparing the results for the two Gaussian distributions (Panels (c) and (d) in Figure 9), we notice that the learning of the Gaussian distribution in the fourth dimension, which has a more complicated correlation structure with the other dimensions, is not as good as the result for the Gaussian distribution in the third dimension.

Bernoulli RBM
We have replicated the previous simulation 50 times 17 and we compare the average of the mean, the standard deviation, the 1 st percentile and the 99 th percentile between training and simulated samples in Table 1. Moreover, in the case of simulated samples, we indicate the confidence interval by ± one standard deviation. These statistics demonstrate that the Bernoulli RBM overestimates the value of standard deviation, the 1 st percentile and the 99 th percentile, especially in the case of Student's t distribution. This means that the Bernoulli RBM is sensitive to extreme values and the learned probability distribution is more heavy-tailed than the empirical probability distribution of the training dataset.
In Figure 10, we show the comparison between the empirical correlation matrix 18 of the training dataset and the average of the correlation matrices computed by the 50 Monte Carlo replications. We notice that the Bernoulli RBM does not capture perfectly the correlation structure since the correlation coefficient values are less significant comparing with the empirical correlation matrix. For instance, the correlation between the first and second dimensions is equal to −57% for the training data, but only −31% for the simulated data on average 19 .

Generated samples
Gaussian RBM Let us now consider a Gaussian RBM with 4 visible layer units and 128 hidden layer units using the same simulated training dataset. After the training process, synthetic samples are always generated using the Gibbs sampling with 1 000 steps between visible and hidden layers. Histograms and QQ-plots of training and generated samples are reported in Figures 11 and 12. We notice that the Gaussian RBM works well for the Gaussian mixture model and the two Gaussian distributions (Panels (a), (c) and (d) in Figure 12). Again, we observe that the Gaussian distribution with the simplest correlation structure is easier to learn than the Gaussian distribution with a more complicated correlation structure. However, the Gaussian RBM fails to learn heavy-tailed distributions such as the Student's t distribution since the model tends to ignore all values in the distribution tails.
Comparing Tables 1 on page 24 for the Bernoulli RBM and Tables 2 on page 26 for the Gaussian RBM, we notice that the Gaussian RBM generally underestimates the value of the standard deviation, the 1 st percentile and the 99 th percentile. This means that this is challenging to generate leptokurtic probability distributions with Gaussian RBMs.
In Figure 13, we also compare the empirical correlation matrix of the training dataset Improving the Robustness of Backtesting with Generative Models

Generated samples
Summary of the results According to our tests, we consider that both Bernoulli and Gaussian well-trained RBMs can transform random noise series to the joint distribution of a training dataset. But they have their own characteristics: • A Bernoulli RBM is very sensitive to extreme values of the training dataset and has a tendency to overestimate the tail of distribution. On the contrary, a Gaussian RBM has a tendency to underestimate the tails of probability distribution and faces difficulties in learning leptokurtic probability distributions.
• Gaussian RBMs may capture more accurately the correlation structure of the training dataset than Bernoulli RBMs.
In practice, we may apply the normal score transformation in order to transform the training dataset, such that each marginal has a standard normal distribution. For that, we rank the values of each dimension from the lowest to the highest, map these ranks to a uniform distribution and apply the inverse Gaussian cumulative distribution function. We then train the RBM with these transformed values. After the training process, the marginals of samples generated by well-trained RBMs will follow a standard normal distribution. By using the inverse transformation, we can generate synthetic samples with the same marginal distributions than those of the training dataset, and a correlation structure that is closed to the empirical correlation matrix of the training data. In this case, we may consider Gaussian RBMs as an alternate method of the bootstrap sampling approach.

Application to financial time series
According to the results in the previous section, Gaussian RBMs capture more accurately the correlation structure of the training dataset and we can avoid its drawbacks by applying the normal score transformation. In this paragraph, we test and compare a Gaussian RBM and a conditional RBM on a multi-dimensional autocorrelated time series in order to check if conditional RBMs may capture at the same time the correlation structure and the time dependence of the training dataset. The RBM is trained on a real financial dataset consisting of 5 354 historical daily returns of the S&P 500 index and the VIX index from December 1998 to May 2018 (see Figure 14). During this period, the S&P 500 index and the VIX index have a remarkable negative correlation that is equal to −71%. In order to reinforce the autocorrelation level in the training samples, we apply a 3-day exponential weighted moving average approach to historical prices before calculating daily returns. Figure 15 shows the autocorrelation of training data. As both daily returns of the S&P 500 and VIX indices are leptokurtic and have some extreme values, we need to apply the normal score transformation on input data before the training process.

VIX index
Normal score transformation In order to verify the learning quality of Gaussian RBMs using the normal score transformation, we run 50 Monte Carlo simulations after the training process and for each simulation, we generate 5 354 simulated observations from different Autocorrelation VIX SPX random noise series. We then calculate the statistics of these samples as we have done in the previous section. According to Table 3, we find that although the Gaussian RBM still underestimates a little the standard deviation, the 1 st percentile and the 99 th percentile, the learning quality is much more improved. In addition, the average correlation coefficient between daily returns of S&P 500 and VIX indices over 50 Monte Carlo simulations is equal to −76%, which is very closed to the figure 71% of the empirical correlation. Therefore, we consider that training Gaussian RBMs with the normal score transformation may be used as an alternate method of bootstrap sampling. Learning time dependence As conditional RBM is an extension of Gaussian RBM, applying the normal score transformation should still work. Moreover, conditional RBM should be able to learn the time-dependent relationship by its design. So, we train a conditional RBM on the same training dataset in order to check whether this model can capture the autocorrelation patterns previously given in Figure 15.
In practice, we choose to use a long memory for the conditional layer in the conditional RBM and the values of the last 20 days will be fed to the model. Therefore, the conditional RBM has 2 visible layer units, 128 hidden layer units and 40 conditional layer units 20 . We then use the normal score transformation in order to transform the training dataset as mentioned before. After having done a training process of 100 000 epochs, we are interested in generating samples of consecutive time steps to verify if the conditional RBM can learn the joint distribution but also capture the autocorrelation of the training dataset. We run 3 Monte Carlo simulations and for each simulation, we generate consecutively samples of 250 observations. In Figures 16 and 17, we compare the autocorrelation function of generated samples by Gaussian and conditional RBMs. These results show clearly that a conditional RBM can capture well the autocorrelation of the training dataset, which is not the case of a traditional Gaussian RBM.

VIX index
Market generator Based on above results, we consider that conditional RBMs can be used as a market generator (Kondratyev and Schwarz, 2019). In the example of S&P 500 and VIX indices, we have trained the models with all historical observations and for each date, we have used the values of the last 20 days as a memory vector for the conditional layer. In this approach, the normal score transformation ensures the learning of marginal distributions, whereas the conditional RBM ensures the learning of the correlation structure and the time dependence. After having calibrated the training process, we may choose a date as the starting point. The values of this date and its last 20 days are then passed to the trained RBM model and, after sufficient steps of Gibbs sampling, we get a sample for the next day. We then feed this new sample and the updated memory vector to the model in order to generate a sample for the following day. Iteratively, the conditional RBM can generate a multi-dimensional time series that has the same patterns as those of the training data in terms of marginal distribution, correlation structure and time dependence.
For instance, Figure 18 shows three time series of 250 trading days starting at the date 22/11/2000, which are generated by the trained conditional RBM. We notice that these three simulated time series have the same characteristics as the historical prices represented by the dashed line. The negative correlation between S&P 500 and VIX indices is also learned and we can see clearly that there exists a positive autocorrelation in each simulated one-dimensional time series.

Application of GANs to generate fake financial time series
In this section, we perform the same tests for GAN models as those we have done for RBMs. The objective is to compare the results between these two types of generative models on the same training dataset. We recall that the first test consists of learning the joint distribution of a simulated multi-dimensional dataset with different marginal distributions and a simple correlation structure driven by a Gaussian copula. In this part, we test a Wasserstein GAN model with only simple dense layers in the generator and in the discriminator. In the second test, we train a conditional Wasserstein GAN model on the historical daily returns of the S&P 500 index and the VIX index from December 1998 to May 2018. As we have done in the previous tests for RBMs, we also apply a three-day exponential weighted moving average approach to historical prices before calculating daily returns. We know that this transformation will reinforce the autocorrelation level in the training samples. As a result, we need to choose a more complex structure for the generator and the discriminator to capture the complex features in the training dataset. For instance, we construct the generator and the discriminator with convolutional neural network as introduced in Section 2.4 on page 18. According to the results of these two tests, we check whether the family of GAN models may also learn the marginal distribution, the correlation structure and the time dependence of the training dataset as RBMs can do.

Simulating multi-dimensional datasets
As we have introduced in Section 2.3 on page 16, Wasserstein GAN models, which use the Wasserstein distance as the loss function, have several advantages comparing the basic traditional GAN models using the cross-entropy loss measure. Therefore, we test in this first study Wasserstein GAN models on the simulated multi-dimensional dataset as we have done for Bernoulli and Gaussian RBMs. We recall that we simulate 10 000 samples of 4dimensional data with different marginal distributions (Gaussian mixture model, Student's t distribution and two Gaussian distributions), as shown in Figure 6 on page 21 and the training dataset has a correlation structure simulated by a Gaussian copula with the correlation matrix given in Figure 7 on page 21. Here, the objective of the test is to check whether a well-trained Wasserstein GAN models can learn the marginal distribution and the copula function of training samples.
Data preprocessing In the practice of GAN models, we need to match the dimension and the range of the output of the generator and the input of discriminator. To address this issue, there are two possible ways: • We can use a complex activation function for the last layer of the generator in order to ensure that the generated samples have the same range as the training samples for the discriminator (Clevert et al., 2015).
• We can modify the range of the training samples by applying a preprocessing function and choose a usual activation function (Wiese et al., 2020).
In our study, we consider the second approach by applying the MinMax scaling function 21 to the training samples for the discriminator. As a result, the input data scaled for the discriminator will take their values in [0, 1] and we may choose the sigmoid activation function for the last layer of the generator to ensure that the outputs of the generator have theirs values in [0, 1]. In addition, we need to simulate the random noise as the inputs for the generator. As the random noise will play the role of latent variables, we have the freedom to choose its distribution. For instance, we use the Gaussian distribution z ∼ N (0, 1) in this study.
Structure and training of Wasserstein GAN models As we have mentioned, we need to match the output of the generator and the input of the discriminator. Therefore, the last layer in the generator should have 4 neurons for our simulated training dataset and we may use sigmoid activation function for each neuron as we have applied MinMax scaling function to the training samples. The discriminator corresponds to a traditional classification problem, so its last layer should have only one neuron. We may choose the sigmoid activation function for this neuron if the training samples are labelled by {0, 1} or the tangent hyperbolic activation function if the training samples are labelled by {−1, 1}.
In our study, we want to construct a Wasserstein GAN model, which can convert a 100-dimensional random noise vector to the data with the same joint distribution as the training dataset. According to our tests, a simple structure as multi-layer perceptrons for the generator and the discriminator is sufficient for learning the distribution of these simulated training samples. More precisely, the generator is fed by 100-dimensional vectors and has 5 layers with the structure S G = {200, 100, 50, 25, 4} where s i ∈ S G represents the numbers of neurons of the i th layer. The discriminator has 4 layers with S D = {100, 50, 10, 1} and takes 4-dimensional vectors as input data. As explained in the previous paragraph, we may choose the sigmoid activation function for the last layer of the generator and the tangent hyperbolic activation function for the last layer of the discriminator. For all other layers of the discriminator and the generator, the activation function f (x) corresponds to a leaky rectified linear unit (or RELU) function with α = 0.5: According to Arjovsky et al. (2017a), the Wasserstein GAN model should be trained with the RMSProp optimizer using a small learning rate. Therefore, we set the learning rate to η = 1.10 −4 and use mini-batch gradient descent with batch size 500. Table 4 summarizes the setting of the Wasserstein GAN model implementation using Python and the TensorFlow library. After the training process, we generate 10 000 samples in order to compare the marginal distributions and the correlation structure with the training dataset. Results Figures 19 and 20 compare the histograms and QQ-plots between training samples and generated samples. According to these figures, we observe that the Wasserstein GAN model can learn very well each marginal distribution of training samples, even better than Bernoulli and Gaussian RBMs. In particular, the Wasserstein GAN model fits more accuracy for heavy-tailed distributions than RBMs as in the case of the Student's t distribution (Panel (b) in Figure 20). As we have done for RBMs, we replicate the previous simulation 50 times and we compare the average of the mean, the standard deviation, the 1 st percentile and the 99 th percentile between training and simulated samples in Table 5. Moreover, in the case of simulated samples, we indicate the confidence interval by ± one standard deviation. These statistics show that the Wasserstein GAN model estimates very well the value of standard deviation, the 1 st percentile and the 99 th percentile for each marginal distribution, especially for Gaussian mixture model and Student's t distribution in the first and second dimensions. One shortcoming is that the Wasserstein GAN model slightly overestimates the value of the 99 th percentile for the two Gaussian distributions. Comparing with the results shown in Table 1 on page 24 for the Bernoulli RBM and Table 2 on page 26 for the Gaussian RBM, we notice that the Wasserstein GAN model learns better the marginal distributions of the training dataset and doesn't have tendency to always underestimate or overestimate the tails of probability distribution as RBMs. However, there exists a small bias between the empirical mean computed with 50 Monte Carlo replications and the mean of the training dataset.
In Figure 21, we also compare the empirical correlation matrix of the training dataset and the average of the correlation matrices computed with 50 Monte Carlo replications. We notice that the Wasserstein GAN model captures very well the correlation structure.  2. We don't need to apply complex transformation to input data for Wasserstein GAN model in data preprocessing. In our study, we just use a MinMax scaling function and we recall that we need to use the binary transformation for Bernoulli RBM and the normal score transformation for Gaussian RBM.
3. The Wasserstein GAN model may capture more accurately the correlation structure of the training dataset than Bernoulli and Gaussian RBMs.

Application to financial time series
According to the results in the previous section, the Wasserstein GAN model preforms very well in the case of multi-dimensional simulated dataset without the time dependence relationship. In this paragraph, we construct a more complex Wasserstein GAN model with convolutional layers in the generator and the discriminator in order to extract more features in the real financial dataset consisting of 5 354 historical daily returns of the S&P 500 index and the VIX index from December 1998 to May 2018. To capture the time dependence of the multi-dimensional autocorrelated time series, we also need to associate the training samples with conditional labels as we have done in the case of the conditional RBM. Mariani et al. (2019) proposed a method to train the model with historical returns over n t consecutive days conditioned by the last n h day historical returns. As a result, the generator of this model should generate a matrix belonging to M nt,nx (R), which means the generator simulate directly daily returns of S&P 500 index and the VIX index for n t days. For reasons of simplicity, this model is named in this paper as the conditional deep convolutional Wasserstein GAN model (or CDCWGAN). In the step of data preprocessing, we need to apply the MinMax scaling function on input data before the training process to ensure that historical returns are scaled into the range [0, 1].
Structure and training of CDCWGAN models In this study on financial time series, the training samples for the discriminator correspond to a n x -dimensional time series representing historical returns over n t days. Therefore, inputs are represented by a matrix belonging to M nt,nx (R). In addition, each training sample is conditioned by the past values of returns over n h days, which are represented by a matrix belonging to M n h ,nx (R). In practice, we concatenate these two matrices to form a matrix belonging to the M n h +nt,nx (R) and this matrix will be fed to the discriminator as inputs. In this study, the inputs of the discriminator will be down sampled using 4 convolutional layers with number of filters {16, 32, 64, 128}. We set the kernel length to 3 and the stride to 2 and we choose to use the leaky RELU function with α = 0.5 for each convolutional layer. For the last layer of the discriminator, we choose always a 1-neuron dense layer and use the tangent hyperbolic activation function for this neuron as in the case of the traditional Wasserstein GAN model.
For the generator, we modify the original method of Mariani et al. (2019) and we borrow the idea of conditional layer in the case of conditional RBM. The generator of our CDCW-GAN model will take two inputs: a 100-dimensional random noise vector and n h past values that are concatenated into a n h × n x -dimensional vector. We then feed these two vectors to the first layer of the generator, which is a dense layer and we will reshape the output to a two-dimensional n t × 256 matrix before passing them to the second layer. We construct the rest of the generator using 3 convolutional layers with number of filters {256, 64, 2}. The kernel length is set to 3 and the stride is set to 2. We also use a leaky RELU function with α = 0.5 for each convolutional layer and since the last convolutional layer will give directly the simulated scenarios, we choose the sigmoid activation function in order to get the value in range [0, 1].
For this conditional deep convolutional Wasserstein GAN model, we use always the RMSProp optimizer with a small learning rate. Indeed, the learning rate is set to η = 1.10 −4 and the batch size is set to 500. Table 6 summarizes the setting of the model implementation using Python and the TensorFlow library. Learning joint distribution and time dependence In this study, we choose to set the length of training data window n t to 5, which means that the generator of the model will generate at each time step the scenarios for the 5 next days. In addition, we also use the values of the last 20 days as a long memory for input data of the model in order to compare with the results of the conditional RBM. After the training process, we generate samples of consecutive time steps to verify the quality of learning the joint distribution of daily returns of the S&P index and the VIX index. We run 50 Monte Carlo simulations and for each simulation, we generate 5 354 simulated observations from different random noise series. According to Table 7, we find that the CDCWGAN model learns also pretty well the joint distribution of the S&P index and the VIX index. In addition, the average correlation coefficient between daily returns of the two indices over 50 Monte Carlo replications is equal to −71%, which is exactly the figure of the empirical correlation. In Figure 22, we observe clearly that a CDCWGAN model can capture well the autocorrelation of the training dataset as in the case of the conditional RBM shown in Figure 17 on page 29.

VIX index
Market generator Based on above results, we consider that the CDCWGAN model can be also used as a market generator. In the example of S&P 500 and VIX indices, we have trained the model with all historical observations using the method proposed by Mariani et al. (2019) and the generator of this method will generate scenarios over several consecutive days. After having calibrated the training process, we may choose a date as the starting point and generate new samples for several days. We then update the memory vector and generate samples for the following days. Iteratively, the CDCWGAN model can generate a multi-dimensional time series that has the same patterns as those of the training data in terms of marginal distribution, correlation structure and time dependence. Comparing with the conditional RBM that we have studied, we don't need the normal score transformation to ensures the learning of marginal distributions and, convolutional layers in the generator and the discriminator may help us to extract more features in the training dataset.
As we have shown in Figure 18 on page 30 for the conditional RBM, Figure 23 shows three time series of 250 trading days starting at the date 22/11/2000, which are generated by the trained CDCWGAN model. We notice that these three simulated time series have clearly the positive autocorrelation in each dimension and they capture very well the negative correlation between S&P 500 and VIX indices. Comparing with the historical prices that are represented by the dashed line, these three simulated time series have the same characteristics.

Managing the out-of-sample robustness of backtesting methods
We now consider an application of generative models in the context of quantitative asset management. Traditionally, we use the historical daily returns of assets to backtest an investment strategy, meaning that only one real time series is available. Therefore, we get only one estimated value for the performance and risk statistics of the strategy, such as the annualized return, the volatility, the Sharpe ratio, the maximum drawdown, etc. Generative models can be used here to manage the out-of-sample robustness of backtesting methods. For instance, we may use a generative model like RBMs and GANs to learn the joint distribution and the time dependence of historical daily returns of assets. We may then use these models to generate new time series. Finally, we may backtest our investment strategy using these new simulated time series in order to construct the probability distribution of the different statistics. In this approach, the estimated value of the statistic obtained with the

VIX index
true real time series becomes one realization of its probability distribution. Suppose for example that the maximum drawdown of the backtest is equal to 10%, and that the live investment strategy has a maximum drawdown of 20%. Does it mean that the process of the investment strategy has been overfitted? Certainly yes if there is a zero probability to experience a maximum drawdown of 20% when the strategy is backtested with generative models. Definitively not if some samples of generative models have produced a maximum drawdown larger than 20%.

Backtest of the risk parity strategy
Our dataset consists of daily returns of six futures contracts on world-wide equity indices such as S&P 500, Eurostoxx 50 and Nikkei 225 indices and 10Y sovereign bonds such as US, Germany and Australia from January 2007 to December 2019. In this study, we build a risk parity strategy on this multi-asset investment universe. Moreover, the strategy is unfunded and we calibrate its leverage in order to obtain a volatility around 3%. In Figure 24, we have reported the cumulative performance of the risk parity strategy. Moreover, the descriptive statistics of performance and risk are given in Table 8, and correspond to the annualized performance µ (x), the volatility σ (x), the Sharpe ratio SR (x), the maximum drawdown 22 MDD (x) and the skew measure ξ (x), which is the ratio between the maximum drawdown and the volatility 23 . The maximum drawdown corresponds to a loss. For example, if the maximum drawdown is equal to 10%, the investor may face a maximum loss of 10%. This is why the maximum drawdown is expressed by a positive value. 23 If ξ (x) is greater than 3, this indicates that the strategy has a high skewness risk.

Market generator of the risk parity strategy
We first train the conditional RBM introduced in Section 3.1.2 on page 27 on the daily returns of the six futures contracts. In this study, the training of the models is done with all historical observations of futures contracts returns from January 2007 to December 2019. Moreover, for each date, we have used the values of the last 20 days as a memory vector for the conditional layer. We also choose a large number for the number of hidden layers units in order to capture the complex correlation structure of financial time series. Therefore, the conditional RBM has 6 visible layer units, 256 hidden layer units and 120 conditional layer units. In addition, we use the normal score transformation to avoid the problem of outliers in the training dataset and improve the learning quality. After having calibrated the training process, we choose the 1 st January 2018 as the starting point and generate new futures contracts returns iteratively for the next two years. We run 500 Monte Carlo simulations with different random noise series and for each simulation, we backtest our risk parity strategy using synthetic time series and calculate the descriptive statistics of the strategy. Finally, we construct probability distributions for these descriptive statistics and we compare them with the results obtained by the traditional method of bootstrap sampling, which consists of using a random sampling of the historical returns with replacement. Conditional RBM Figure 25 shows the cumulative performance of ten risk parity backtests using synthetic time series generated by the two methods. Among these simulations, we notice that the cumulative performance of the risk parity strategy using synthetic time series generated by the bootstrap sampling method is more centered around the real time series of backtesting. This phenomenon corresponds to the drawback of traditional bootstrap sampling method, since it cannot capture the time dependence in risk parity strategy. As shown in Figure 26, our risk parity strategy has a first-lag positive autocorrelation of 20%, and we compare the average autocorrelation function over 500 Monte Carlo simulations generated by the two approaches. We notice that the bootstrap sampling method cannot replicate this positive autocorrelation as the conditional RBM can do. This property is very important when we backtest meta-strategies, which means a strategy of an existing strategy. For instance, if we want to design a stop-loss strategy for our risk parity strategy, we should tune several parameters for implementing this stop-loss strategy. If we use the real historical time series to calibrate the stop-loss parameters, we might fall into the trap of overfitting because we only have one real historical time series. Therefore, if we use the time series generated by the bootstrap sampling method, we will not find the appropriate values of the stop-loss parameters, since the autocorrelation plays a key role in the stop-loss strategy as explained by Kaminski and Lo (2014). Using the conditional RBM as a market generator to generate time series is then a better way to find the appropriate parameters for the meta-strategy and manage the out-of-sample robustness.

Building the probability distribution of backtest statistics
In our study, we are interested in building the probability distribution of the maximum drawdown MDD (x) and the skew measure ξ (x) for the risk parity strategy. Figure 27 shows the distributions of the maximum drawdown generated by the bootstrap sampling and conditional RBM methods. We notice that the distribution of the maximum drawdown generated by the conditional RBM are more centered and have less severe scenarios. For instance, the value of the maximum drawdown of the real risk parity strategy from January 2018 to December 2019 is 3.69% and this value corresponds respectively to the 61.5% quantile in the probability distribution generated by the bootstrap sampling method and the 69.1% quantile in the probability distribution generated by the conditional RBM. We consider that the difference comes from the quality of learning about the time dependence of the dataset.
Although we use a volatility targeting method to control the volatility of risk parity strategies, we cannot ensure that the strategy volatility is exactly equal to 3%. To avoid the influence of the difference between realized volatility, we plot in Figure 28 the probability distribution of the skew measure. In this case, we notice that the tail of the probability distribution generated by the conditional RBM is less fat than that generated by the bootstrap sampling method but has several extreme severe scenarios. The skew measure of the real risk parity strategy from January 2018 to December 2019 is equal to 1.05, and this value corresponds respectively to the 66.3% quantile in the probability distribution generated by the bootstrap sampling method and the 70.1% quantile in the probability distribution generated by the conditional RBM.

Comparison with Wasserstein GAN models
For the purpose of illustration, we train also a conditional Wasserstein GAN model with a simple structure as multi-layer perceptrons for the generator and the discriminator. In order to make a comparison with the results of the conditional RBM, we also use the values of the last 20 days as a long memory for input data of the model. More precisely, the generator will take two inputs: a 100-dimensional random noise vector and 20 past values that are concatenated into a 20 × 6-dimensional vector. We then construct the generator using 4 dense layers with the structure S G = {100, 50, 10, 6}. We use a leaky RELU function with α = 0.5 for each dense layer and a sigmoid activation function for the output layer. The discriminator also takes two inputs: a 6-dimensional vector of current daily returns of futures contracts and a 20 × 6-dimensional vector of past values. We then construct the discriminator using 4 dense layers with the structure S D = {100, 50, 10, 1}. For the first three dense layers, the activation function corresponds to a leaky RELU function with α = 0.5, whereas the activation function of the output layer is a tanh activation function.
To avoid the problem of outliers, we also use the normal score transformation as in the case of conditional RBMs before applying the MinMax scaling function on input data. In Figure 29, we compare the distributions of the skew measure of the risk parity strategy using synthetic time series generated by the conditional RBM and conditional Wasserstein GAN models. We notice that these two distributions are similar, except that we have several extreme severe scenarios in the case of the conditional RBM. We recall that the value of the skew measure of the real risk parity strategy from January 2018 to December 2019 is equal to 1.05 and this value corresponds respectively to the 70.1% quantile in the probability distribution generated by the conditional RBM and the 72.9% quantile in the probability distribution generated by the conditional Wasserstein GAN.

Augmenting the investment universe with market regime indicators
In the real world of finance, the correlation structure is very complex and learning the joint distribution of daily returns of futures contracts using only time series themselves is not sufficient. For instance, if we train the generative models by including the time series of the VIX index, it's sure that we will obtain different results. Indeed, since the daily returns of the VIX index has a very negative correlation with equity futures contracts and a positive correlation with bond futures contracts and its distribution is very leptokurtic, the generative models trained using historical returns of futures contracts and VIX index will generate more distinct and more severe scenarios. Figure 30 shows the histograms of the skew measure of the risk parity strategy using synthetic time series generated by two conditional RBMs. The first one is trained with only historical daily returns of futures contracts and the other is trained with historical daily returns of futures contracts and VIX index. In this figure, we notice that the probability distribution generated by the model trained with the VIX index has a fatter tail and more extreme scenarios. Therefore, we consider that this model may generate more realistic data than the model trained using only the time series of the future contracts. In order to have a high-quality market generator, we believe that we should train a conditional RBM or Wasserstein GAN model with not only the time series of assets that compose the investment portfolio, but also those of the several market regime indicators such as the VIX index.

Conclusion
In this article, we explore the use of generative models for improving the robustness of trading strategies. We consider two approaches for simulating financial time series. The first one is based on restricted Boltzmann machines, whereas the second approach corresponds to the family of generative adversarial networks, including Wasserstein distance models. Given an historical sample of financial data, we show how to generate new samples of financial data using these techniques. These new samples of financial data are called synthetic or fake financial time series, and the objective is to preserve the statistical properties of the original sample. By statistical properties, we mean the statistical moments of each univariate time series, the stochastic dependence between the different variables that compose the multidimensional random vector, and also the time dependence between the observations. If we consider the financial times series as a multi-dimensional data matrix, the challenge is then to model both the row and column stochastic structures.
There are few satisfactory methods for simulating non-gaussian multi-dimensional financial time series. For instance, the bootstrap method does not preserve the cross-correlation between the different variables. The copula method is better, but it must use the techniques of conditional augmented data in order to reproduce the autocorrelation functions. Restricted Boltzmann machines and generative adversarial networks have been successful for generating complex data with non-linear dependence structures. By applying them to financial market data, our first results are encouraging and show that these new alternative techniques may help for simulating non-gaussian multi-dimensional financial time series. This is particularly true when we consider the backtesting of trading strategies. In this case, RBMs and GANs may be used for estimating the probability distribution of performance and risk statistics of the backtest. This opens the door to a new field of research for improving the risk management of quantitative investment strategies. over G. This means that there exists a set of strictly positive functions {ψ C , C ∈ C}, such that the joint probability distribution is given by a product of factors: where ψ C is the potential function for the clique C, v C are all the nodes belonging to the clique C and Z is the partition function given by: The partition function Z is the normalization constant what ensures the overall distribution sums to 1.

A.1.4 Energy function and Boltzmann distribution
The Hammersley-Clifford theorem is valid only if each potential function ψ C is strictly positive. Thus, we can introduce a new function E in order to rewrite the probability distribution P (x): where E (x) = − C∈C log ψ C (x (v C )) is called the energy function. Because natural exponential function is always positive, this guarantees that the energy function will result in a positive probability for any state. In addition, a large value of energy indicates a low probability of the state. According to LeCun et al. (2007), models of this form are called energy-based models. Using the energy function described above, the strictly positive probability distribution of a Markov random field can be expressed in the form P (x) = Z −1 e −E(x) . This form of distribution is also called Boltzmann (or Gibbs) distribution for a system in statistical physics. It is defined as: kT where p i is the probability of state i of the system, E i is the energy of state i, k is the Boltzmann constant, and T is the temperature of the system. Since N is the number of states accessible to the system, the normalization constant is Z T = N i=1 e −Ei/(kT ) . If we set kT to 1, we find that the Boltzmann distribution and the probability distribution of a Markov random field have the same formula. For this reason, many energy-based models are called Boltzmann machines.

A.1.5 The example of restricted Boltzmann machines
The restricted Boltzmann machine introduced in Section 2.1 on page 3 is a Markov random field associated with a bipartite undirected graph G = (V, E). In other words, all visible layer units and hidden layer units can be considered as nodes in an undirected graph 25 : . . , V m , H 1 , H 2 , . . . , H n } 25 For simplicity reasons, we make a little abuse of notation here to use V i and H j to represent not only nodes in the graph but also random variables associated with these nodes and we use v i and h j to denote respectively the possible values of the variables associated with the i th visible unit and j th hidden unit. and all connections between visible layer and hidden layer are edges of G. In the case of RBMs, we know that there are only cliques of size 1 (one visible unit or one hidden unit) and cliques of size 2 (a pair of one visible unit and one hidden unit) in the graph G. In addition, it is easy to show that all these cliques are maximal. Let C 1 and C 2 be respectively the set of all the cliques of size 1 and the set of all the cliques of size 2. We obtain: According to the Hammersley-Clifford theorem, the probability distribution of an RBM is given by: (v,h) where: In the case of Bernoulli RBMs introduced in Section 2.1.1 on page 4, we defined It follows that the energy function of a Bernoulli RBM is equal to: where a i and b j are bias terms associated with the visible and hidden variables V i and H j , and w i,j is the weight associated with the edge between V i and H j .

A.2.2 Bernoulli RBMs and neural networks
Following Fischer and Igel (2014), we divide the energy function E (v, h) into two parts: one collecting all terms involving v i and one collecting all the other terms v −i : where: The Bayes theorem gives: We deduce that the conditional probability P (v i = 1 | h) is equal to: where σ (x) is the sigmoid function: Similarly, we can divide the energy function E (v, h) into two parts: one collecting all terms involving h j and one collecting all the other terms h −j : where: Using the same approach as previously, we can show that:

A.2.3 Gradient of the Bernoulli RBM log-likelihood function
We have: Fischer and Igel (2014) computed the log-likelihood gradient ∇ θ (v) = ∂ θ (θ | v): We have 26 : and: We recall that ∂ ai E (v, h) = −v i and ∂ bj E (v, h) = −h j . It follows that: For the gradient with respect to W , we have ∂ wi,j E (v, h) = −v i h j and: Using the Bayes theorem, the conditional probability distribution P (h | v) is equal to: We notice that we have to sum over 2 m possible combinations of the visible variables when calculating the second term ∇ (2) θ (v). Therefore, we generally approximate the expectation P (v ) by sampling from the model distribution.

A.2.4 Gradient of the contrastive divergence
The contrastive divergence function is equal to: where: We have: When P does not depend on the parameter set θ, we have ∂ θ P (v) = 0 and the derivative reduces to 28 : we also notice that another expression of Equation (25) is: Finally, we deduce that: In Equation (26), it is often possible to compute the exact values for the first two terms, but not for the third term. However, Hinton (2002) showed experimentally that this last term is so small compared the other two terms that it can be ignored. Thus, we can consider the following approximation: (s) is the s th sample of the training set and v (k) (s) is the associated sample after running k steps of Gibbs sampling 29 . By noticing that log P θ (v) = (θ | v) and using Equation (24), we obtain: Therefore, we can compute the derivatives with respect to the parameters a i , b j and w i,j :

A.2.5 Gradient of the Gaussian-Bernoulli RBM log-likelihood function
By updating Equation (24), we can easily show that the gradient of the log-likelihood function is equal to: where: (s) is the starting value of the Gibbs sampler.
We have: We deduce that the derivative of (θ | v) with respect to the weight w i,j is given by: Similarly, we can compute the other derivatives. For a i and b j , we obtain: Finally, we obtain for the parameter σ i :

A.2.6 Gradient of the conditional RBM log-likelihood function
The log-likelihood function for the conditional RBM is: where the set of parameters becomes θ = (a, b, W, P, Q). We can then compute the partial derivative of the energy function by using the chain rule: and we obtain:

Improving the Robustness of Backtesting with Generative Models
Similarly, we have: and: The calculation for the other derivatives remains unchanged and are the same as those obtained for the Gaussian-Bernoulli RBM.

A.3 A unified approach of GAN models
Generative models are trained to perform a mapping from a latent space to some specified data manifold, which is generally represented by the empirical distribution of real data. The problem consists then in finding a mapping function that best matches the target data in the sense of a certain discrepancy measure. For comparing the theoretical distribution with the empirical distribution, the two important metrics used in the context of generative modeling are divergence measures (φ-divergence) and integral probability metrics (IPM), which give two different families of GANs. In this section, we show that all models share many significant common points. Let X be a topological space that is compact, complete and separable. Let us assume that there are two probability measures P and Q which can be defined on X . GAN optimization relies on the fact that both φ-divergence and IPM can be written as: where F is the class of functions ϕ defined on X and ∆ : F → X is a discrepancy operator.
In the following, we show that the choice of the class F and the discrepancy operator ∆ lead to different GAN models.

A.3.1 φ-GAN models
Variational estimation of φ-divergences The very first model of GAN is often presented as a minmax optimization problem (Goodfellow et al., 2014). However, it is possible to find a direct correspondence between Goodfellow's saddle point problem and divergence minimization. Let us recall the definition of a divergence measure. We assume that there are two probability measures P and Q that can be defined on X . Moreover, P must be absolutely continuous 30 with respect to Q, which is denoted P Q. The φ-divergence D φ is defined as: where φ : R + → R∪{+∞} is a convex, lower-semi-continuous function such that 31 φ (1) = 0. Looking at the closed-form solution of φ-divergence, we note that it can be interpreted as the likelihood ratio between two probability distributions 32 .
In order to establish the link with GAN optimization problems, Nguyen et al. (2010) proposed to compute a variational characterization of these φ-divergence measures by looking at the convex dual. For that, we need to introduce the Fenchel conjugate, which is a 30 This means that if we consider a σ-field A ⊆ X such that Q (A) = 0, then P (A) = 0. 31 This last condition ensures that D φ (P Q) = 0 if P = Q. 32 The condition that imposes that P must be absolutely continuous with respect to Q is related to the Radon-Nikodym theorem. It states that if P Q, then there is a function f : X → R that satisfies P (A) = A f (x) Q (dx) for all σ-field A. The function f is often denoted by dP dQ .
However, evaluating this quantity is impossible because the distribution function P is unknown. Therefore, ϕ should be flexible enough to approximate the derivative φ everywhere. This is why GAN models use deep neural networks to estimate it. This leads us to introduce the parameter θ d that will be optimized during the training process. In this context, we write the parameterized function D (x, θ d ), which aims to estimate the function ϕ (x): Consequently, Nowozin et al. (2016) proposed to use the resulting lower bound in order to train GANs. G (z, θ g ) represents the generative model that is also a neural network and allows us to build the probability distribution P model , which should estimate the given probability distribution P data . Thus, we obtain the saddle point problem: where: D (x, θ d ) ∈ F and P model = G (z, θ g ). Therefore, we have constrained the generator to be in a smaller class of functions belonging to F even if the neural network can approximate any function.
Remark 8. In order to obtain the minimax problem of Goodfellow et al. (2014), Nowozin et al. (2016) considered the function: φ (x) = x log x − (x + 1) log (x + 1) We deduce that: φ (x) = log x x + 1 and 34 : φ * (x) = − log (1 − e x ) 34 We have: φ * (t) = sup x∈dom φ x, t − φ (x) It follows that h (x) = xt − x log x + (x + 1) log (x + 1) and: The supremum is reached at the point h (x ) = 0, implying that: x = e t 1 − e t Finally, we obtain: Using Equation (30), we deduce that: If we set ϕ (X) = log D (X, θ d ), we find the minimax function of Goodfellow et al. (2014): Another representation of φ-divergences GAN training can be viewed as a process of successively estimating the optimal function ϕ and minimizing the φ-divergence. Chu et al. (2019) proposed a more general view of this process. Previously, the study has been done on the space R thanks to the dual form of the given function φ associated with the φ-divergence. Alternatively, Chu et al. (2019) proposed to directly work on the probability space X . This imposes to define a probability functional J P (Q) : B (X ) −→ R, where B (X ) is the space of Borel probability measures on X . In the context of GAN optimization problems, we would like to show that: min D φ (P Q) = min However, contrary to the previous case, probability functionals take value in probability spaces, where functional derivatives need to be defined. Moreover, the dimension of the space is potentially infinite. Therefore, the difficulty is to transform functional optimization into a convex optimization problem that can be solved using traditional numerical algorithms such as the gradient descent. In order to define the derivative in B (X ), Chu et al. (2019) used the Gâteaux derivative of the functional J P at Q in the direction of δ = P − Q: This definition initially comes from Von Mises calculus, and the Gâteaux derivative is also called the 'Volterra' derivative (Fernholz, 2012). Chu et al. (2019) recalled that the Gâteaux derivative has an integral representation J P (δ) = X g (x) δ (dx) where the function g : X → R is called the 'influence function' or the influence curve 35 . It follows that: We are now able to compute the derivative in order to recover the optimal function that satisfy the minimum of the φ-divergence: We deduce that the influence function for the φ-divergence is equal to: The influence function can then be associated with the optimal function φ introduced before because we have g φ = ϕ ∈ F. We conclude that the GAN discriminator will try to estimate the influence function. While Nguyen et al. (2010) focused on the dual form of the function φ, Chu et al. (2019) used a more general approach by considering the dual form of the probability functional in order to recover the Goodfellow's saddle point problem.
In the case of probability functionals which take values in B (X ), the dual space can be defined as the space L 1 (X ) of all real-valued Lipschitz function that takes values in X (Laschos et al., 2019). According to the Riesz representation theorem, it is possible to identify the space B (X ) to it dual. We consider the product space B (X ) × L 1 (X ) with the scalar product defined as: where ϕ ∈ L 1 (X ) and Q ∈ B (X ) × R. Thus, we can rewrite J P in term of its convex transform: To bridge the gap between the two approaches, we have to demonstrate that J * P (ϕ) = E [φ * (ϕ (X)) | X ∼ P] for a given estimate of the function ϕ ∈ F. For that, Chu et al. (2019) considered the Fenchel conjugate of the Jensen-Shannon divergence 36 :

A.3.2 IPM-GAN models
In the case of Wasserstein generative adversarial networks introduced by Arjovsky et al. (2017a,b), φ-divergences have to be replaced by integral probability metrics (IPMs) in order to compare two different probability distributions. Let F be a class of functions defined on X . Müller (1997) defined an IPM I F between P and Q in the following way: IPMs are looking for a critic function that maximizes the average discrepancy between the two distributions P and Q. Contrary to φ-GAN models, where we have to find the variational form of the divergence φ (x), the definition of an IPM directly gives the discrepancy operator: for all ϕ ∈ F. The Wasserstein GAN proposed by Arjovsky et al. (2017a,b) considers the function class F such that ϕ (x) is a 1-Lipschitz function. In order to present different choices for the function class F, we consider the Lebesgue norm on the measurable space Ω = (X , P): f 2 2 = X f 2 (x) P (dx). Let us denote the normed space by L 2 (Ω) = {f : X → R | f 2 < +∞} and the unit ball by B 1 = f ∈ L 2 (Ω) | f 2 ≤ 1 . Therefore, choosing the class function such that F = B 1 is called a Fisher GAN by . In a similar way, Mroueh et al. (2018) proposed to define Sobolev GAN models by considering the following class of functions F = f ∈ L 2 (Ω) | ∇ x f 2 ≤ 1 .

Improving the Robustness of Backtesting with Generative Models
The convex conjugate of J JS is then: Remark 9. In order to retrieve Equation (40), we have to interchange P and Q, because we need to compute J P (Q) and not J Q (P).

A.5 Derivation of the minimax cost function
The cost function can be viewed as a binary cross-entropy measure. Let Y andŶ be two random variables with probability mass function p and q. The cross-entropy function is equal to: For discrete probability distributions, we obtain: In a binary classification problem, for a given observation x i , we have p = y i ∈ {0, 1}, which is the true label and q =ŷ i ∈ [0, 1] which is the predicted probability of the current model. We can use binary cross entropy to get a measure of dissimilarity between y i andŷ i : In the case of m samples, the loss function is then given by: Under the GAN framework, we have m samples of x 0 and m samples of x 1 , which serve as the input data of the discriminator model. We note x = {x 0 , x 1 } the set of the two samples. Sinceŷ i = D (x i ; θ d ), the formula above can be written as: If x i comes from the sample x 0 , y i takes the value 0, otherwise it takes the value 1. It follows that: Therefore, the loss function L (θ g , θ d ) corresponds to the average of the cross-entropy when considering several observations: We notice that C (θ g , θ d ) is equal to −2 · L (θ g , θ d ). Minimizing the loss function is then equivalent to maximize C (θ g , θ d ) with respect to θ d .