Anomaly Detection on Financial Time Series by Principal Component Analysis and Neural Networks

A major concern when dealing with financial time series involving a wide variety ofmarket risk factors is the presence of anomalies. These induce a miscalibration of the models used toquantify and manage risk, resulting in potential erroneous risk measures. We propose an approachthat aims to improve anomaly detection in financial time series, overcoming most of the inherentdifficulties. Valuable features are extracted from the time series by compressing and reconstructingthe data through principal component analysis. We then define an anomaly score using a feedforwardneural network. A time series is considered to be contaminated when its anomaly score exceeds agiven cutoff value. This cutoff value is not a hand-set parameter but rather is calibrated as a neuralnetwork parameter throughout the minimization of a customized loss function. The efficiency of theproposed approach compared to several well-known anomaly detection algorithms is numericallydemonstrated on both synthetic and real data sets, with high and stable performance being achievedwith the PCA NN approach. We show that value-at-risk estimation errors are reduced when theproposed anomaly detection model is used with a basic imputation approach to correct the anomaly.


Introduction
In the context of financial risk management, financial risk models are of utmost importance in order to quantify and manage financial risk. Their outputs, risk measurements, can both help in the process of decision making and ensure that the regulatory requirements are met (Basel Committee on Banking Supervision, 2013). Financial management thus heavily relies on financial risk models and the interpretation of their outputs. The data usually consist of time series representing a wide variety of market risk factors. A major issue of such data is the presence of anomalies. A value of the time series is considered to be abnormal whenever its behaviour is significantly different from the behaviour of the rest of the time series (Hawkins, 1980). In this work, we focus on the detection of abnormal observations in a market risk factor time series used to calibrate financial risk models. Indeed, erroneous input data may wrongly impact the calibrated model parameters. A typical example is the estimation of the covariance matrix of a bank market risk factors. This covariance matrix is involved in the computation of value-at risks (VaR) or expected shortfalls (expected losses above the VaR level). Since the true covariance is unknown, it has to be estimated from the data. However, the presence of anomalies in the data might have an impact on this estimation. For this specific case, robust methods, less sensitive to anomalies, can be used. Yet, existing robust estimators are computationally expensive, with a polynomial or even exponential time complexity in terms of the number of market risk factors. A faster approach was suggested by Cheng et al. (2019), but this algorithm only applies to the estimation of the covariance matrix in the case of Gaussian distribution. Financial risk models used by banks are widespread and various. Therefore, instead of seeking for a robust version for each of them, we propose to detect anomalies directly on the time series.
Contributions of the Paper Our PCA NN methodology identifies anomalies on time series using a principal component analysis (PCA) and neural networks (NN) as underlying models. This methodology overcomes common pitfalls associated with anomaly detection. PCA is first used as a features extractor on (augmented if needed) data. Anomaly detection is then performed in two steps. The first step identifies the time series with anomalies evaluating the propensity of the time series to be contaminated, as reflected by its so-called anomaly score. Toward this end, we calibrate a feedforward neural network through the minimization of a customized loss. This customized loss allows us to calibrate the cut-off value on the anomaly scores, without resorting to expert judgement. In this way, we remove the expert bias. The second step localizes the anomaly among the observed values of the identified contaminated time series.
Outline The PCA NN approach is detailed in Section 2. Section 3 describes the methodology used for generating the data on which our approach is thoroughly tested in Sections 4, 5 and 6, exploiting the knowledge of the data generating process for benchmarking purposes. Section 7 illustrates the benefit of using our approach as a data cleaning preprocessing stage on a downstream task, namely value-at-risk computations. This is completed by further numerical experiments on real data sets in Section 8. Section 9 concludes. Section A and B provide reviews of the anomaly detection literature and algorithms. Section C addresses the data stationarity issue.

Notations
Let X = X 1 , X 2 , . . . , X i , . . . , X n , where the column-vector X i = x i t 1 , x i t 2 , . . . , x i t j . . . , x i tp ∈ R p corresponds to the i-th observed time series, and x i t j to the value observed at time t j in the i-th time series. Our method fits into the supervised framework. We thus assume that the data matrix X comes along with two label vectors. The first label vector A = A 1 , . . . , A n ∈ {0, 1} n identifies the time series containing anomalies referred to as contaminated time series. For i = 1, . . . , n, we define the identification labels as A i = 1 if there exists j such that x i t j is an anomaly. 0 otherwise. (2) The second label vector L concerns solely the contaminated time series. Its coefficients correspond to the localization labels, i.e. time stamps at which an anomaly occurs. For j = 1, . . . , p and i ∈ I c = {i : A i = 1} (the set of contaminated time series), L i =  when the anomaly occurs at time t  for the i-th contaminated time series, X, A and L are general notations. As the model involves a learning phase, we denote by X T rain , A T rain and L T rain the data used for calibration and, by X T est , A T est and L T est the independent data sets on which the model performance is evaluated.

Methodology Overview
The anomaly detection model we propose is a two step supervised-learning approach. The contaminated time series identification step aims at identifying among the observed time series the ones with potential anomalies. The anomaly localization step consists in finding the location of the anomaly in each contaminated time series identified during the first step.
The model used in the first step falls under the scope of binary classification models, assigning to each time series a predicted label A i = 1 if X i is considered as contaminated by the identification model. 0 otherwise.
The second step uses a multi-classification model. For each contaminated time series identified in the first step, the model predicts a unique time stamp L i at which the anomaly has occurred. Formally, the observed value at the time stamp ti of the i-th contaminated time series is considered abnormal, meaning that the observation x i t i is abnormal. Thus, Our two-step anomaly detection model localizes only one anomaly per run, if any. Several iterations allow removing all anomalies. Indeed, as long as the time series is identified as contaminated, by the first step of the method, the time series can go through the second step. Once all anomalies have been localized the final stage of the approach is to remove the abnormal values and suggest imputation values (see Section 6.1).

Theoretical Basics of Principal Component Analysis
A fundamental intuition behind our anomaly detection algorithm is the existence of a lower dimensional subspace where normal and abnormal observations are easily distinguishable. Since the selected principal components explain a given level of the variance of the data, they represent the common and essential characteristics to all observations. As one might reasonably expect, these characteristics mostly represent the normal observations. When the inverse transformation is applied, the model successfully reconstructs the normal observations with the help of the extracted characteristics, while it fails into reconstructing the anomalies.
In this section, we recall the theoretical tenets of the principal component analysis (PCA) useful for our purpose. The idea of PCA is to project the data from the observation space of dimension p into a latent space of dimension k, with k < p. This latent space is generated by the k directions for which the variance retained under projection is maximal.
The first principal component U is given by U = w X, where w is a to-be-determined vector of weights. Since the PCA aims at maximizing the variance along its principal components, w is sought for as arg max w Σw = λ. Following the same process with an additional constraint regarding the orthogonality between the principal components, one can show that the k first principal components correspond to the eigenvectors associated with the k dominant eigenvalues of Σ. The transfer matrix Ω X ∈ R k×p is defined by the eigenvectors of Σ associated with its k largest eigenvalues. This transfer matrix Ω X applied to any observation in the original observation space projects it into the latent space. Since the transformation is linear, reconstructing the initial observations from their equivalent in the latent space is straightforward. Using PCA on the data matrix X, the transfer matrix Ω X can be inferred. The projection of the observations into the latent space Z ∈ R n×k is then given by Z = XΩ X , while their reconstructed values are given by X = ZΩ X .
The reconstruction errors ε i for each time series defined by referred in the sequel as the new features, are our new representation of data, given as inputs to our models. Algorithm 1 summarizes how to compute these reconstruction errors. The two steps of the model use the same inputs, namely these reconstruction errors. Nevertheless, the processing these inputs undergo in each step differs.
Algorithm 1 Features extraction with PCA inputs : number of principal components k out of sample data X T est training data if not trained yet X T rain if not trained yet then Ω X T rain ← PCA(X T rain , k) end ε T rain = X T rain Ω X T rain Ω X T rain − I p ε T est = X T est Ω X T rain Ω X T rain − I p outputs : ε T rain reconstruction errors for the train set ε T est reconstruction errors for the test set

