Justin Sirignano

Rama Cont

Universal features of price formation in financial markets: perspectives from Deep Learning

Using a large-scale Deep Learning approach applied to a high-frequency database containing billions of electronic market quotes and transactions for US equities, we uncover nonparametric evidence for the existence of a universal and stationary price formation mechanism relating the dynamics of supply and demand for a stock, as revealed through the order book, to subsequent variations in its market price. We assess the model by testing its out-of-sample predictions for the direction of price moves given the history of price and order flow, across a wide range of stocks and time periods. The universal price formation model exhibits a remarkably stable out-of-sample prediction accuracy across time, for a wide range of stocks from different sectors. Interestingly, these results also hold for stocks which are not part of the training sample, showing that the relations captured by the model are universal and not asset-specific.

The universal model -trained on data from all stocks -outperforms, in terms of out-of-sample prediction accuracy, asset-specific linear and nonlinear models trained on time series of any given stock, showing that the universal nature of price formation weighs in favour of pooling together financial data from various stocks, rather than designing asset-or sector-specific models as commonly done. Standard data normalizations based on volatility, price level or average spread, or partitioning the training data into sectors or categories such as large/small tick stocks, do not improve training results. On the other hand, inclusion of price and order flow history over many past observations improves forecasting performance, showing evidence of path-dependence in price dynamics.

1 Price formation: how market prices react to supply and demand

The computerization of financial markets and the availability of detailed electronic records of order flow and price dynamics in financial markets over the last decade has unleashed TeraBytes of high frequency data on transactions, order flow and order book dynamics in listed markets, which provide us with a detailed view of the high-frequency dynamics of supply, demand and price in these markets [START_REF] Cont | Statistical modeling of high frequency financial data: Facts, models and challenges[END_REF]. This data may be put to use to explore the nature of the price formation mechanism which describes how market prices react to fluctuations in supply and demand. At a high level, a 'price formation mechanism' is a map which represents the relationship between the market price and variables such as price history and order flow:

Price(t + ∆t) = F Price history(0...t), Order Flow(0...t), Other Information = F (X t , t ), where X t is a set of state variables (e.g., lagged values of price, volatility, and order flow), endowed with some dynamics, and t is a random 'noise' or innovation term representing the arrival of new information and other effects not captured entirely by the state variables. Empirical and theoretical market microstructure models, stochastic models and machine learning price prediction models can all be viewed as different ways of representing this map F , at various time resolutions ∆t. One question, which has been implicit in the literature, is the degree to which this map F is universal (i.e., independent of the specific asset being considered). The generic, as opposed to asset-specific, formulation of market microstructure models seems to implicitly assume such a universality. Empirical evidence on the universality of certain stylized facts [START_REF] Cont | Empirical properties of asset returns: stylized facts and statistical issues[END_REF] and scaling relations [START_REF] Benzaquen | Unravelling the trading invariance hypothesis[END_REF][START_REF] Andersen | Intraday Trading Invariance in the E-Mini S&P 500 Futures Market[END_REF][START_REF] Kyle | Market microstructure invariance: Empirical hypotheses[END_REF][START_REF] Mandelbrot | The multifractal model of asset returns[END_REF] seems to support the universality hypothesis. Yet, the practice of statistical modeling of financial time series has remained asset specific: when building a model for the returns of a given asset, market practitioners and econometricians only use data from the same asset. For example, a model for Microsoft would only be estimated using Microsoft data, and would not use data from other stocks. Furthermore, the data used for estimation is often limited to a recent time window, reflecting the belief that financial data can be 'non-stationary' and prone to regime changes which may render older data less relevant for prediction.

Due to such considerations, models considered in financial econometrics, trading and risk management applications are asset-specific and their parameters are (re)estimated over time using a time window of recent data. That is, for asset i at time t the model assumes the form Price i (t + ∆t) = F X i 0:t , t | θ i (t) , where the model parameter θ i (t) is periodically updated using recent data on price and other state variables related to asset i. As a result, data sets are fragmented across assets and time and, even in the high frequency realm, the size of data sets used for model estimation and training are orders of magnitude smaller than those encountered in other fields where Big Data analytics have been successfully applied. This is one of the reasons why, except in a few instances [START_REF] Buhler | Deep hedging[END_REF][START_REF] Dixon | Sequence Classification of the Limit Order Book using Recurrent Neural Networks[END_REF][START_REF] Kolanovic | Big Data and AI Strategies: Machine Learning and Alternative Data Approach to Investing[END_REF][START_REF] Sirignano | Stochastic Gradient Descent in Continuous Time[END_REF][START_REF] Sirignano | Deep Learning for Limit Order Books[END_REF][START_REF] Sirignano | Deep Learning for Mortgage Risk[END_REF], large-scale learning methods such as Deep Learning [START_REF] Goodfellow | Deep Learning[END_REF] have not been deployed for quantitative modeling in finance. In particular, the non-stationarity argument is sometimes invoked to warn against their use.

On the other hand, if the relation between these variables were universal and stationary, i.e. if the parameter θ i (t) varies neither with the asset i nor with time t, then one could potentially pool data across different assets and time periods and use a much richer data set to estimate/ train the model. For instance, data on a flash crash episode in one asset market could provide insights into how the price of another asset would react to severe imbalances in order flow, whether or not such an episode has occurred in its history.

In this work, we provide evidence for the existence of such a universal, stationary relation between order flow and market price fluctuations, using a nonparametric approach based on Deep Learning. Deep learning can estimate nonlinear relations between variables using 'deep' multilayer neural networks which are trained on large data sets using 'supervised learning' methods [START_REF] Goodfellow | Deep Learning[END_REF].

Using a deep neural network architecture trained on a high-frequency database containing billions of electronic market transactions and quotes for US equities, we uncover nonparametric evidence for the existence of a universal and stationary price formation mechanism relating the dynamics of supply and demand for a stock, as revealed through the order book, to subsequent variations in its market price. We assess the model by testing its out-of-sample predictions for the direction of price moves given the history of price and order flow, across a wide range of stocks and time periods. The universal price formation model exhibits a remarkably stable out-of-sample prediction accuracy across time, for a wide range of stocks from different sectors. Interestingly, these results also hold for stocks which are not part of the training sample, showing that the relations captured by the model are universal and not asset-specific. We observe that the neural network thus trained outperforms linear models, pointing to the presence of nonlinear relationships between order flow and price changes.

Our paper provides quantitative evidence for the existence of a universal price formation mechanism in financial markets. The universal nature of the price formation mechanism is reflected by the fact that a model trained on data from all stocks outperforms, in terms of out-of-sample prediction accuracy, stock-specific linear and nonlinear models trained on time series of any given stock. This shows that the universal nature of price formation weighs in favour of pooling together financial data from various stocks, rather than designing stock-or sector-specific models as commonly done. Also, we observe that standard data transformations such as normalizations based on volatility or average spread, or partitioning the training data into sectors or categories such as large/small tick stocks, do not improve training results. On the other hand, inclusion of price and order flow history over many past observations improves forecasting performance, showing evidence of path-dependence in price dynamics.

Remarkably, the universal model is able to extrapolate, or generalize, to stocks not within the training set. The universal model is able to perform well on completely new stocks whose historical data the model was never trained on. This implies that the universal model captures features of the price formation mechanism which are robust across stocks and sectors. This feature alone is quite interesting for applications in finance where missing data problems and newly issued securities often complicate model estimation.

Outline Section 2 describes the dataset and the supervised learning approach used to extract information about the price formation mechanism. Section 3 provides evidence for the existence of a universal and stationary relationship linking order flow and price history to price variations. Section 4 summarizes our main findings and discusses some implications.

A data-driven model of price formation via Deep Learning

Applications such as image, text, and speech recognition have been revolutionized by the advent of 'Deep Learning' -the use of multilayer ('deep') neural networks trained on large data sets to uncover complex nonlinear relations between high-dimensional inputs ('features') and outputs [START_REF] Goodfellow | Deep Learning[END_REF]. At an abstract level, a deep neural network represents a functional relation y = f (x) between a high-dimensional input vector x and an output y through iterations ('layers') consisting of weighted sums followed by the application of nonlinear 'activation' functions. Each iteration corresponds to a 'hidden layer' and a deep neural network can have many hidden layers. Neural networks can be used as 'universal approximators' for complex nonlinear relationships [START_REF] Hornik | Multilayer Feedforward Networks are Universal Approximators[END_REF], by appropriately choosing the weights in each layer.

In supervised learning approaches, network weights are estimated by optimizing a regularized cost function reflecting in-sample discrepancy between the network output and desired outputs. In a deep neural network, this represents a high-dimensional optimization over hundreds of thousands (or even millions) of parameters. This optimization is computationally intensive due to the large number of parameters and large amount of data. Stochastic gradient descent algorithms (e.g., RMSprop or ADAM) are used for training neural networks, and training is parallelized on Graphics Processing Units (GPUs).

We apply this approach to learn the relation between supply and demand on an electronic exchange -captured in the history of the order book for each stock -and the subsequent variation of the market price. Our data set is a high-frequency record of all orders, transactions and order cancellations for approximately 1000 stocks traded on the NASDAQ between Figure 1: The limit order book represents a snapshot of the supply and demand for a stock on an electronic exchange. The 'ask' side represents sell orders and the 'bid' side, buy orders. The size represents the number of shares available for sale/purchase at a given price. The difference between the lowest sell price (ask) and the highest buy price (bid) is the 'spread' (in this example, 1 ¢).

Electronic buy and sell orders are continuously submitted, cancelled and executed through the exchange's order book. A 'limit order' is a buy or sell order for a stock at a certain price and will appear in the order book at that price and remain there until cancelled or executed. The 'limit order book' is a snapshot of all outstanding limit orders and thus represents the visible supply and demand for the stock (see Figure 1). In US stock markets, orders can be submitted at prices occurring at multiples of 1 cent. The 'best ask price' is the lowest sell order and the 'best bid price' is the highest bid price. The best ask price and best bid price are the prices at which the stock can be immediately bought or sold. The 'mid-price' is the average of the best ask price and best bid price. The order book evolves over time as new orders are submitted, existing orders are cancelled, and trades are executed.

In electronic markets such as the NASDAQ, new orders may arrive at high frequencysometimes every microsecond -and order books of certain stocks can update millions of times per day. This leads to TeraBytes of data, which we put to use to build a data-driven model of the price formation process.

When the input data is a time series, causality constraints require that the relation between input and output respects the ordering in time. Only the past may affect the present. A network architecture which reflects this constraint is a recurrent network (see an example in Figure 2) based on Long Short-Term Memory (LSTM) units [START_REF] Gers | Learning to Forget: Continual Prediction with LSTM[END_REF].

Each LSTM unit has an internal state which maintains a nonlinear representation of all past data. This internal state is updated as new data arrives. Our network has 3 layers of LSTM units followed by a final feed-forward layer of rectified linear units (ReLUs). A probability distribution for the next price move is produced by applying a softmax activation function. LSTM units are specially designed to efficiently encode the temporal sequence of data [START_REF] Gers | Learning to Forget: Continual Prediction with LSTM[END_REF][START_REF] Goodfellow | Deep Learning[END_REF].

h 1 t-1 h n-m t-1 h n t-1 Y t-1 X t-1 h 1 t-1(1) h n-m t-1(1) h n t-1(1) h 1 t-1(2) h n-m t-1(2) h n t-1(2) h 1 t h n-m t h n t Y t X t h 1 t(1) h n-m t(1) h n t(1) h 1 t(2) h n-m t(2) h n t(2) h 1 t+1 h n-m t+1 h n t+1 Y t+1 X t+1 . . .
We train the network to forecast the next price move from a vector of state variables, which encode the history of the order book over many observation lags. The index t represents the number of price changes. At a high level, the LSTM network is of the form

(Y t , h t ) = f (X t , h t-1 ; θ). (2.1)
Y t is the prediction for the next price move, X t is the state of the order book at time t, h t is the internal state of the deep learning model, reprenting information extracted from the history of X up to t, and θ designates the model parameters, which correspond to the weights in the neural network. At each time point t the model uses the current value of state variables X t (i.e. the current order book) and the nonlinear representation of all previous data h t-1 , which summarizes relevant features of the history of order flow, to predict the next price move. In principle, this allows for arbitrary history-dependence: the history of the state variables (X s , s ≤ t) may affect the evolution of the system, in particular price dynamics, at all future times T ≥ t in a nonlinear way. Alternative modeling approaches typically do not allow the flexibility of blending nonlinearity and history-dependence in this manner.

A supervised learning approach is used to learn the value of the (high-dimensional) parameter θ by minimizing a regularized negative log-likelihood objective function using a stochastic gradient descent algorithm [START_REF] Goodfellow | Deep Learning[END_REF]. The parameter θ is assumed to be constant across time, so it affects the output at all times in a recursive manner. A stochastic gradient descent step at time t requires calculating the sensitivity of the output to θ, via a chain rule, back through the previous times t -1, t -2, . . . , t -T (commonly referred to as 'backpropagation through time'). In theory, backpropagation should occur back to time 0 (i.e., T = t). However, this is computationally impractical and we truncate the backpropagation at some lag T . In Section 3.4, we discuss the impact of the past history of the order book and the 'long memory' of the market.

The resulting LSTM network involves up to hundreds of thousands of parameters. This is relatively small compared to networks used for instance in image or speech recognition, but it is huge compared to econometric models traditionally used in finance. Previous literature has been almost entirely devoted to linear models or stochastic models with a very small number of parameters. It is commonly believed that financial data is far too noisy to build such large models without overfitting; our results show that this is not necessarily true.

Given the size of the data and the large number of network parameters to be learned, significant computational resources are required both for pre-processing the data and training the network. Training of deep neural networks can be highly parallelized on GPUs. Each GPU has thousands of cores, and training is typically 10× faster on a GPU than a standard CPU. The NASDAQ data was filtered to create training and test sets. This data processing is parallelized over approximately 500 compute nodes. Training of asset-specific models was also parallelized, with each stock assigned to a single GPU node. Approximately 500 GPU nodes are used to train the stock-specific models.

These asset-specific models, trained on the data related to a single stock, were then compared to a 'universal model' trained on the combined data from all the stocks in the dataset. Data from various stocks were pooled together for this purposes without any specific normalization.

Due to the large amount of data, we distributed the training of the universal model across 25 GPU nodes using asynchronous stochastic gradient descent (Figure 3). Each node loads a batch of data (selected at random from all stocks in the dataset), computes gradients of the model on the GPU, and then updates the model. Updates occur asynchronously, meaning node j updates the model without waiting for nodes i = j to finish their computations.

Results

We split the universe of stocks into two groups of roughly 500 stocks; training is done on transactions and quotes for stocks from the first group. We distinguish:

• stock-specific models, trained using data on all transactions and quotes for a specific stock.

• the 'universal model', trained using data on all transactions and quotes for all stocks in the training set.

All models are trained for predicting the direction of the next price move. Specifically, if τ 1 , τ 2 , . . . are the times at which the mid-price P t changes, we estimate P[P τ k+1 -P τ k > 0|X τ 0:k ] Figure 3: Asynchronous stochastic gradient descent for training the neural network. The dataset, which is too large to be held in the nodes' memories, is stored on the Online Storage system. Batches of data are randomly selected from all stocks and sent to the GPU nodes.

Gradients are calculated on the GPUs and then the model is asynchronously updated.

and P[P τ k+1 -P τ k < 0|X τ 0:k ] where X t is the state of the limit order book at time t. The models therefore predict whether the next price move is up or down. The events are irregularly spaced in time. The time interval τ k+1 -τ k between price moves can vary considerably from a fraction of a second to seconds. We measure the forecast accuracy of a model for a given stock via the proportion of observations for which it correctly predicts the direction of the next price move. This can be estimated using the empirical estimator

A i =
Number of price changes where model correctly predicts price direction for stock i Total number of price changes × 100%.

All results are out-of-sample in time. That is, the accuracy is evaluated on time periods outside of the training set. Model accuracy is reported via the cross-sectional distribution of the accuracy score A i across stocks in the testing sample, and models are compared by comparing their accuracy scores.

In addition, we evaluate the accuracy of the universal model for stocks outside the training set. Importantly, this means we assess forecast accuracy for stock i using a model which is trained without any data on stock i. This tests whether the universal model can generalize to completely new stocks.

Typically, the out-of-sample dataset is a 3-month time period. In the context of highfrequency data, 3 months corresponds to millions of observations and therefore provides a lot of scope for testing model performance and estimating model accuracy. In a data set with no stationary trend (as in the case at such high frequencies), a random forecast ('coin-flip') would yield an expected score of 50%. Given the large size of the data set, even a small deviation (i.e. 1%) from this 50% benchmark is statistically significant.

The main findings of our data-driven approach may be summarized as follows:

• Nonlinearity: Data-driven models trained using deep learning substantially outperform linear models in terms of forecasting accuracy (Section 3.1).

• Universality: The model uncovers universal features that are common across all stocks (Section 3.2). These features generalize well: they are also observed to hold for stocks which are not part of the training sample.

• Stationarity: model performance in terms of price forecasting accuracy is remarkably stable across time, even a year out of sample. This shows evidence for the existence of a stationary relationship between order flow and price changes (Section 3.3), which is stable over long time periods.

• Path-dependence and long-range dependence: inclusion of price and order flow history is shown to substantially increase the forecast accuracy. This provides evidence that price dynamics depend not only on the current or recent state of the limit order book but on its history, possibly over long time scales (Section 3.4).

Our results show that there is far more common structure across data from different financial instruments than previously thought. Providing a suitably flexible model is used which allows for nonlinearity and history-dependence, data from various assets may be pooled together to yield a data set large enough for deep learning.

Deep Learning versus Linear Models