Contaminated Time Series Identification
The first step of the approach aims at identifying the time series with anomalies based on the reconstruction errors obtained using PCA and Algorithm 1. In this section, we detail the process these new features undergo. We assign to each reconstruction error an anomaly score which computation is handled by a neural network (NN) (see Goodfellow et al., 2016). We also numerically demonstrate that our NN has an advantage over a naive approach in terms of being able to assess more accurately the propensity of a time series of being contaminated.
Naive approach A naive and natural approach to define the anomaly score is to consider the 2 -norm of the reconstruction errors ε i = ε i 2 , i = 1, . . . , n. This anomaly score represents the propensity of a time series to be contaminated or not. The scores are first split into two sets depending on whether the corresponding time series is contaminated or not. We denote by ε c and ε u the set of contaminated and uncontaminated time series anomaly scores, i.e.
The density distribution function of each class of time series f u and f c , i.e. for both ε c and ε u , is then estimated using kernel density estimation (Węglarczyk, 2018). For l = {c, u}, let ε l be an i.i.d sample of n l observations from a population with unknown density f l . The corresponding kernel estimator is given byf where K is a kernel function and h is a smoothing parameter. The cut-off value s is then chosen as the intersection between the two empirical density functions, i.e.ŝ = arg min where η is a to be tuned precision level. The selected cut-off valueŝ represents the value of the score for which the area under the curve of the density of uncontaminated time series aboveŝ and the area under the curve of contaminated time series belowŝ are as small as possible. We expect it to correspond to a value exceeded by a relative low number of ε i with i / ∈ I c asf u (s) is expected to be flat around the cut-off value exceedance region, and exceeded only by few ε i with i ∈ I c .
As shown in Figure 1, the naive approach results in a non-negligible overlapping region between the two densities, both on the train set (left) and for the test set (right). This region represents the anomaly score associated with time series that could be either contaminated or uncontaminated. The uncertainty regarding the nature of the observation when its anomaly score belongs to this region is relatively high in comparison with observations whose anomaly scores lie on the extreme left-hand or right-hand side of the calibrated cut-off value. Moreover, because we are not able to provide a clear separation between the scores of uncontaminated and contaminated time series, we may expect the model to mislabel future observations, resulting in a high rate of false positives and true negatives.
Neural network approach In view of getting a clearer separation of the densities, we propose an alternative approach for the computation of the anomaly scores built upon the reconstruction errors. Let F denote the function that associates with each reconstruction error ε i its anomaly score. In the naive approach, F corresponds to the 2 -norm. Hereafter, F instead represents the outputs of a feedforward Neural Network (NN). The training of the corresponding NN aims at minimizing a loss function that reflects our ambition to construct a function F giving accurate anomaly scores. The network used for computing the anomaly scores is defined by where h i (ε) = W i · ε + b i + represents the computation carried in the i-th layer of the neural network with n i h hidden units. H stands for the number of layers of the NN, with ReLU activation applied element-wise. We finally estimate the labels from the outputs followinĝ The idea of the learning is to calibrate the weights W : . , H such that the NN is able to accurately assess to which extent a time series is contaminated, with a clear distinction between the anomaly scores assigned to uncontaminated time series and those assigned to the contaminated ones, while integrating the calibration of the cut-off value s as part of the learning. We denote by Θ := {W , b, s} the set of parameters to be calibrated through the NN training. To meet these needs, the loss we minimize during the learning is given by where is the binary cross-entropy, a well-known loss function classically used for classification problems.
To this first component of our loss function, we add two components that aim at downsizing the overlapping region between the density of anomaly score of both types of observations. In order to have a control on the latter, we consider AUCDensity u and AUCDentsity c , which correspond to the area under the curve of the probability density function of anomaly scores the model assigns to contaminated and uncontaminated observations, respectively, i.e The bounds of these integrals depend on the cut-off value s and define the region for which we want the probability density functions to be as small as possible, which allows us to estimatê s. The estimated probability density functionsf u A,F (ε) andf c A,F (ε) depend on F (ε), the scores assigned to contaminated time series identification model, and of the identification labels A. We describe with Algorithm 2, the scoring and cut-off value calibration achieved by the calibration of a feedforward network with the loss (5). Each update of Θ, AdamStep in Algorithm 2, is carried following the Adam optimization algorithm of (Kingma and Ba, 2014).
Algorithm 2 Scoring and cut-off calibration inputs : learning rate lr number of maximum iterations K kernel density estimator parameters K,h training data ε T rain , A T rain Initialize parameter Θ, L = ∞ and count index k = 0 while k < K do // Density estimation following (4) f u A T rain ,scores u = KernelDensityEstimator(K, h, scores u ) f c A T rain ,scores c = KernelDensityEstimator(K, h, scores c ) // Loss evaluation following (5) and (6)  Table 1 shows that with the NN approach the densities of scores assigned to each type of time series display the expected behaviours on the left-hand (right-hand) side of the cut-off value for contaminated (uncontaminated) time series, on both the train and test sets. Actually, the lower AUCDensity u (AUCDensity c ), the lower the number of uncontaminated (contaminated) time series to which are assigned anomaly scores above (below) the cut-off value, which prevents mislabelling.
Once the features ε and the calibrated NN are provided, the contaminated times series identification model is ready for use. This step is described by Algorithm 3. Algorithm 3 Contaminated time series identification model inputs : time series to analyze X i , calibrated model parameter Θ, ε i ← PCAFeaturesExtraction(X i ) ; // cf. Algorithm 1 in Section 2.3 score i ← F Θ (ε i ) A i ← 1 score i >ŝ outputs : identification labelÂ i .