Linear state space models, such as Vector Autoregressive (VAR) models, have been widely used in the modeling of high frequency data and in empirical market microstructure research [START_REF] Hasbrouck | Empirical Market Microstructure: The Institutions[END_REF] and provide a natural benchmark for evaluating the performance of a forecast. Linear models are easy to estimate and capture in a simple way the trends, linear correlations and autocorrelations in the state variables.

The results in Figure 4 show that the deep learning models substantially outperform linear models. Given the large sample size, an increase of 1% in accuracy is considered significant in the context of high-frequency modeling.

The linear (VAR) model may be formulated as follows: at each observation we update a vector of linear features h t and then use a probit model for the conditional probability of an upward price move given the state variables:

h t = Ah t-1 + BX t , Y t = P(∆P t > 0|X t , h t ) = G(CX t + Dh t ). (3.-1)
where G depends on the distributional assumptions on the innovations in the linear model. For example, if we use a logistic distribution for the innovations in the linear model, then the probability distribution of the next price move is given by softmax (logistic) function applied to a linear function of the current order book and linear features:

P(∆P t > 0|X t , h t ) = Softmax(CX t + Dh t ).
We compare the neural network against a linear model for approximately 500 stocks. To compare models we report the difference in accuracy scores across the same test data set. Let

• L i be the accuracy of the stock-specific linear model g θ i for asset i estimated on data only from stock i,

• Âi be the accuracy of the stock-specific deep learning model f θ i trained on data only from stock i, and

• A i be the accuracy for asset i of the universal deep learning model f θ trained on a pooled data set of all quotes and transactions for all stocks.

The left plot in Figure 4 reports the cross-sectional distribution for the increase in accuracy Âi -L i when moving from the stock-specific linear model to the stock-specific deep learning model. We observe a substantial increase in accuracy, between 5% to 10% for most stocks, when incorporating nonlinear effects using the neural networks. The right plot in Figure 4 displays histograms of A i (red) and L i (blue). We clearly observe that moving from a stock-specific linear model to the universal nonlinear model trained on all stocks substantially improves the forecasting accuracy by around 10%. The deep neural network outperforms the linear model since it is able to estimate nonlinear relationships between the price dynamics and the order book, which represents the visible supply and demand for the stock. This is consistent with an abundant empirical and econometric literature documenting nonlinear effects in financial time series, but the large amplitude of this improvement can be attributed to the flexibility of the neural network in representing nonlinearities.

More specifically, sensitivity analysis of our data-driven model uncovers stable nonlinear relations between state variables and price moves, i.e. nonlinear features which are useful for forecasting. Figure 5 presents an examples of such a feature: the relation between the depth on the bid and ask sides of the order book and the probability of a price decrease. Such relations have been studied in queueing models of limit order book dynamics [START_REF] Cont | A stochastic model for order book dynamics[END_REF][START_REF] Cont | Price dynamics in a Markovian limit order market[END_REF]. In particular, it was shown in [START_REF] Cont | Price dynamics in a Markovian limit order market[END_REF] that when the order flow is symmetric then there exists a 'universal' relation -not dependent on model parameters -between bid depth, ask depth and the probability of a price decrease at the next price move. However, the derivations in these models hinge on many statistical assumptions which may or may not hold, and the universality of such relations remained to be empirically verified.

Our analysis shows that there is indeed evidence for such a universal relation, across a wide range of assets and time periods. Figure 5 (left) displays the probability of a price decrease as a function of the depth (the number of shares) at the best bid/ask price. The larger the best ask size, the more likely the next price prove will be downwards. The probability is approximately constant along the center diagonal where the bid/ask imbalance is zero. However, as observed in queueing models [START_REF] Cont | A stochastic model for order book dynamics[END_REF][START_REF] Cont | Price dynamics in a Markovian limit order market[END_REF][START_REF] Figueroa-Lopez | One-level limit order book model with memory and variable spread[END_REF], even under simplifying assumptions, the relation between this probability and various measures of the bid/ask imbalance is not linear. Furthermore, such queueing models typically focus on the influence of depth at the top of the order book and it is more difficult to extract information from deeper levels of the order book. The right contour plot in Figure 5 displays the influence of limit orders deeper in the order book (here: total size aggregated across levels 5 to 10) on the probability of a price decrease. We see that the influence is less than the depth at the top of the book, as illustrated by the tighter range of predicted probabilities, but still significant.

Universality across assets

A striking aspect of our results is the stability across stocks of the features uncovered by the deep learning model, and its ability to extrapolate ('generalize') to stocks which it was not trained on. This may be illustrated by comparing forecasting accuracy of stock-specific models, trained only on data of a given stock, to a universal model trained on a pooled data set of 500 stocks, a much larger but extremely heterogeneous data set. As shown in Figure 6, which plots A i -Âi , the universal model consistently outperforms the stock-specific models. This indicates there are common features, relevant to forecasting, across all stocks. Features extracted from data on stock A may be relevant to forecasting of price moves for stock B.

Given the heterogeneity of the data, one might imagine that time series from different stocks should be first normalized (by average daily volume, average price or volatility etc.) before pooling them. Surprisingly, this appears not to be the case: we have observed that standard data transformations such as normalizations based on average volume, volatility or average spread, or partitioning the training data into sectors or categories such as large/small tick stocks do not improve training results. For example, a deep learning model trained on small tick stocks does not outperform the universal model in terms of forecasting price moves for small tick stocks. It appears that the model arrives at its own data-driven normalization of inputs based on statistical features of the data rather than ad hoc criteria. The source of the universal model's outperformance is well-demonstrated by Figure 7. The universal model most strongly outperforms the stock-specific models on stocks with less data. The stock-specific model is more exposed to overfitting due to the smaller dataset while the universal model is able to generalize by interpolating across the rich scenario space of the pooled data set and therefore is less exposed to overfitting. So, the existence of these common features seems to argue for pooling the data from different stocks, notwithstanding their heterogeneity, leading to a much richer and larger set of training scenarios. Using 1 year of the pooled data set is roughly equivalent to using 500 years (!) of data for training a single-stock model and the richness of the scenario space is actually enhanced by the diversity and heterogeneity of behavior across stocks.

Due to the large amount of data, very large universal models can be estimated without overfitting. Figure 8 shows the increase in accuracy for a universal model with 150 units per layer (which amounts to several hundred thousand parameters) versus a universal model with 50 units per layer. Remarkably, the universal model is even able to generalize to stocks which were not part of the training sample: if the model is only trained on data from stocks {1, . . . , N }, its forecast accuracy is similar for stock N +1. This implies that the universal model is capturing features in the relation between order flow and price variastions which are common to all stocks. Table 1 illustrates the forecast accuracy of a universal model trained only on stocks 1-464 (for January 2014-May 2015), and tested on stocks 465-489 (for June-August 2015). This universal model outperforms stock-specific models for stocks 465 -489, even though the universal model has never seen data from these stocks in the training set. The universal model trained only on stocks 1 -464 performs roughly the same for stocks 465 -489 as the universal model trained on the entire dataset of stocks 1 -489. Results are reported in Table 1.

Figure 10 displays the accuracy of the universal model for 500 completely new stocks, which are not part of the training sample. The universal model achieves a high accuracy on these new stocks, demonstrating that it is able to generalize to assets that are not included in the training data. This is especially relevant for applications, where missing data issues, stock splits, new listings and corporate events constantly modify the universe of stocks.

Stationarity

The relationships uncovered by the deep learning model are not only stable across stocks but also stationary in time. This is illustrated by examining how forecast accuracy behaves when the training period and test period are separated in time.

Figure 10 shows the accuracy of the universal model on 500 stocks which were not part of the training sample. The left histogram displays the accuracy in June-August, 2015, shortly after the training period (January 2014-May 2015), while the right plot displays the cross-sectional distribution of accuracy for the same model in January-March, 2017, 18 months after the training period. Interestingly, even one year after the training period, the forecasting accuracy is stable, without any adjustments.

Such stability contrasts with the common practice of 'recalibrating' models based on a moving window of recent data due to perceived non-stationarity. If the data were nonstationary, accuracy would decrease with the time span separating the training set and the prediction period and it would be better to train models only on recent periods immediately before the test set. However, we observe that this is not the case: Table 2 reports forecast results for models trained over periods extending up to 1, 3, 6, and 19 months before the test set. Model accuracy consistently increases as the length of the training set is increased. The message is simple: use all available data, rather than an arbitrarily chosen time window.

Note that these results are not incompatible with the data itself being non-stationary. The stability we refer to is the stability of the relation between the inputs (order flow and price history) and outputs (forecasts). If the inputs themselves are non-stationary, the output will be non-stationary but that does not contradict our point in any way.

Path-dependence

Statistical modeling of financial time series has been dominated by Markovian models which, for reasons of analytical tractability, assume that the evolution of the price and other state variables only depends on their current value and there is no added value to including their history beyond one lag. There is a trove of empirical evidence going against this hypothesis, and pointing to long-range dependence in financial time series [START_REF] Bacry | Continuous cascade models for asset returns[END_REF][START_REF] Lillo | The long memory of the efficient market[END_REF][START_REF] Mandelbrot | The multifractal model of asset returns[END_REF]. Our results are consistent with these findings: we find that the history of the limit order book contains significant additional information beyond that contained in its current state.