Anomaly Localization Step
Once the time series containing an anomaly have been identified, the second step of our approach aims at localizing an abnormal observation among each contaminated time series. Again we use the reconstruction errors defined in (3) as inputs of this second step. The difference lies on the transformation these reconstruction errors undergo before labelling the different time stamp observations. Namely, we now consider the following transformation of the reconstruction errors: with |·| meant component wise. The model input F (ε i ) ∈ R p + is thus assigned to each time series X i ∈ R p . The time stamp of occurrence of the anomaly is given by Algorithm 4 recaps the anomaly localization model.
Localization εL Imputation X,L X no yes Figure 3: Flow chart of our two step anomaly detection model (PCA NN), depicting the process a time series X goes through.
Algorithm 4 Anomaly Localization Model inputs : time series X i to analyze.
Note that anomalies are not necessarily extrema. For this reason, the above PCA features extraction step is necessary. Building on the reconstruction errors, the model identifies these extrema-anomalies, but also abnormal observations which are not necessarily extrema. This subtlety of the nature of the observations underlines the importance of going further than taking the index of the highest observed value of the contaminated time series X i .
The flow chart of Figure 3 along with Algorithm 5 summarize our approach. When a time series X is given to our model, we suggest a new representation of X, namely ε, through a features engineering step involving PCA. The resulting representation ε i feeds the first component of the model, i.e. the identification step, which evaluates the time series propensity of being contaminated. The optimal parameter Θ of the identification model solves whereÂ = 1 F ((ΩΩ −Ip)X)>s and L Θ is the loss function defined in (5). Then, if the model considers the time series as contaminated, the localization model takes over to localize the abnormal value in X. For this second step of the model, the time stamp of occurrence of the anomaly is given by Finally, the anomaly is imputed. As our approach integrates the computation of a reconstruction of the time series, we could replace the anomaly with the corresponding reconstructed value. The model will be then able to detect the anomaly and suggest an imputation value. However, numerical tests described in Section 6.1 show that imputations with naive approaches perform better.
Algorithm 5 Anomaly Detection Model inputs : time series X i to analyze, calibrated identification model parameter Θ anomaly localization L i .

Data Generation Process
Since anomalies are rare by definition, real data sets are very imbalanced: the proportion of anomalies compare to normal observations is very small, making the learning phase of the model difficult. This has led us to first consider synthetic data for the model calibration. The corresponding data set is obtained through a three-step process including time series simulations, contamination and data augmentation. We point out that, in the data simulation, care was taken to ensure that the generated data sets stay realistic. In particular, only few anomalies were added to time series as described in Section 3.2. To this extent, we still face the problem of scarcity of anomalies in this synthetic framework, and we provide some preprocessing steps to sidestep this issue as well.

Data Simulation
In this section, we describe the model used to simulate the data. Recall that our primary motivation is to detect anomalies in financial time series. For that purpose, we consider share price sample paths generated through the Black and Scholes model, i.e. geometric Brownian motions. Under this framework, the share price S t is defined by where W represents a standard Brownian motion, µ the drift, σ the volatility of the stock and S 0 the initial stock price. Let N stocks S 1 , . . . , S N simulated simultaneously from this model, where is the time series of length T representing the (time discretized) path diffusion of the i-th stock. Each stock S i has its own drift µ i , volatility σ i and initial value S i t 0 . The paths parameters are selected randomly according to For the sake of realism, the Brownian motions driving the N -stocks are correlated. The contamination of the obtained time series by anomalies is described in the next section.

Time Series Contamination
We apply a quite naive approach to introduce anomalies into our time series. We introduce the same fixed number of anomalies n anom to each time series by applying a shock on some original values of the observed time series. Formally, the -th added anomaly is characterized by its location t  corresponding to the time stamp at which the anomaly has occurred, its shock δ  , the amplitude of the shock is given by |δ  | and its sign by sgn (δ  ). We denote by S a,i the time series resulting from the contamination of the i-th clean time series S i . For i ∈ {1, . . . , N }, let J i be the set of indices of the time stamps at which an anomaly occurs for the i-th time series. For  ∈ J i , the abnormal values are The location, sign and amplitude of the shocks are generated randomly according to uniform distributions: where ρ is an upper bound on the shock amplitude. We define the anomaly mask matrix R N ×T by setting, for i = 1, . . . , N , and j = 1, . . . , T , Under this framework, when we incorporate the anomalies to the clean time series driven by the geometric Brownian motion, we assign labels to the time series observations according to whether the values correspond to anomalies or normal observations. We thus provide the labels Y i t j associated with each value S a,i t j . Hence, for i = 1, . . . , N , and j = 1, . . . , T , Algorithm 6 recaps the time series contamination procedure, where AnomalyMask and GetLabels refer to the operators defined by (10) and (11).
Algorithm 6 Data Contamination inputs : amplitude range ρ, number of anomalies to add in time series n anom , outputs : time series with anomalies S a , labels associated with each value of the time series Y .

Data Augmentation
By definition, anomalies are rare events and thus represent only a low fraction of the data set. Yet the suggested approach needs an important training set for an efficient learning. To overcome this issue, we apply a sliding window data augmentation technique (Le Guennec et al., 2016). This method not only extends the number of anomalies within the data set, but also allows the model to learn that anomalies could be located anywhere in the time series. It consists in extracting N p = T − p + 1 sub-time series of length p of the initial observed time series of length T . Note that we have to split the data into train and test sets before augmentation to guarantee that the time series considered in the training set do not share any observation with the ones we use for the model evaluation. In view of simplification, we introduce the data augmentation process, without loss of generality, for S and Y . In practice this process has to be applied to S T rain ,Y T rain and S T est , Y T est separately.
For i = 1, . . . , N and q ∈ {1, . . . , N p }, the sub-time series S i,q and the associated labels Y i,q are defined by (see Algorithm 7) Each sub-time series S i,q+1 thus results from the shift forward in time of one observation of the previous sub-time series S i,q . The final data set X and the labels s Y are then defined as the matrices whose rows correspond to the sub-time series S i,q q=1,...,Np;i=1,...,N , and Y i,q q=1,...,Np;i=1,...,N , respectively, i.e.
With this data configuration, an observation refers to a time series S i,q obtained through the sliding window technique. Two observations may represent the same stock but on different time intervals.
Algorithm 7 Data Augmentation with sliding window technique inputs : window size p, data and labels S,Y ; Initialize empty slided time series and labels matrices X, s Y ; end outputs : resulting slided time series and labels X, s Y While the fact that the observations share the same values may be argued to wrongly impact the learning process, we point out that a real benefit can be drawn from this situation. Indeed, thanks to this sliding window technique, the number of anomalies is considerably increased. This technique also allows the model to learn that anomalies could be located anywhere in the time series, reducing the dependency on the event location (Um et al., 2017).
However, once we apply the sliding window technique, we do not only extend the number of contaminated time series: the number of uncontaminated time series is also increased. But the minority class (contaminated time series) has, at least, a significant number of instances, denoted by N c . In order to get a more balanced data set for the training set, we perform an undersampling, selecting randomly N c observations from the N u uncontaminated time series without any anomalies (RandomSampling in Algorithm 8). The resulting retained number of observations 2N c is now sufficient to train the model. The test set in turn is imbalanced. To sharpen the imbalanced characteristic of the test set, we specify a contamination rate r c , which corresponds to the rate of contaminated time series in the data set. The construction of the test set is described in Algorithm 8.

Algorithm 8 Time series selection
Inputs : slided data and labels X, s Y , contamination rate r c (for test set) outputs : time series X with at most one anomaly, corresponding labels for identification task A, corresponding labels for localization task L.
As mentioned in Section 2.2, Algorithm 5 is designed to predict the localization of only one anomaly, (if there is more than one anomaly in the time series, it should be run iteratively as detailed in Section 2.2). Here we only consider time series with at most one anomaly (without loss of generality).
Assumption 1 A time series S i,q contains at most one anomaly among all its observed values.
We assign the identification label A i,q (see Section 2.1) to each S i,q following the rule 1 if there is an anomaly among the observed values of S i,q 0 otherwise.
Regarding the localization labels, we recall that they only concern the time series with an anomaly, therefore L i,q is defined following With Algorithm 9, we give a rundown of the construction process of the identification and localization labels, namely A and L departing from s Y .
The supervised learning framework is adopted herein, since we have at our disposal labelled data.
The resulting matrices with the observed values of the time series in X the associated identification labels in A and localization labels in L constitute the data set used for our model calibration and evaluation.
Algorithm 9 Time series labelling Inputs : slided labels s Y ; outputs : corresponding labels for identification task A, corresponding labels for localization task L.

Model Evaluation: Setting the Stage
After introducing the relevant performance metrics used to assess the PCA NN performance, we describe the synthetic data used for the numerical experiments of Sections 5 to 7. We then briefly describe the process of the latent space dimension calibration.

Performance Metrics
Common methods to assess the performance of binary classifiers include true positive and true negative rates, as well as ROC (Receiver Operating Characteristics) curves, displaying the true positive rate against the false positive rate. These methods, however, are uninformative when the classes are severely imbalanced. In this context, F 1 -score and Precision-Recall curves (PRC) should be used (Brownlee, 2020;Saito and Rehmsmeier, 2015). Our F 1 -score (Chinchor and Sundheim, 1993;Van Rijsbergen, 1979) combines these two measures in a single index defined as The closer the F 1 -score to 1, the better the prediction model. A PRC displays the values of Precision and Recall as the cut-off s varies from 0 to 1. The PRC of a skillful model bows towards the point with coordinates (1, 1). The curve of a no-skill classifier is an horizontal line at some y-level proportional to the proportion of Positives in the data set. For a balanced data set this proportion is just 0.5 (Brownlee, 2020).
PRC and F 1 -score are complementary in our approach. The F 1 -score is used on the anomaly scores outcomes of the models, to identify the best configuration and the best model. A PRC is used to select the best cut-off value used in the prediction of the two classes, for each model.

Synthetic Data Set
The anomaly detection task is performed on N = 20 stocks simultaneously. The stock prices are diffused according to the Black and Scholes model. Each stock has its own drift and volatility and the 20 stocks are correlated, as described in Section 3. Each time series represents T = 1, 500 daily stock prices, split into two sets. The first 1, 000 observations, corresponding to the train set, are used to learn the model parameters, i.e. the PCA transfer matrix and the NN weights and biases. The last 500 observations, corresponding to the test set, are used to assess the quality of the estimated parameters when applied to unseen samples. For the application of the sliding window technique, we set the length of the resulting time series p to be 206.  We recall that both steps of the approach are preceeded by a features extraction step, for which the latent space dimension k needs to be calibrated. As shown by the numerical results provided in Section C, the features extraction also guarantees the stationarity of the time series used for the anomaly detection task, ε.

Calibration of the Latent Space Dimension
When performing PCA, k is determined through a scree-plot, which is the representation of the proportion of variance explained by each component. The optimal k corresponds to the number of principal components explaining a given level of the variance of the original data. However, this method has its limitation as stated in (Linting et al., 2007). More importantly, it is not suitable for our approach. Indeed, in our case, the number of selected principal components must achieve a trade-off between information and noise in the latent space. If we consider a too low number of principal components, we may lose information regarding the normal observations, leading to false alarms. If a too high number of principal components is retained, we may include components representing noise, which prevents the model from detecting some anomalies.
To select the optimal dimension of the latent space, we consider the distribution of anomaly scores obtained through the naive approach of Section 2.4, empirically calibrating a cut-offŝ through the non parametric estimation of the distribution of the anomaly scores of uncontaminated and contaminated time series. We chose to calibrate the dimension of the latent space with the naive approach, because using the NN to this end would be very costly. Hence, for each k = 5, 10, . . . , 200, we construct a PCA model from which we infer reconstruction errors which are then converted into anomaly scores. Based on these anomaly scores, we tune the cut-off valueŝ thanks to the distributions and we finally convert the scores into labels. We evaluate the predictions of the naive approach for each value of k. The results on the evaluation metrics on the train set, as legitimate to choose k, are represented on Figure 4. The highest values are reached for k ∈ {40, . . . , 145}. The performance seems to be stable in terms of F 1 -score and accuracy. Therefore, for computational reasons, we choose the optimal k to be 40. Figure 4 also shows, as expected, a downward trend of the scores when represented against the highest values of k. This demonstrates that when a high level of variance is explained, it become much harder to perform anomaly detection based on the reconstruction errors.

Model Evaluation: Main Results
We evaluate the performance of the identification and localization stages of our approach on synthetic data, using appropriate performance metrics. We demonstrate numerically the efficiency of the PCA NN over baseline anomaly detection algorithms.

Contaminated Time Series Identification Step
For the features extraction step, we considered a latent space dimension k = 40. The NN built to compute the anomaly scores and convert them into labels was calibrated on the train set. The result of this calibration is shown in Table 3.   When an observation deviates significantly from the rest of the time series values, the model is able to recognise that the concerned time series contains an abnormal observation.

Data set Accuracy Precision
In Figure 6 with uncontaminated time series, we see that even when there is a local upward trend in the time series values, the model is able to make the distinction between this market move and the occurrence of an anomaly.
The last set of time series displays some limits of the model. The top graphs in Figure 7 represent the situation where the model did not manage to identify the contaminated time series. One explanation could be that the size of the anomaly is not significantly large enough to be spotted by the model. Indeed, looking at the graph, the time series does not seem to contain any anomaly. In contrast, the bottom graphs in Figure 7 show a stock price path predicted to be contaminated whereas it is not. We illustrate the robustness of the approach by assessing the model predictions on 100 distinct data sets.

Anomaly Localization Step
Regarding the localization step, the dummy approach consists in taking the argmax of time series observations as the anomaly location. We distinguish between two cases, anomalies which are extrema and anomalies which are not. The quality of the detection of our model is assessed for both types of anomalies. The results are displayed in Tables 5 and 6.   Numerical tests show the necessity of the features extraction step for the anomaly localization task, as applying the dummy approach alone is not enough when the anomaly is not an extreme value. Figure 8 represents a stock price time series with the true and predicted anomaly location.
In these examples, the locations are accurately predicted.
To guarantee the robustness of the model on the anomaly location, we evaluate the model prediction on 100 data sets. Table 7 shows the mean and standard deviation for the localization step of all types of anomaly.