Figure 11 shows the increase in accuracy when using an LSTM network, which is a function of the history of the order book, as compared with a feedforward neural network, which is only a function of the most recent observation (a Markovian model). The LSTM network, which incorporates temporal dependence, significantly outperforms the Markovian model.

The accuracy of the forecast also increases when the network is provided with a longer history as input. Figure 12 displays the accuracy of the LSTM network on a 5, 000-step sequence minus the accuracy of the LSTM network on a 100-step sequence. Recall that a step ∆ k = τ k+1 -τ k is on average 1.7 seconds in the dataset so 5000 lags corresponds to 2 hours on average. There is a significant increase in accuracy, indicating that the deep learning model is able to find relationships between order flow and price change events over long time periods.

Our results show that there is significant gain in model performance from including many lagged values of the observations in the input of the neural network, a signature of significant -and exploitable -temporal dependence in order book dynamics.

Discussion

Using a Deep Learning approach applied to a large dataset of billions of orders and transactions for 1000 US stocks, we have uncovered evidence of a universal price formation mechanism relating history of the order book for a stock to the (next) price variation for that stock. More importantly, we are able to learn this mechanism through supervised training of a deep neural network on a high frequency time series of the limit order book. The resulting model displays several interesting features: Figure 11: Comparison of out-of-sample forecast accuracy of a LSTM network with a feedforward neural network trained to forecast the direction of next price move based on the current state of the limit order book. Cross-sectional results for 500 stocks for test period June-August, 2015.

Figure 12: Out-of-sample increase in accuracy when using a 5000-step sequence versus a 100-step sequence, across 1, 000 stocks. Test period : June-August 2015.

• Universality: the model is stable across stocks and sectors, and the model trained on all stocks outperforms stock-specific models, even for stocks not in the training sample, showing that features captured are not stock-specific.

• Stationarity: model performance is stable across time, even a year out of sample.

• Evidence of 'long memory' in price formation: including order flow history as input, even up to several hours, improves prediction performance.

• Generalization: the model extrapolates well to stocks not included in the training sample. This is especially useful since it demonstrates its applicability to recently listed instruments or those with incomplete or short data histories.

Our results illustrate the applicability and usefulness of Deep Learning methods for modeling of intraday behavior of financial markets. In addition to the fundamental insights they provide on the nature of price formation in financial markets, these findings have practical implications for model estimation and design. Training a single universal model is orders of magnitude less complex and costly than training or estimating thousands of single-asset models. Since the universal model can generalize to new stocks (without training on their historical data), it can also be applied to newly issued stocks or stocks with shorter data histories.

Figure 2 :

 2 Figure 2: Architecture of a recurrent neural network.

Figure 4 :

 4 Figure 4: Comparison of a deep neural network with linear models. Models are trained to predict the direction {-1, +1} of next mid-price move. Comparison for approximately 500 stocks and out-of-sample results reported for June-August, 2015. Left-hand figure: increase in accuracy of stock-specific deep neural networks versus stock-specific linear models. Righthand figure: accuracy of a universal deep neural network (red) compared to stock-specific linear models (blue).

Figure 5 :

 5 Figure 5: Left: relation between depth at the bid, depth at the ask and the probability of a price decrease. The x-axis and y-axis display the quantile level corresponding to the observed bid and ask depth. Right: Contour plot displaying the influence of levels deeper in the order book (5 to 10) on the probability of a price decrease.

Figure 6 :

 6 Figure 6: Out-of-sample forecasting accuracy of the universal model compared with stockspecific models. Both are deep neural networks with 3 LSTM layers followed by a ReLU layer. All layers have 50 units. Models are trained to predict the direction of the next move. Comparison across 489 stocks, June-August, 2015.

Figure 7 :

 7 Figure 7: Increase in out-of-sample forecast accuracy (in %) of the universal model compared to stock-specific model, as a function of size of training set for stock specific model (normalized by total sample size, N = 24.1 million). Models are trained to predict the direction of next price move. Comparison across 500 stocks, June-August, 2015.

Figure 8 :

 8 Figure 8: Comparison of two universal models: a 150 unit per layer model versus 50 unit per layer model. Models are trained to predict direction {-1, +1} of next mid-price move. Out-of-sample prediction accuracy for direction of next price move, across approximately 500 stocks (June-August, 2015).

Figure 9 :

 9 Figure 9: Performance on approximately 500 new stocks which the model has never seen before. Out-of-sample accuracy reported for June-August, 2015. Universal model trained during time period January 2014-May 2015.

Figure 10 :

 10 Figure 10: Performance on 500 new stocks which the model has never seen before. Left: outof-sample accuracy reported for June-August, 2015. Right: out-of-sample accuracy reported for January-March, 2017. Universal model trained on data from January 2014-May 2015.

Table 1 :

 1 Comparison of universal model trained on stocks 1-464 versus (1) stock-specific models for stocks 465-489 and (2) universal model trained on all stocks 1-489. Models are trained to predict direction of next mid-price move. Second column shows the fraction of stocks where the universal model trained only on stocks 1-464 outperforms models (1) and (2). The third column shows the average increase in accuracy. Comparison for 25 stocks and out-of-sample results reported for June-August, 2015.

	Model	Comparison Average increase in accuracy
	Stock-specific	25/25	1.45%
	Universal	4/25	-0.15%

Table 2 :

 2 Out-of-sample forecast accuracy of deep learning models trained on entire training set (19 months) vs. deep learning models trained for shorter time periods immediately preceding the test period, across 50 stocks Aug 2015. Models are trained to predict the direction of next price move. Second column shows the fraction of stocks where the 19-th month model outperforms models trained on shorter time periods. The third column shows the average increase in accuracy across all stocks.

Historical order book data was reconstructed from NASDAQ Level III data using the LOBSTER data engine[START_REF] Huang | LOBSTER: Limit Order Book Reconstruction System[END_REF].

the London Quant Summit 2018, JP Morgan and Princeton University for their comments. Computations for this paper were performed using a grant from the CFM-Imperial Institute of Quantitative Finance and the Blue Waters supercomputer grant "Distributed Learning with Neural Networks".


Anne Blavette
email: anne.blavette@ens-rennes.fr

Thibaut Kovaltchouk
email: thibaut.kovaltchouk@ac-reims.fr

François Rongère
email: francois.rongere@ec-nantes.fr

Marilou Jourdain De Thieulloy

Paul Leahy
email: paul.leahy@ucc.ie

Bernard Multon
email: bernard.multon@ens-rennes.fr

H Ben Ahmed
email: hamid.benahmed@ens-rennes.fr

Marilou Jourdain

Hamid Ben

Influence of the wave dispersion phenomenon on the flicker generated by a wave farm

Keywords: Flicker, aggregation effect, hydrodynamic modelling, time delay-based approach

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Influence of the wave dispersion phenomenon on the flicker generated by a wave farm I. INTRODUCTION The inherently fluctuating nature of waves may be reflected to some extent in the power output of wave energy converters (WECs). These fluctuations can induce voltage fluctuations which can potentially generate flicker [START_REF] Molinas | Power Smoothing by Aggregation of Wave Energy Converters for Minimizing Electrical Energy Storage Requirements[END_REF]- [START_REF] Kovaltchouk | Wave farm flicker severity: Comparative analysis and solutions[END_REF] Hence, wave farm managers will be required to demonstrate that their farm is compliant with grid codes and similar regulations, in order to be granted grid connection. This is usually performed through grid impact assessment studies by means of numerical power system simulators such as DIgSILENT PowerFactory [START_REF]DIgSILENT PowerFactory[END_REF], PSS®E [START_REF] Pss®e | [END_REF], etc. Hence, numerical models of the considered WEC(s) are necessary.

These models can be based on experimental data in the form of electrical power output time series [START_REF] Blavette | Impact of a Medium-Size Wave Farm on Grids of Different Strength Levels[END_REF] or on socalled "wave-to-wire" models which compute this type of data from the (usually simulated) sea surface elevation. A significant number of such "wave-to-wire" models have been developed, as reviewed in [START_REF] Penalba | A Review of Wave-to-Wire Models for Wave Energy Converters[END_REF]. However, few of them have considered arrays of WECs as described in [START_REF] Forehand | A Fully Coupled Wave-to-Wire Model of an Array of Wave Energy Converters[END_REF]. Regarding these latter, different approaches have been used. The most comprehensive approach consists in simulating the sea surface elevation at each node of the wave farm where a WEC is located, taking into account the wave dispersion phenomenon, as well as the hydrodynamic interactions between WECs due to radiation and diffraction. Using this approach is very heavy from a computational perspective and should be restricted to simulating the output power of a wave farm whose WECs are closely located. In the case where the WECs are sufficiently far away from each other so that their hydrodynamic interactions can be considered as negligible, a second approach should be used which consists in calculating the sea surface elevation at each node of the wave farm where a WEC is located without taking into account the radiation and diffraction due to the neighbouring WECs. Finally, a third simplified approach has been widely used in the electrical engineering community in flicker studies focussing on wave farms. This approach consists in calculating the output power of an entire farm based on the power profile of a single WEC. This power profile serves as a reference to which a random time delay is applied for each WEC in order to model the device aggregation effect [1]- [START_REF] Kovaltchouk | Wave farm flicker severity: Comparative analysis and solutions[END_REF], which will be described in Section C.2. The computational effort regarding this latter approach is extremely light with respect to the other two approaches. However, questions remain concerning its physical validity, as this approach does not take into account the wave dispersion phenomenon.