Data set Accuracy Precision
Recall F 1 -score Train set 89.58% (4.3%) 89.58% (4.3%) 89.96% (4.1%) 89.67% (4.3%) Test set 89.49% (4.8%) 89.49% (4.8%) 90.00% (4.5%) 89.58% (4.7%) Table 7: Mean (standard deviation) of performance metrics of suggested model over multiple runs for localization step on all type of anomalies. Table 8 displays the results on the localization of non extrema anomalies. These results show the importance of the features extraction step. If we only consider the maximal observed value of the time series to be the anomaly, we would not be able to localize any non-extrema anomaly. Applying Algorithm 4 instead leads to satisfying anomaly detection rates.

Numerical Results Against Benchmark Models
To assess the performance, the suggested approach is compared numerically with well-known machine learning algorithms for anomaly detection, that is isolation forest (IF), local outlier factor (LOF), density based clustering of applications with noise (DBSCAN), k-nearest neighbors (KNN), and support vector machine (SVM), all reviewed in Section B. In addition to these state of the art techniques, we consider the recent anomaly detection technique proposed by Akyildirim et al. (2022), referred to as sig-IF, combining a features extraction step through path signatures computation with IF. Tables 9 and 11 summarize the performance on the train and test sets of Section 3 of each model, for both the contaminated time series identification and the anomaly localization steps.   The following paragraphs provide similar conclusions drawn from the results on the identification and localization steps (see Tables 9 and 11).
For the unsupervised learning methods, since there is not a proper learning step, even though Tables 9 and 11 show the scores on the train and test set, these scores should be seen as the ones obtained by testing the models on independent data sets. For IF, LOF and DBSCAN models, the performance across the sets is not stable, as a significant difference could be observed between the scores on the train and test sets. The poor performance of unsupervised learning algorithms could be explained by the difficulty to estimate the contamination rate for IF and the LOF algorithms. For DBSCAN the poor performance is rather due to the high dimensionality of the data. As for the sig-IF approach, its performance in this specific case is not as overwhelming as when used to detect pump and dump operations in (Akyildirim et al., 2022). This difference in performance is driven by the fact that, in our case, signatures are computed on the stock price path only, since no additional information is available to describe the price path we are analyzing. For pump/dump detection, instead, signatures are computed on a set of variables path including the price path and additional variables paths helpful in the identification of the pumps/dumps attempts.
Regarding the supervised approaches, although KNN is outperforming the suggested approach on the training set, there is a non-negligible decrease of these scores on the test set. This may suggest an over-fitting on the training data, which makes KNN unable to generalize what it has learnt to unseen samples. The same trend is seen on the scores with our approach, however the loss is no as harsh as in the KNN case.
Hence, according to the results, our PCA based approach seems to be the most suitable method for the problem at stake of anomaly detection on time series. Its satisfying performance in terms of accuracy and F 1 -score, as well as its low computational cost for both steps (see Tables  10 and 12), make the approach stand out.

Model Evaluation: Additional Results
To take our approach a leap beyond classical anomaly detection algorithms, we evaluate the anomalies imputation suggested by the PCA NN approach. We also test the robustness of the calibrated cut-off value and examine the sensitivity of the PCA NN performance to the amplitude of the anomaly.