The objective of this paper is to tackle this question by comparing the flicker level obtained from the second and the third approaches described in this section. Sufficiently distantly located WECs will be considered in order to neglect the inter-WECs hydrodynamic interactions. The farm output power is then injected in a local electrical grid model developed under PowerFactory to compute the corresponding voltage profile at the Point of Common Coupling (PCC). This voltage profile serves as input to a flickermeter from which the associated short-term flicker level is computed. The modelling hypotheses will be described in Section II and the results in Section III. In Section IV, the conclusions will be detailed.

The results of the comparative study will contribute in defining the required level of hydrodynamic detail necessary for simulating the output power of a wave farm when it is to be used for flicker analyses.

II. MODELLING HYPOTHESES

A. Hydrodynamic simulation

The hydrodynamic model is based on linear wave theory and simulates wave field from a superposition of Airy waves obtained through discretising a JONSWAP spectrum and using random phases. Contrary to the time delay-based method, the wave dispersion phenomenon is taken into account here.

By discretising the wave spectrum 𝑆(𝜔) using 𝑛 regularly spaced frequency components, the amplitude of each elementary wave component is given by [START_REF] Faltinsen | Sea loads on ships and offshore structures[END_REF] as:

𝑎(𝜔 𝑖 ) = 𝑎 𝑖 = √2𝑆(𝜔 𝑖 )δ𝜔 (1)
Each wave component then represents a complex elementary free surface elevation at horizontal position (𝑥, 𝑦) and time 𝑡: 𝜂 ̃𝑖(𝑥, 𝑦, 𝑡) = 𝑎 𝑖 𝑒 𝑗[𝑘 𝑖 (𝑥 cos 𝜃 ̅ +𝑦 sin 𝜃 ̅ ) -𝜔 𝑖 𝑡 + 𝜑 𝑖 ]

(2

)
where 𝜃 ̅ is the mean wave direction of the mono-directional wave field, 𝜑 𝑖 ∈ [0, 2𝜋[ is the random phase of the 𝑖 th wave component chosen at wave field initialisation and 𝑘 𝑖 is the wave number, which is solution of the dispersion relation given by

𝜔 𝑖 2 = 𝑔𝑘 𝑖 tanh 𝑘 𝑖 ℎ (3)
with ℎ being the mean water depth at position (𝑥, 𝑦). If the water depth can be considered as infinite, the relation degenerates to 𝜔 𝑖 2 = 𝑔𝑘 𝑖 . Summing the 𝑛 contributions gives the free surface elevation: This linear wave field modelling is then integrated into a linear framework for wave structure interactions that relies on hydrodynamics coefficients obtained from the linear potential flow theory. Note that in this study, no hydrodynamic interactions between the scattered wave fields are taken into account so that only coefficients for one isolated body are to be calculated using the seakeeping software NEMOH [START_REF] Babarit | Theoretical and numerical aspects of the open source BEM solver NEMOH[END_REF] which is based on the Boundary Element Method (BEM).

The time domain linear excitation force applying to a body having position (𝑥, 𝑦) is then obtained by the superposition of the excitation generated by each wave component as:

𝐹 𝑒𝑥 (𝑥, 𝑦, 𝑡) = 𝑅𝑒 {∑ 𝐹 ̃𝑒𝑥 (𝜔 𝑖 )𝜂 ̃𝑖(𝑥, 𝑦, 𝑡) 𝑛 𝑖=1 } (5)
Sea-states were simulated for significant heights equal to 1 m and 3 m, as well as for peak periods equal to 7 s, 9 s, 10 s and 12 s.

B. Wave device

The wave farm is composed of identical heaving buoys controlled passively and described in a previous paper [START_REF] Kovaltchouk | Influence of control strategy on the global efficiency of a Direct Wave Energy Converter with Electric Power Take-Off[END_REF]. As the focus of this paper is on the comparison of two methods for modelling the device aggregation effect on flicker, a simple, passive control strategy was adopted for the WEC. It consists of the application of a constant damping factor as a function of the sea-state characteristics (significant wave height 𝐻 𝑠 and 𝑇 𝑝 ). This damping factor is optimised with respect to a given sea-state during a preliminary offline study. For the sake of realism, levelling is applied on the power takeoff (PTO) force, which is limited to 1 MN, and on the output electrical power, which is limited to 1 MW. Each WEC is connected to the offshore grid through a fully rated back-toback power electronic converter.

C. Wave farm 1) Wave farm layout

The wave farm considered in this study is considered to be composed of 24 of the devices described in the previous section. All these devices are deemed identical in terms of hydromechanical and electro-mechanical properties. They are placed at a distance 𝑑 of each other, on 3 rows and 8 columns facing the incoming waves, as shown in Fig. 1. The inter-WEC distance 𝑑 is supposed to be sufficient so that the hydrodynamic interactions between the devices can be considered as negligible. In this research work, it was assumed equal to 600 m in the full hydrodynamic approach, while it is made approximately equal to 600 m in the timedelay based approach, as it will be described in Section C.4 [START_REF] Babarit | Impact of long separating distances on the energy production of two interacting wave energy converters[END_REF].

2) Introduction to the two approaches studied here

The wave farm power output is computed as the sum of the power output of all the WECs composing the farm. As the temporal profile of the sea surface elevation at a given node in the wave farm is not expected to be identical to this corresponding to another node, it is not expected either that two power output temporal profiles from two different WECs could be identical. Hence, the power output of a wave farm cannot be computed as the product of the power output of a single device times the number of WECs composing the farm. Also, the fact that different WECs achieve peak power at different times leads to a reduced peak-to-average ratio of the wave farm power output compared to that of a single device. This is illustrated in Fig. 2 which shows the temporal power output profile for a single WEC and for the wave farm composed of 24 devices normalised by their respective average value. While the peak-to-average ratio is equal to 3.6 for the single WEC (even though its power output is limited to 1 MW), it is equal to 1.9 for the wave farm. This decrease in the peak-to-average ratio implies that the temporal power output profile is "smoother" in the case of a wave farm than in the case of a single WEC which is usually referred to as the device aggregation effect. The objective of this paper was to determine whether a full hydrodynamic simulation was required to model this effect on flicker, or whether a simplified, time delay-based method was sufficient. These two approaches are described in the next sections. Fig. 2 Temporal power output profile (over 100s) for a single WEC (blue) and for the wave farm composed of 24 devices (pink) for significant wave height 𝐻 𝑠 =3 m and peak period 𝑇 𝑝 =7 s. The profiles are normalised with respect to their average value.

2.1) Full hydrodynamic simulation

In this approach, the wave excitation force at each node of the wave farm where a WEC is located is computed by means of the code described in Section A. Then, the power output of each WEC is computed based on its corresponding excitation force temporal profile. Hence, the power output 𝑃 𝑓𝑎𝑟𝑚 of the wave farm corresponds to the algebraic sum of the power output 𝑃 𝑊𝐸𝐶𝑖 of each WEC 𝑖, such as:

𝑃 𝑓𝑎𝑟𝑚 (𝑡) = ∑ 𝑃 𝑊𝐸𝐶𝑖 (𝑡) 24 𝑖=1 (7)

2.2) Time delay-based method

The time delay-based method requires only a single power output temporal profile 𝑃 𝑊𝐸𝐶1 of a single WEC, to which different uniformly distributed random time delays ∆𝑡 𝑖 are applied to represent the power output of the other devices composing the farm. Hence, the wave farm output power can be expressed as: 𝑃 𝑓𝑎𝑟𝑚 (𝑡) = ∑ 𝑃 𝑊𝐸𝐶1 (𝑡 + ∆𝑡 𝑖 ) 24 𝑖=1 [START_REF]DIgSILENT PowerFactory[END_REF] where ∆𝑡 𝑖 models the fictive propagation of a wave group whose envelope characteristics are independent of the travelled distance. This means that the wave dispersion phenomenon is not taken into account here. This physical effect implies that in dispersive media such as water, the travel speed of a sine wave is linked to its frequency through the dispersion relationship which, for deep water waves, can be expressed as [START_REF] Falnes | Ocean waves and oscillating systems[END_REF]:

𝜔 2 = 𝑔𝑘 (9)
as mentioned earlier. Term 𝛥𝑡 𝑖 is assumed equal to:

𝛥𝑡 𝑖 = 𝑑 𝑡𝑑 𝑣 𝑔 = 4𝜋𝑑 𝑡𝑑 𝑔𝑇 𝑝 ( 10
)
where 𝑑 𝑡𝑑 is distance between a reference WEC (whose time delay is equal to zero) and given WEC 𝑖, and 𝑣 𝑔 is the group speed which is defined as equal to 𝑔𝑇 𝑝 /4𝜋 here, where 𝑔=9.81 m.s² is the gravity of Earth. Given that the incoming waves are simulated as mono-directional waves, the distance 𝑑 𝑡𝑑 taken into account here is equal to the distance along the axis parallel to the wave front propagation direction. In the full hydrodynamic approach, the distance between the WECs was assumed to be equal to a fixed distance 𝑑 =600 m. However, if this constant distance were used in the time delay-based approach (i.e. 𝑑 𝑡𝑑 = 𝑑), given that the excitation force temporal profile is similar for all WECs, then all the devices located on a given row of the wave farm would present the same power output at any time 𝑡, thus resulting in coincident power profiles for 8 WECs, which is unrealistic. Hence, in order to avoid this situation, an additional uniformly distributed random distance 𝑑 𝑟𝑎𝑛𝑑 , arbitrarily selected as ranging between -50 m and +50 m (in order to represent WECs linear drift), is added to the fixed inter-WEC distance 𝑑 such as:

𝑑 𝑡𝑑 = 𝑑 + 𝑑 𝑟𝑎𝑛𝑑 where 𝑑 𝑟𝑎𝑛𝑑 ∈ [-50; 50] m (11)
Ten time delay sets were used in this study.

D. Electrical grid

An electrical grid model was developed under the power system simulator PowerFactory and is shown in Fig. 3. This model is inspired from the Atlantic Marine Energy Test Site (AMETS) [START_REF] Amets | SEAI website[END_REF] located off Belmullet, Ireland for the onshore local grid part. It is composed of a 10/20kV transformer whose impedance is equal to 2.10 -4 +𝑗0.06 pu (where 𝑗 is the imaginary unit) and of a 0.1 MW load representing the onshore substation connected to the rest of the national network through a 5 km-long overhead line of impedance 0.09+𝑗0.3 Ω/km. On the 20 kV bus (which is the Point of Common Coupling (PCC)), a VAr compensator maintains power factor at unity. Then, a 20/38 kV transformer (of impedance equal to 2.10 -4 +𝑗0.06 pu) connects to the farm to the local (national) grid where a 2 MW load (representing the consumption of a local town) is also connected. The rest of the national grid is modelled by means of a 38 kV voltage source in series with an impedance. This impedance magnitude is selected to be equal to Z=20 Ω (i.e. equal to a short-circuit level of 72 MVA), and its angle is selected to be equal to 30°, which corresponds to a "weak grid" and constitutes thus a worst case scenario in which relatively high flicker levels can be expected. The offshore grid is composed of a 20 km-long submarine cable of series impedance equal to 0.07+ 𝑗 0.11 Ω/km and capacitance equal to 0.31 µF/km, as described in [START_REF]Nexans data sheet[END_REF]. The cable distance was selected according to the values observed for two planned or already existing wave energy test sites [START_REF] Amets | SEAI website[END_REF], [START_REF][END_REF]. The offshore network is also composed of a 0.4/10 kV transformer (of impedance equal to 2.10 -4 +𝑗0.06 pu) and of 24 wave devices. The influence on the study results of the internal network between the WECs and the 0.4/10 kV transformer was deemed negligible and was therefore not included in the model.

III. RESULTS

A. Flicker level with respect to time delays

As mentioned earlier, ten different time delay sets were used in the time delay-based approach. The corresponding minimum, maximum, and average short-term flicker levels 𝑃 𝑠𝑡 are shown in Table I and Table II for the two significant wave heights considered here (𝐻 𝑠 =1 m and 𝐻 𝑠 =3 m). The standard deviation, also shown in these tables, indicates that for most cases the deviation from the average value is relatively small, compared to the allowed flicker limits which range usually between 0.35 and 1 [START_REF] Blavette | Impact of a Medium-Size Wave Farm on Grids of Different Strength Levels[END_REF]. However, some higher values of the standard deviation indicate that the time delay set can have a non-negligible influence on flicker, and that it is therefore important to average the flicker level corresponding to several time delay sets in order to obtain a reasonable estimation of the flicker which would have been obtained through the more realistic, full hydrodynamic approach. As the average flicker level is mostly representative of the order of magnitude of the flicker level corresponding to the 10 different time delay sets, it will be used for the comparative study between the time delay-based approach and the full hydrodynamic approach, as described in the following section.

B. Comparison of the two approaches

It is shown in Fig. 4 that both approaches generate similar results with a difference which is generally negligible in comparison with the usual maximum allowed flicker limits ranging between 0.35 and 1. This observation applies to both the low-energy and the mild sea-states (𝐻 𝑠 =1 m and 𝐻 𝑠 =3 m respectively).

Hence, it can be concluded that the flicker level generated by a wave energy farm can be estimated with a relatively high level of accuracy in most cases by means of the average flicker level corresponding to several time delay sets (here, ten time delay sets were used). In other words, this means flicker can be estimated from a single WEC power output, without further requirement for modelling the hydrodynamic conditions at each node in the farm where a WEC is expected to be located. This means also, in physical terms, that the wave dispersion phenomenon can usually be considered as negligible when it comes to flicker studies under the conditions considered in this study. Fig. 4 Short-term flicker level 𝑃 𝑠𝑡 as a function of the sea-state peak period 𝑇 𝑝 for 𝐻 𝑠 =1 m and 𝐻 𝑠 =3 m, and for the two considered approaches IV. CONCLUSIONS This paper has described a comparative study between a time delay-based approach and a more realistic, full hydrodynamic approach for determining the flicker level generated by a wave farm composed of 24 devices. The results have shown that in most cases, using the average flicker level corresponding to 10 different time delay sets leads to a negligible error compared to the full hydrodynamic approach. This means that the wave dispersion phenomenon has a limited impact on flicker. However, some non-negligible values for the flicker level error in some rare cases suggest that the time delay-based approach should be restricted to estimating flicker at a first stage before more refined studies based on the full hydrodynamic approach are conducted. Future work will focus on the comparative analysis, in terms of flicker level, between the two approaches described in this paper and a more comprehensive hydrodynamic approach including the hydrodynamic interactions between WECs. It will also investigate the influence of several parameters such as inter-WEC distance, WEC spatial arrangement, device number, etc.

Fig. 1

 1 Fig. 1 Wave farm spatial layout

Fig. 3

 3 Fig. 3 Electrical grid model developed under DIgSILENT PowerFactory

Table I

 I Short-term flicker levels 𝑃 𝑠𝑡 for the ten time delay sets (𝐻 𝑠 =1 m)

				Peak period 𝑇 𝑝 (s)
			7	9	10	12
	level 𝑃 𝑠𝑡 Short-term flicker	Average Standard deviation Minimum Maximum	0.09 0.04 0.04 0.14	0.40 0.12 0.26 0.68	0.50 0.14 0.31 0.84	0.76 0.16 0.60 1.14

Table II Short

 II -term flicker levels 𝑃 𝑠𝑡 for the ten time delay sets (𝐻 𝑠 =3 m)

					Peak period 𝑇 𝑝 (s)
				7	9	10	12
	Short-term	level 𝑃 𝑠𝑡 flicker	Average Standard deviation Minimum Maximum	0.63 0.12 0.43 0.80	0.80 0.09 0.67 0.94	0.96 0.21 0.73 1.40	0.83 0.10 0.73 1.09

ACKNOWLEDGMENT

The research work presented in this paper was partly conducted in the frame of the QUALIPHE project (ANR-11-PRGE-0013) funded by the French National Agency of Research (ANR) which is gratefully acknowledged.


Abdelghani Matine

Charles-Henri Bonnard

Anne Blavette

Salvy Bourguet

François Rongère

Thibaut Kovaltchouk

Emmanuel Schaeffer

Optimal sizing of submarine cables from an electro-thermal perspective

Keywords: Submarine cables, optimal sizing, electro-thermal, wave energy converter (WEC), finite element analysis (FEA)

à la diffusion de documents scientifiques de niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Optimal sizing of submarine cables from an electro-thermal perspective Abdelghani Matine #1 , Charles-Henri Bonnard #2 , Anne Blavette* 3 , Salvy Bourguet #4 , François Rongère +5 , Thibaut Kovaltchouk°6, Emmanuel Schaeffer #7 # IREENA, Université de Nantes, 37 boulevard de l'Université, 44602 Saint-Nazaire, France 1 abdelghani.matine@univ-nantes.fr , 2 charles-henri.bonnard@univ-nantes.fr, 4 salvy.bourguet@univ-nantes.fr, 7 emmanuel.schaeffer@univ-nantes.fr * SATIE (UMR 8029), Ecole Normale Supérieure de Rennes, avenue Robert Schuman, 35170 Bruz, France 3 anne.blavette@ens-rennes.fr + LHEEA (UMR 6598), Ecole Centrale de Nantes, 1 Rue de la Noë, 44300 Nantes, France 5 francois.rongere@ec-nantes.fr °Lycée F. Roosevelt, 10 rue du Président Roosevelt, 51100 Reims, France 6 thibaut.kovaltchouk@ac-reims.fr