Anomalies Imputation
To assess the imputation values suggested by our approach, we start from time series simulated without anomalies, then we randomly select a time stamp of the time series and add a noise to the corresponding value following the methodology described in Section 3. The imputation value the PCA NN suggests, is the reconstructed observation. This imputation technique is compared to naive methods of missing values imputation, namely backward fill (BF), consisting in replacing the anomaly by the previous value, and linear interpolation (LI). The choice of these imputation methods is motivated by their low computational cost. In fact, the PCA NN imputation comes at no additional cost as stated before. Therefore, we only challenge its performance using methods with similar complexity. To assess the quality of the imputation value suggested by each approach, we consider the following metrics. The imputation errors are computed on each time series i as n anom , 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6 6WRFN ,PSXWDWLRQ(UURU 3&$ $QR where S i refers to as the i-th path price with imputed values. The error on the covariance matrix is computed as where · Frob is the Frobenius norm, Σ the sample covariance matrix, and Σ the covariance matrix estimated on the data after the anomalies have been replaced by their respective imputation values. Each stock path was diffused and contaminated 100 times. The anomalies were then imputed following two baseline approaches (BF and LI). Figure 9 clearly shows that imputation using the reconstructed values with PCA reduces the imputation error but not as much as basic imputation techniques (LI and BF) as it is shown on Figure 10. Table 13 shows that the baseline imputation approaches achieve even lower errors on the estimation of the covariance matrix. Figure 11 shows a higher mean but a lower variance of the errors for the PCA-based imputation approach. In conclusion, it would be rather recommended to replace the flagged anomalies with approaches such as backward fill or linear interpolation.

Method
Ano PCA LI BF Clean Mean 0.004208 0.001575 0.001554 0.001554 0.001552 Standard deviation 0.000272 0.000165 0.000170 0.000172 0.000173 Table 13: Mean and standard deviation error on covariance matrix after imputation of anomalies with baseline methods.
A natural explanation to the relative inefficiency of the reconstruction value as imputation value comes from the fact that, by construction, the reconstructed value integrates in its computation the abnormal value, which is not the case when using other imputation techniques. Therefore, even if the reconstructed value is closer to the true value than it is to the abnormal value, the spread between the imputation and the true value is still significant. Figures 12 and  13 illustrate the latter with a plot of a reconstructed and original price path. The two additional graphs represent the region around anomalies, where the cross and circle respectively show the true value of the stock and its value after imputation of the anomalies using the reconstructed value suggested by the PCA NN approach.

Cut-off Value Robustness
In most anomaly detection models, the cut-off value is a hand-set parameter. In our approach, the cut-off value is a model parameter and as such it is calibrated through the learning. Testing the robustness of the cut-off given by the model is therefore a must. The identification step is the unique step of the approach being concerned by the robustness, since it is the only step involving the cut-off calibration. Robustness is checked by shocking the suggested cut-off with different level of noise and observing the impact of these shocks on the model performance. We consider several shock amplitude γ ∈ ±{10 −4 , 10 −3 , 10 −2 , 10 −1 , 1, 2}.  Excluding the extreme cases where the shock amplitude is ±{1, 2}, one could see that the accuracy is barely impacted by the shocked cut-off values, both on the train and test sets. The interpretation of the remaining results is split in two and is applicable for both the training and test sets.
When negative shocks are applied, the model predicts more anomalies and less normal observations (compared to the predicted numbers with the calibrated cut-off value). Therefore, the model is able to identify anomalies which were missed initially. Hence, applying negative shocks increases the recall. As for the precision, i.e. the rate of correctly identified contaminated time series, since on left hand side of the cut-off value we are in the density region where contaminated time series represent the dominant class, the new position of the cut-off value leads to a misidentification of this type of observations. This entails a deterioration of the precision. Opposite behaviours of the precision and recall are observed when positive shocks are applied, since some abnormal time series are missed by the model. Figure 14 shows the precision and recall when cut-off values other than the one we calibrated are used in our model, on the train and test sets. These scores are also compared with the performance of the no-skill model 1 . We note that our approach edges out the no-skill model on both sets. The area under the curve (AUC) score 2 on the train set and the test set are respectively 0.97 and 0.87. This shows that we are better performing on the train set, which was expected, but the performance on test set is just as satisfying. These results reinforce our conclusions regarding the robustness of our approach. One can see that the calibrated cut-off value represents, for the train set, the point where an equilibrium is being reached between precision and recall. The calibrated cut-off allows to achieve high precision and recall scores, simultaneously. We conclude that the cut-off given by the model is suitable for the training samples and for unseen samples as well.
To briefly sum up this section, although we mentioned that some shocks on the cut-off value induce higher scores, the improvement over the scores with the calibrated cut-off value is not significant (unless high amplitude shocks are considered). The numerical tests and the precision-recall curves are consistent with the robustness of the cut-off value suggested by the approach.

Sensitivity to Anomaly Amplitude
Time series are manually contaminated as described in Section 3.2. The abnormal value S a,i t of the i-th time series results from a shock of the initial value of the time series denoted by S i t as stated in (9). We recall that the shock is represented by δ, whereas its amplitude |δ| is uniformly drawn from [0, ρ]. Hence, ρ is the parameter that ultimately controls the amplitude of the anomaly. Since ρ is fixed by the user, it is interesting to investigate the sensitivity of the PCA NN approach performance to the amplitude of the shock.
The PCA NN approach evaluated in Sections 5.1-5.2 is calibrated on time series which were contaminated with shock amplitude drawn from U ([0, ρ]) with ρ = 0.04. For the data sets on which the model was calibrated and then evaluated, the anomalies were grouped according to the shock amplitude they result from. For the identification step, we distinguish four groups of contaminated times. For instance, the first line of Table 15 defines the first group of contaminated time series, for which the anomalies results from shock amplitude ∈ [0.3091e − 2, 1.46e − 2[. For this first group with this specific range of amplitude, 77% of the contaminated time series were identified during the PCA NN identification step.
Similarly, for the localization step, we consider four groups of anomalies. The first line of Table 16 represents the anomalies with shock amplitude ∈ [0.309e − 2, 1.31e − 2[. 86% of the anomalies belonging to this first group were correctly localized by the PCA NN localization step.  The bounds defining each group are chosen to be the minimal value, 25%-quantile, 50%quantile, 75%-quantile, and the maximum value over the shock amplitude. From the results of Tables 15 and 16, one could see that, except for the first group which represents the lowest shock amplitude, the detection ratio of the remaining groups are similar.
The calibrated PCA NN approach performance is then tested on new data sets contaminated with distinct values of ρ. Figures 15 and 16 show that the F 1 scores are almost similar across the high values of ρ, both for the train and test sets. Moreover, Figure 16 clearly states similar and high performance of the localization stage on both data sets, regardless of the value of ρ. To conclude on the sensitivity of the model performance to the amplitude of the anomaly shocks, although it is clear that a lower performance is observed on the identification and localization of smaller amplitude shock anomalies, the detection ratio is still satisfactory for this category of anomalies.

Application to a Downstream Task: Value-at-Risk Computations
In this section, we illustrate the benefit drawn when applying the PCA NN approach as a pre-processing before value-at-risk computations. Given a random variable representing the loss in portfolio position over a time horizon h, its value-at-risk at the confidence level α ∈ ( 1 2 , 1), VaR α ( ), is defined by the quantile of level α of the loss distribution, i.e. P( ≤ VaR α ( )) = α (assuming atomless for simplicity). Let (S t ) t=t 1 ,...,t T be a path price diffusion distributed according to the Black-Scholes model (7). The logarithmic returns to maturity are distributed according to the Gaussian distribution with µ ∈ R and σ > 0. We assume the vector R of log-returns on our stocks to be joint-normal, i,R and Σ R is the covariance matrix. We consider a portfolio on N stocks, which return is given by P = Q R, where Q ∈ R N defines the composition of the portfolio. Hence, P ∼ N µ P , σ 2 P , where µ P = Q µ R and σ 2 P = Q Σ R Q. The value-at-risk for the time horizon h at level α is where q α is the α-quantile of a standard normal distribution and the parameters µ R and Σ R are estimated from different types of time series to evaluate the impact of localizing and removing anomalies following our approach. Under the adopted framework, the true VaR α (P ), VaR theo is thus known and can be computed using the diffusion parameters. An estimation VaR α (P ) can be obtained by replacing the parameters in (15) by their estimates µ P and σ P computed from the time series with anomalies, or after imputation of anomalies. Then, absolute errors and relative errors are computed as AbsoluteError VaR = VaR α (P ) − VaR α (P ) ,  To conduct this analysis, we generate new stock path samples that we assume to be clean, fixing the diffusion parameters µ R and Σ R . We then add the anomalies following the procedure described in Section 3. We apply our model to localize the anomalies and replace the localized anomalies using the backward fill (BF) approach (shown to be the more efficient imputation technique in Section 6.1). For each stock and each run of simulation, we obtain four estimates of the distribution parameters of the associated log-returns.   Figure 17 and Table 18 summarize the distribution of the VaR estimates for α = 0.99 and h = 1(day), over several simulation runs. The boxplots show the dispersion of the portfolio VaR estimates on several diffusions. The green square represents the mean of the VaR estimates. For the four first boxplots, the means are approximately on the same level, which is confirmed by the results of Table 18. The anomalies present among the time series observed values have a non-negligible impact on the distribution parameters estimation, which ultimately causes a wrong estimation of the VaR. Thanks to the localization of the anomalies by the suggested model and their imputation as per Section 6.1, we are able to get a more accurate estimation of the VaR. VaR loc,true and VaR loc,pred are quite similar, which shows that the model accurately localizes the anomalies.

VaR
VaR clean VaR loc,true VaR loc,pred VaR anom Absolute Error 0.007995 0.008596 0.008622 0.02235 Relative Error 0.014620 0.015720 0.015767 0.04087 Table 19: Mean absolute and relative error on VaR estimations for P .
We also evaluate the error on the VaR estimation using the mean absolute error and the mean relative error, taking the VaR theo as our benchmark. As one can tell from Table 19, even when the distribution parameters are estimated from the clean time series, the VaR computed with these parameters is not exactly the one computed with the theoretical parameters. This can be explained by the historical size of the observed values used to estimate the parameters. This table shows that by removing anomalies we can reduce by a factor two the error on VaR estimation.
Additionally, we assess the impact on VaR theo , VaR clean , VaR anom , VaR loc,true and VaR loc,pred of increasing n anom . To this end, we perform 50 simulation runs of stock paths for n anom and for each of those scenarios we estimate the VaR on the portfolio. We summarize the results on Figure  18, where each curve represents the mean VaR estimation with respect to n anom along with a representation of the uncertainty around each evaluated point through a confidence interval.
When we compute the VaR using the abnormal time series, we notice that the difference between VaR anom and VaR theo increases with n anom , which is natural to expect. However, when the time series are cleaned prior to VaR estimation, the curve representing the VaR estimates are much closer to the ones representing VaR theo and VaR clean , showing an undeniable improvement in the accuracy of the VaR estimation over the estimation based on abnormal time series. Furthermore, VaR estimation after the imputation following the model prediction or knowing the true localization of the anomalies seem to be quite similar for low n anom , while some discrepancies between the two become more significant as n anom increases. A natural explanation could be that when the number of anomalies increases and the model suggests wrong anomalies localization, normal values are being replaced while true anomalies remain among the observed values, which wrongly impacts the VaR estimation. However, the results of Tables 20 and 21 indicate that the anomaly localizations suggested by the model are overall correct and allow removing most of the anomalies, as the relative error of VaR loc,pred is, regardless of n anom , always lower than the relative error of VaR anom (e.g a relative error of 0.0248 for VaR loc,pred , against 0.0658 for VaR anom , when there are 76 anomalies among the 1, 500 observed values of the time series).

Numerical Results on Real Data
We consider a labelled real data set including stock prices, bonds yields, CDS spreads, FX rates, and volatilities. These data were collected from a financial data provider for the period between 2018 and 2020. These data sets were labelled by experts. They are provided on https: //github.com/MadharNisrine/PCANN in the form of 132,000 time series (after augmentation). The train set is balanced, while the test set is imbalanced with 20% of contaminated time series.
To ensure the independence of the train and test data sets, we calibrate the PCA NN approach considering the time series of 2018 and 2019, while the 2020 time series are dedicated to model evaluation. The performance evaluation of the PCA NN approach after its calibration on real data sets is shown in Table 22. As visible from the upper panel, the PCA NN performs quite well for the identification step. Once the contaminated time series are identified, the model is able to localize the anomaly with high accuracy, as reflected by the scores on the test set in the lower panel.

PCA NN against State of the Art Models on Real Data
Even if the F 1 -score on the identification step is not that high, the results show that the related performance remains better than the one of alternative state of the art approaches. Tables 23 and 24 compare the performance of the PCA NN approach against the state of the art models of Section B on our real data set. PCA NN outperforms the benchmark models in both steps, with overwhelming results for the anomaly localization step.

Model
Accuracy

Conclusion
We propose a two step approach for detecting anomalies on a panel of time series that can reflect a wide variety of market risk factors. The first step aims at identifying the contaminated time series, i.e. time series with anomalies. The second step focuses on the localization of the anomaly among the observed values of the identified contaminated time series. As preprocessing, our methodology integrates the extraction of features from the time series with PCA. This part of the method proves to be essential, as it provides the models with inputs on which the distinction between abnormal/contaminated and normal instances is eased, while also ensuring the stationarity of the model (time series) inputs. Another key point of the approach is the calibration of the cut-off value, the key parameter in the identification of contaminated time series, by means of a feedforward neural network with a customized loss function. The proposed approach suggests an imputation value, however this value is strongly influenced by the abnormal value. Therefore, basic imputation approaches with similar complexity are preferred. Our numerical experiments show not only that our approach outperforms baseline anomaly detection models, but also show the real benefit that could be gained from applying it as a data cleaning step preliminary to VaR computations. Future research could focus on the replacement of PCA by partial least squares (PLS) or deep PLS (Polson et al., 2021) for endogeniging the features extraction stage.
Regarding downstream tasks, our approach might be of special interest for reverse stress tests (Eichhorn et al., 2021).

A Literature Review
We start with a review of the anomaly detection literature. See also Section B for a more technical presentation of some of the below-mentioned algorithms.

A.1 Baseline Algorithms
Anomaly detection aims at finding an "observation that deviates so much from other observations as to arouse suspicion that it was generated by a different mechanism" (Hawkins, 1980). The baseline anomaly detection algorithms, described in (Chandola, 2009), struggle to identify anomalies in time series, mainly because their assumptions are invalidated. If we consider models built for spatial data, a major assumption of these models is that observations are independent, whereas, for time series, high dependency exists between different time stamps.
Clustering-based approaches, like the density based spatial clustering with noise (DBSCAN) method, are particularly impacted by this aspect: if an anomaly occurs at a given time stamp and is followed by incorrect values, clustering-based approaches consider that the observations of the time series belong to two different clusters and thus fail to identify the anomaly. Another limitation when considering this type of techniques is the choice of the similarity metric used for data clustering. This task, although being a crucial pillar of these approaches, is not trivial and becomes very challenging for high dimensional problems.

A.2 Statistical Approach to Anomaly Detection
The statistical techniques for anomaly detection can be split in two families: statistical tests and predictive models. Both suffer from the curse of dimensionality and model/data mismatch. When anomaly detection relies on hypothesis tests, it usually tests whether the observations are drawn from a known distribution (Zhang and Paschalidis, 2017), supposing that the user knows the probability distribution of the normal observations. Such a parametric framework narrows down the scope of applicability of hypothesis tests, as the data does not always coincide with the assumed distribution. Moreover, the tests provided in the literature are not suitable in multivariate settings (Kurt et al., 2020). The statistical techniques relying on fitting a predictive model to each time series also require strong assumptions on the data. Predictive models are usually autoregressive (AR), moving average (MA), or ARMA models. Anomalies are then detected relatively to the forecasts suggested by the model (Chandola, 2009). In these parametric approaches, some parameters have to be specified again, starting with the order of the models. Selecting the optimal model parameters with respect to an information criterion is not always possible. Additionally, these models assume that the time series are homogeneous, i.e. drawn from the same distribution . This is not always satisfied in the financial risk management case where several types of market risk factors are treated simultaneously.

A.3 Score Based Anomaly Detection Models
Additional anomaly detection challenges are of general concern. Most of anomaly detection algorithms are score-based, in the sense that these approaches return an anomaly score reflecting to which extent the observation is considered to be abnormal by the model. In order to decide whether an observation is abnormal or not, a cut-off value of the score has to be selected. Empirical approaches are often used, consisting in setting the cut-off value as a quantile or elbow point of the distribution of the anomaly score. However, the selected cut-off value according to such methods remains arbitrary. (Gao and Tan, 2006) propose to rely on a cost of misclassification based on a weighted classification accuracy. Approaches to calculate "optimal" weights are described in (Lu et al., 2019), but they involve a heuristic grid search technique. Another alternative is to determine the cut-off value by cross validation on the training data (Saha et al., 2009). Finally, some methods do not select any cut-off value, but are based instead on a contamination rate. However, fixing a cut-off value or deciding on a contamination rate is not so different.

A.4 Scarcity of Anomalies and Data Augmentation
The scarcity of anomalies within the data sets is another typical problem in anomaly detection. Anomalies are, by definition, rare events, therefore they are under-represented in the data set used to fit the models. This under-representation is not helping in the design of a reliable model able to identify anomalies. Classical methods to overcome this issue consider data augmentation. These techniques aim at producing new synthetic samples that will ultimately enhance model performance, leveraging on a better representation of the features space. As reported in , time series can be augmented by using a simple transformation in time domain (Cui et al., 2016), frequency domain , or more advanced generative approaches involving deep learning techniques, such as recurrent generative adversarial networks (Esteban et al., 2017). In practice, the use of generative models for data augmentation of time series with anomalies presents two limitations. First, training such models requires a large number of samples to guarantee a satisfactory performance. While restricted Boltzman machines do not exhibit this problem, (Kondratyev et al., 2020) have shown that they fail into fitting multivariate complex distributions with nonlinear dependence structure. A more fundamental limitation affects the very idea behind generative models. Such networks are trained to learn a given distribution. However, by definition, anomalies are different from each other and therefore there is not a distribution that characterises them.

A.5 Supervised vs. Unsupervised learning
Since anomalies are the realisations of atypical events for which the distribution is unknown, it seems quite natural to use unsupervised algorithms. However, these approaches are deemed more suitable for learning complex patterns and are task specific. Moreover, (Görnitz et al., 2013) pointed out that they often do not present a high prediction performance, in particular in high dimensional settings (Ruff et al., 2019). Indeed, the performance of unsupervised shallow anomaly detection algorithms depends upon a features engineering step. (Akyildirim et al., 2022) use signatures to extract features which then feed algorithms such as isolation forest. This combination of techniques is shown to over-perform benchmark approaches. However, the designed model is task specific (detection of pump and dumps attacks) and the features extraction step is only efficient when at least one explanatory variable is considered in the analysis. Moreover, the unavailability of labelled data makes the model building and its evaluation even more complex. As for supervised methods, the only limitation on which the literature tends to agree is their incapacity to generalize the learned patterns to new samples, which is the consequence of misrepresentation of anomalies among the training samples (Zhao and Hryniewicki, 2018). However, the scarcity of labelled data can be sidestepped through the use of data augmentation techniques on the fraction of available labelled data. For these reasons the supervised learning framework is to be preferred even when only a small set of labelled data is available.

A.6 Anomaly Detection on Time Series
Usually, anomaly detection models on time series have two main components. The first component aims at extracting a parsimonious yet expressive representation of the time series. Several approaches are suggested in the literature to deal with such features extraction. Recently, deep neural networks have been shown to suffer from overparametrization and to be often computationally expensive (Dempster et al., 2020;Akyildirim et al., 2022). Path signatures are also computationally demanding. This could perhaps be alleviated by the random signatures (Compagnoni et al., 2022). However, the information extracted with signatures is of most interest when the considered paths are characterized by several variables. The resulting representation is then transformed into a (typically continuous) anomaly score which is in turn converted into a binary label (Braei and Wagner, 2020).

A.7 PCA and Anomaly Detection
In the literature, PCA is usually used in anomaly detection for its dimension reduction properties. Anomalies are identified on the latent space, either by applying some anomaly detection algorithm or by assuming a given distribution on the principal component and identifying the anomalies relatively to a quantile (Shyu et al., 2006). Assuming that a normal subspace representation of the data set can be constructed with the first k-components (Ringberg et al., 2007), anomaly detection is achieved by looking at the observations that cannot be expressed in terms of the first k-components (Bin et al., 2016). Note that the particular power-fullness of auto-encoder, a non-linear PCA, as data compressor is not desirable herein, since auto-encoders compress all patterns including abnormal ones. While (Ding and Tian, 2016) claimed that the PCA-based models are stable with respect to their parameters, such as the number of principal components k spanning the subspace or the cut-off level, (Ringberg et al., 2007) found instead that PCA-based anomaly detection is sensitive to these parameters and to the amplitude of the anomalies. Indeed, the latter may undermine the construction of the normal subspace representation, in turn leading to misidentification of anomalies. In light of that, we take an extra care regarding these aspects, and appropriate tests are conducted to show that the proposed approach is not subject to these issues. In the end, with a relatively low number of features given by PCA, we are able to accurately describe the dynamics of market risk factors represented by times series. The reason standing behind that is the high correlation structure displayed by market risk factors.

B State of the Art Anomaly Detection Models
There are mainly two categories of machine learning anomaly detection models : density-based models, where a distribution is used to fit the data and anomalies are defined relatively to this distribution, and depth-based models, which, instead of modeling the normal behaviour, isolate anomalies.

B.1 Density-Based Models
Density based spatial clustering with noise (DBSCAN) DBSCAN (Ester et al., 1996;Schubert et al., 2017) is an unsupervised clustering methodology that groups together comparable observations based on a similarity metric. Clusters are high density regions and are defined by the ε-neighbourhood of observations and by MinPts, the minimum number of points required to be in a radius of ε from an observation to form a dense region. An anomaly is any observation which has not M inP ts in its ε-neighbourhood nor appears in the ε-neighbourhood of other observations.
K-nearest neighbours (KNN) Usually used for classification purposes, KNN (Hand, 2007) is a supervised algorithm that can be used for anomaly detection as well. For each observation it selects its closest K observations generally in terms of distance, but other similarity metrics can be considered. The anomaly score of an observation is computed as function of its distances to K-nearest neighbours, weighted average of the distances for instance. The points with the highest anomaly scores are considered as anomalies.

Support Vector Machines (SVM)
The ultimate aim of SVM (Cortes and Vapnik, 1995) is to define a hyperplan that separates the data. The specificity of this hyperplan is that it maximizes the distances to the set of features representing each class. When the data is not linearly separable, a map φ is applied to the initial features vector so the data became linearly separable in the new space in which it was projected.

B.2 Depth-Based Models
Isolation Forest (IF) The IF method (Hariri et al., 2019) applies a depth approach to detect anomalies. The algorithm is based on the idea that anomalies are easier to isolate and thus will be isolated closer to the root, while normal observations are isolated much further from the root. The algorithm uses random decision trees to separate the observations. The anomaly score is calculated as the path length to isolate the observation. As it is defined in the algorithm, the anomaly score will be closer to 1 for anomalies and 1 for normal observations. This allows ranking the observations, from the most abnormal observation to the most regular one. However, the choice of an explicit cut-off between this two type of instances is not obvious.
Local Outlier Factor (LOF) The LOF algorithm (Breunig et al., 2000;Alghushairy et al., 2020) tries to assess the isolation of one observation relatively to the rest of the data set, before flagging it as an anomaly. This model relies on the concept of local density. The anomaly score of an observation x is its local outlier factor, which quantifies how dense is the location area of x compared to the one of its neighbours. Hence, for each observation, LOF k (x) ≈ 1 means that the density of observations around x is similar to the one of its neighbours, therefore x could not be considered as an isolated observation. LOF k (x) 1, instead, shows that the density of x is lower compared to its neighbours, hence x should be flagged as anomaly

C The Data Stationary Issue
Assume that a process ε i t satisfies the following representation ∆ε i t = γε i t−1 + θ 1 ∆ε i t−1 + . . . + θ p−1 ∆ε i t−p+1 + z t , where p is the lag order, ∆ is the difference operator, i.e ∆ε t = ε t − ε t−1 , and z t is a white noise. The stationarity of (ε i t ) is shown using the augmented Dickey-Fuller test (Fuller, 2009). Namely, the process ε i is stationary if there is a unit root, i.e. γ = 0. Therefore, for each ε i the test is carried under the null hypothesis H 0 : γ = 0 against H 1 : γ < 0.  If the p-values are lower than the significance level, then the null hypothesis is rejected for all the reconstruction errors and we conclude that the time series we work with do not suffer from the non-stationarity issue.