Abstract-In similar fashion to most renewables, wave energy is capital expenditure (CapEx)-intensive: the cost of wave energy converters (WECs), infrastructure, and installation is estimated to represent 60-80% of the final energy cost. In particular, grid connection and cable infrastructure is expected to represent up to a significant 10% of the total CapEx. However, substantial economical savings could be realised by further optimising the electrical design of the farm, and in particular by optimising the submarine cable sizing. This paper will present the results of an electro-thermal study focussing on submarine cables temperature response to fluctuating electrical current profiles as generated by wave device arrays, and obtained through the finite element analysis tool COMSOL. This study investigates the maximum fluctuating current loading which can be injected through a submarine cable compared to its current rating which is usually defined under steady-state conditions, and is therefore irrelevant in the case of wave energy. Hence, using this value for design optimisation studies in the specific context of wave energy is expected to lead to useless oversizing of the cables, thus hindering the economic competitiveness of this renewable energy source.

I. INTRODUCTION

In similar fashion to most renewables, harnessing wave energy is capital expenditure (CapEx)-intensive: the cost of wave energy converters (WECs), infrastructure, and installation is estimated to represent 60-80% of the final energy cost [START_REF]Ocean Energy Strategic Roadmap 2016 -Building ocean energy for Europe[END_REF]. In particular, cable costs (excluding installation costs) are expected to represent up to a significant 10% of the total CapEx, based on the offshore wind energy experience [START_REF] On | Wind Turbine Technology and Operations Factbook[END_REF][START_REF] Macaskill | Offshore wind-an overview[END_REF]. However, substantial economical savings could be realised by further optimising the electrical design of the farm, and in particular by optimising the submarine cable sizing. Several studies have focussed on optimising the electrical network composed of the wave farm offshore network and/or of the local onshore network [START_REF] Nambiar | Optimising power transmission options for marine energy converter farms[END_REF][START_REF] Blavette | Dimensioning the Equipment of a Wave Farm: Energy Storage and Cables[END_REF]. In [START_REF] Nambiar | Optimising power transmission options for marine energy converter farms[END_REF], a techno-economic analysis was conducted on maximising the real power transfer between a wave energy farm and the grid by varying three design parameters of the considered array, including the export cable length, but excluding its current rating which was thus not considered as an optimisation variable. It seems that this latter parameter was selected as equal to the maximum current which may theoretically flow through the cable. This theoretical value was obtained from the maximum theoretical power output supposedly reached for a given sea-state and which was extracted from a power matrix. Following this, the corresponding scalar value for the maximum current for a given sea-state was obtained by means of load flow calculations based on a pre-defined electrical network model. Hence, the cable rating was assumed to be equal to the maximum value among the different current scalar values corresponding to several sea-states. In similar fashion, in [START_REF] Beels | A methodology for production and cost assessment of a farm of wave energy converters[END_REF] where a power transfer maximisation is also conducted, the cable rating is determined as well based on the maximum power output by a wave device array for different sea-states. In both these papers, there is no limitation on the considered sea-states from which energy can be extracted, apart from the operational limits of the wave energy device itself in [START_REF] Nambiar | Optimising power transmission options for marine energy converter farms[END_REF]. In other words, as long as the sea-state characteristics are compatible with the wave energy device operational limits in terms of significant wave height and period, it is considered that wave energy is harnessed. Another approach challenging this idea was proposed in [START_REF] Sharkey | Maximising value of electrical networks for wave energy converter arrays[END_REF] based on the offshore wind energy experience [START_REF] Crown | Round 3 Offshore Wind Farm Connection Study[END_REF]. This approach consists in rating the submarine export cable to a current level less than the maximum current which could be theoretically harnessed when the wave device operational constraints only are considered. The rationale underpinning this approach consists in considering that the most highly energetic seastates contribute to a negligible fraction of the total amount of energy harnessed every year. Hence, this corresponds to a negligible part of the annual revenue. However, harnessing wave energy during these highly energetic sea-states leads to an increased required current rating for the export cable whose associated cost is expected to be significantly greater than the corresponding revenue. Consequently, it seems more reasonable, from a profit maximisation perspective, to decrease the export cable current rating, even if it means shedding a part of the harnessable wave energy. However, in similar fashion to the papers mentioned previously, current is calculated in this paper as a scalar value representing the maximum level which can be reached during a given sea-state. In other words, the fluctuating nature of the current profile during a sea-state is not considered. However, the maximum current value, from which the cable current rating is usually calculated, flows in the cables during only a fraction of the sea-state duration. Based on a very simple model, it was shown in [START_REF] Blavette | Dimensioning the Equipment of a Wave Farm: Energy Storage and Cables[END_REF] that the slow thermal response of the cable (relatively to the fast current fluctuations generated from the waves) [START_REF] Adapa | Dynamic thermal ratings: monitors and calculation methods[END_REF][START_REF] Hosek | Dynamic thermal rating of power transmission lines and renewable resource[END_REF] leads to temperature fluctuations of limited relative amplitude compared to the current fluctuations. Hence, this implies that it could be feasible to inject a current profile whose maximum value is greater than the current rating without exceeding the conductor maximum allowed temperature, which is usually equal to 90°C for XLPE cable [START_REF] Nexans | Submarine cable 10kV[END_REF][START_REF] Abb | XLPE Submarine Cable Systems Attachment to XLPE Land Cable Systems -User's Guide[END_REF][START_REF]Nexans data sheet[END_REF]. Downrating submarine cables in this manner, compared to rating them with respect to the maximum but transient, current level flowing through it, could lead to significant savings from a CapEx point of view. In this perspective, this paper presents a detailed study on the thermal response of a submarine cable subject to a fluctuating current as generated by a wave device array. Section II will describe the development of a finite element analysis (FEA) based on a 2D thermal model of a 20kV XLPE submarine cable and performed using commercial FEA software COMSOL. The thermal response of a submarine cable to the injection of a fluctuating current profile as generated by wave energy arrays is analysed in Section III. The objective of this study is to determine the maximum current loading which can be injected through a submarine cable without exceeding the conductor thermal limit (equal to 90°C here) and compare this value to the cable rating. As mentioned earlier, this latter value is usually defined under static conditions which are irrelevant in the case of wave energy. Hence, its use for design optimisation studies in this specific context is expected to lead to useless oversizing of the cables, thus hindering the economic competitiveness of wave energy.

II. THERMAL MODELLING OF THE SUBMARINE POWER CABLE

A. Cable design and characteristics

This study considers a 20 kV XLPE insulated power cable containing three copper conductors, each with a cross section of 95 mm² and having each a copper screen, as shown in I. The static current carrying capacity of the considered cable is equal to 290A and is calculated according to IEC standards 60287-1-1 [START_REF]Calculation of the current rating: Part 1-1 Current rating equations (100% load factor) and calculation of losses[END_REF] and 60287-2-1 [START_REF]Electric cables -Calculation of the current rating -Part 2-1: Calculation of Thermal Resistance[END_REF]. It is based on the following assumptions:

B. Thermal Model

This section describes the development of a 2D finite element analysis (FEA) of the submarine cable thermal model using commercial software COMSOL. In order to predict the temperature distribution within the cable, the heat transfer equation of thermal conduction in transient state is applied [START_REF] Long | Essential Heat Transfer. s.l[END_REF]:

( )
where ρ is the mass density (kg.m -3 ), is the specific heat capacity (J.kg -1 .K -1 ), T is the cable absolute temperature (K), K is the thermal conductivity (W.m -1 .K -1 ) and Q is a heat source (W.m -3 ).

The heat sources in cable installations can be divided into two generic groups: heat generated in conductors and heat generated in insulators. The losses in metallic elements are the most significant losses in a cable. They are caused by Joule losses due to impressed currents, circulating currents or induced currents (also referred to as "eddy currents").

The heat produced by the cable metallic components (namely conductors, sheath and armour) can be calculated based on equations provided in IEC standard 60287-1-1 [START_REF]Calculation of the current rating: Part 1-1 Current rating equations (100% load factor) and calculation of losses[END_REF]. First, the Joule losses W c of the conductor can be calculated by using the following formula:

( ( )) ( )
where R 20°C is the resistance of the cable conductor at 20°C (W/km), α is the constant mass temperature coefficient at 20°C (K -1 ), and I c (A) is the current density. Terms y p and y s are the skin effect factor and the proximity factor respectively. The sheath and armour losses (W s and W a respectively) can be calculated such as: where U 0 is the applied voltage, δ is the loss angle, and ω is the angular frequency. However, the heat produced in the insulating layers expected to be significant, compared to the heat produced by the metallic components, under certain high voltage conditions only. Finally, the boundary conditions for this model are illustrated in Fig. 2. The modelled region is 7 m deep (H=7 m). The two side boundaries A and B are placed sufficiently far away from the cable so that there is no appreciable change in temperature with distance in the xdirection close to the boundaries. The cable is placed in the middle of the modelled region with respect to the x-direction, and its length is equal to 10 m (this value was proved sufficient in [START_REF] Swaffield | Methods for rating directly buried high voltage cable circuits[END_REF] for meeting the zero heat flux boundary conditions for sides A and B).

The soil surface C is assumed to be at constant ambient temperature of 12 ºC. The vertical sides A and B of the model are assumed to have a zero net heat flux across them due to their distance from the heat source. The equation in the side A and B is defined as:

( )

where n is the unit vector normal to the surface. On side D, the thermal exchanges via convection between the sea bed and the seawater must be taken into account. The heat convection exchange is defined as:

( ) ( )
where h is the heat transfer coefficient, T out is the sea temperature and T the temperature of the upper boundary of the sea bed (side D).

C. Current temporal profiles injected in the cable

The wave farm is composed of 15 to 20 identical heaving buoys controlled passively and described in another paper [START_REF] Kovaltchouk | Influence of control strategy on the global efficiency of a Direct Wave Energy Converter with Electric Power Take-Off[END_REF]. Different power output temporal profiles for a single WEC were computed by means of several combined simulation programmes described in [START_REF] Blavette | Influence of the wave dispersion phenomenon on the flicker generated by a wave farm[END_REF]. The first programme computes the wave excitation force at a single location in the wave farm. Then, the temporal profile of the excitation force is injected into a wave device model to obtain the corresponding electrical power output profile from which the wave farm power output is calculated. In order to model the device aggregation effect, the power output profiles of the other WECs composing the farm are computed by shifting the power profile of a single device by a random time delay, as described in a paper mentioned earlier [START_REF] Blavette | Influence of the wave dispersion phenomenon on the flicker generated by a wave farm[END_REF]. These power profiles are then injected into an electrical grid numerical model which is shown in Fig. 3. This model has been developed under the power system simulator PowerFactory [START_REF]DIgSILENT PowerFactory[END_REF] and is described in more detail in [START_REF] Blavette | Influence of the wave dispersion phenomenon on the flicker generated by a wave farm[END_REF]. The components of the offshore grid are highlighted in [START_REF] Blavette | Influence of the wave dispersion phenomenon on the flicker generated by a wave farm[END_REF].

A. Model validation under steady-state conditions

A steady-state load current equal to 290A (i.e. equal to the steady-state capacity of the considered cable) is used for the analysis. If the model is valid, the conductor temperature should remain below its maximum allowed limit which is equal to 90°C. Table III summarizes the calculation of the different heat sources needed to solve the heat transfer problem according to the equations given in IEC standard 60287-1-1 [START_REF]Calculation of the current rating: Part 1-1 Current rating equations (100% load factor) and calculation of losses[END_REF]. The same heat fluxes are used as source terms for the FEA model, which allows to compare the calculated temperature with the two methods. Fig. 8 shows the meshing of the submarine cable and its surrounding area that need finer elements because these areas are the most important sections of the presented analysis. Then, the areas which are far away from the cables can be modelled with a coarser mesh. The steady-state temperature field distributions in the cable and in its environment are shown in Fig. 9 and Fig. 10 respectively. The calculation from the IEC standard leads to a temperature of 67°C for the copper cores and 39°C for the external sheath. It can be seen that the maximum temperature of the copper conductor resulting from the FEA is 75°C while the external sheath reached 44°C, a little bit higher but in the same order of magnitude than the IEC standards. Note that both methods return copper core temperature below the critical temperature of 90°C, which provide a safety margin with a normal current load of 290 A. Despite the higher temperatures resulting from the FEA, one can see this results comparison as a form of validation of the model, especially considering that IEC standard uses a simplified model to calculate the temperature, i.e. an electricequivalent circuit composed of thermal resistors and current sources. Hence, it leads us to conclude that the FEA model can be used to calculate the temperature of a submarine cable under fluctuating current as generated by a group of WECs.

B. Thermal response under a fluctuating current profile

This section describes the transient thermal response of the submarine cable to different current profiles as generated by an array of wave energy devices considering several sea states, as described in Table II. The objective of this study is to investigate the levels of current which can be transmitted through a submarine cable without the conductor exceeding the thermal limit of 90°C. For each simulated case, we consider the maximum value of the current and its percentage with respect to the continuous current rating of the cable, i.e. 290 A. Table IV shows that the maximum current of each current profile. It is important to highlight that these maximum currents can be far above the continuous current rating in all cases. Simulation results of such a thermal problem depend on the initial conditions. Hence, it is important to accurately define the initial thermal conditions of the surrounding soil. The simplest initial condition which can be defined is a uniform temperature field. The value of this initial thermal condition should correspond to the case where the cable is subject to a current load equal to the average of the fluctuating current profile to be applied afterwards. The role of this first phase of the simulation is to quickly bring the cable temperature close to the expected range within which it is expected to vary once the fluctuating current profile is applied, thus reducing the simulation time. We used enough sequential repetitions of the current depicted in Figs 4 to 7 to reach a simulation time of 100 ks, i.e. a duration that is necessary to reach close-toequilibrium conditions for the thermal problem.

Figs. 11, 12, 13 and 14 show the conductor thermal response versus time, for Cases 1 to 4 respectively. In these cases, the current maximal values are equal to 97% to 287% of the cable capacity. The maximum temperature does not exceed the allowed limit of 90°C for the first two cases. In other words, as the temperature is below the allowed limit of 90°C, the cable could be considered as overrated with respect to the considered current profiles. The third case presents good agreement between the sea state and the sizing of the cable (Fig. 13), as the temperature is close to the maximum allowed value of 90°C. Fig. 14 shows the simulation results for Case 4. In this case, the current maximal value is up to 287 % of the cable capacity but the temperature exceeds 90°C. Hence, the cable can be considered as underrated here. In summary, under the conditions considered in this study, it appears that the cable is able to carry a fluctuating current profile whose maximum value is approximately equal to two and a half times the cable steady-state current rating.

CONCLUSIONS

This paper describes the results of a study focusing on the electrothermal response of a submarine cable to fluctuating current profiles as generated by a wave farm under different conditions, in particular regarding sea-state characteristics and the number of devices composing the farm. It was shown that the cable temperature remained below the allowed limit (equal to 90°C for the considered cable) when the maximum value of the injected current profile is as high as about two and a half times the (steady-state) rated current. Hence, it could be possible to size a cable to be included in a wave farm to the half of the maximum current, rather than to 100% of this current while keeping a good margin of safety. This could lead to significant savings in the CapEx of wave farm projects, thus contributing to improve the economic competitiveness of wave energy.

ACKNOWLEDGMENT

The research work presented in this paper was conducted in the frame of the EMODI project (ANR-14-CE05-0032) funded by the French National Agency of Research (ANR) which is gratefully acknowledged. This project is also supported by the S2E2 innovation centre ("pôle de compétitivité").

Fig. 1

 1 Fig. 1 Cross section of the considered three-phase export cable as modelled under COMSOL Fig. 1. The sheaths are made of polyethylene and the bedding of polypropylene yarn. The surrounding armour is made of galvanized steel. This cable is currently installed in the SEM-REV test site located off Le Croisic, France and managed by Ecole Centrale de Nantes. Usual values regarding the cable material thermal properties, as provided in IEC standard 60853-2 [15], are presented in TableI.

Fig. 2

 2 Fig. 2 Illustration of the cable environment and of its boundary conditions.-Maximal allowed conductor temperature at continuous current load : 90°C -Current frequency : 50 Hz -Ambient temperature : 12°C -Cable burial depth : 1.5 m -Thermal resistivity of surroundings: 0.7 K.m/W

Fig. 3

 3 Fig. 3 Wave farm electrical network model as developed under PowerFactory λ λ where λ 1 and λ 2 are dissipation factors in the sheath and in the armour respectively. Insulating materials also produce heat. Dielectric losses in the insulation are given by:

Fig. 4 Fig. 5 Fig. 6 Fig. 7 Fig. 8

 45678 Fig. 4 Current profile flowing through the cable (Case 1)

Fig. 9

 9 Fig. 9 Steady-state temperature field simulation (°C) of the submarine cable under normal load conditions.

Fig. 10

 10 Fig. 10 Steady-state temperature field simulation (°C) of the cable environment. The blue arrows represent the heat flux.

Fig. 13

 13 Fig.13 Cable temperature versus time (Case 3)

  It is composed of 15 to 20 WECs, of a 0.4/10 kV transformer, of a 6.5 km-long cable, and of a 10/20 kV transformer. The rest of the network is located onshore is shaded in blue in the figure. Four current temporal profiles were simulated for different sea-state characteristics (the significant wave height H s and the peak period T p ) and device numbers, as detailed in TableII. They are shown in Fig.4to 7.

TABLE III HEAT

 III LOSSES IN THE SUBMARINE CABLE

			Numerical
	Losses type	Formula	value
			(W/m)
	Conductor Joule losses W c		21.025
	Sheath losses W s		0.2
	Armour losses W a		0.027
	Dielectric losses W d		0.028