{'fraction_non_alphanumeric': 0.03628052547049376, 'fraction_numerical': 0.01639006712514928, 'mean_word_length': 4.4829532625874915, 'pattern_counts': {'":': 0, '<': 3, '<?xml version=': 0, '>': 0, 'https://': 0, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 4, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Modeling Multivariate Cyber Risks: Deep Learning Dating Extreme Value Theory

Modeling cyber risks has been an important but challenging task in the domain of cyber security. It is mainly because of the high dimensionality and heavy tails of risk patterns. Those obstacles have hindered the development of statistical modeling of the multivariate cyber risks. In this work, we propose a novel approach for modeling the multivariate cyber risks which relies on the deep learning and extreme value theory. The proposed model not only enjoys the high accurate point predictions via deep learning but also can provide the satisfactory high quantile prediction via extreme value theory. The simulation study shows that the proposed model can model the multivariate cyber risks very well and provide satisfactory prediction performances. The empirical evidence based on real honeypot attack data also shows that the proposed model has very satisfactory prediction performances.

things worse is that cyber systems have become underlying pillars for various infrastructures and physical systems. For example, there are growing concerns about cyber threats to critical infrastructures because of the linkage between cyber and physical systems [21]. In addition, cyber attacks causes billions of dollars in losses each year based on the research by the Ponemon Institute and published by IBM Security [11]. In light of the risk and potential consequences of cyber risks, securing cyberspace has become not only a mission of homeland security but also an emerging task for private companies. This clearly calls for research into cyber risk modeling, of which risk prediction is an important task.
In this work, we investigate a particular challenge that is encountered when modeling and analyzing cyber risk, namely the modeling of multivariate cyber risks. This challenge is imposed by the high-dimensional dependence among cyber attacks and extreme attacks during a short time period. The dependence among cyber attacks is very common which is the nature of cyber risks. For example, the Distributed Denial of Service (DDoS) attack, a large-scale DoS attack where the perpetrator often uses thousands of hosts infected with malware to launch attacks to a target or multi targets (e.g., ports, computers, severs), is in nature high-dimensional and dependent. In particular, the DDos attacks are often extreme during a short time period. The defender often needs to estimate or predict the cyber risks for the purpose of adjusting the defense posture in practice. However, this is very challenging because of the high-dimensionality and heavy tails (i.e., extreme attacks during a short period) exhibited by cyber attack data [28,29,27,22], which hinders the development of modeling multivariate cyber risks. Therefore, in the literature of statistical modelings, there are only few studies investigating multivariate cyber risks. For example, Xu et al. [27] proposed a vine copula approach for modeling the dependence among the time series of the number of cyber attacks, and the dependence between the time series of the number of attacked computers. Peng et al. [22] developed a copula-GARCH model, which also uses vine copulas to model the multivariate dependence among cyber attacks. Ling et al. [18] proposed a vector autoregression (VAR) approach to identifying a geospatial and temporal patterns in the cyberattacks by considering the long range dependence.
In recent years, there is a growing interest in developing deep learning models for time series forecasting [17,5,23]. For multivariate time series forecasting, Goel et al. [7] proposed a hybrid R2N2 (residual recurrent neural networks) model, which first models the time series with a simple linear model (e.g., VAR) and then models its residual errors using recurrent neural networks (RNNs). Lai et al. [16] developed a long-and short-term time-series network (LSTNet) to model the multivariate time series. The main idea is to use the convolution neural network and the recurrent neural network (RNN) to model short-term local dependence and long-term trends for time series. Che et al. [1] studied a deep learning model based on gated recurrent unit (GRU) for modeling multivariate time series with missing data. Wang et al. [25] presented an end-toend framework called deep prediction interval and point estimation, which can simultaneously perform point estimation and uncertainty quantification for multivariate time series forecasting. The main approach is to model the loss function by penalizing the loss of point estimation and the loss of prediction interval based on the deep learning model of long short-term memory network. To the best of our knowledge, there is no deep learning framework for modeling the multivariate cyber attack data in the domain of cyber security.
In this work, we propose a novel framework which is different from those in the literature for modeling and predicting the multivariate cyber attack data. Specifically, we develop a hybrid model by combining deep learning and extreme value theory [4,20]. That is, we first model the multivariate attack time series via a deep learning model which aims to capture the highdimensional dependence via the deep learning network. We then model the residuals exhibiting heavy tail via the extreme value theory. The proposed approach not only can provide accurate point prediction but also capture the tails very well (i.e., predicting the extreme attacks). This novel framework is particularly useful in practice as it can guide the defender to prepare the resource for both the regular attack (in terms of mean prediction) and worst attack (in terms of tail prediction) scenarios.
The rest of the paper is organized as follows. In Section 2, we present the preliminaries on the deep learning model and extreme value theory. In Section 3, we discuss the proposed the framework for model fitting and prediction. Section 4 assesses the proposed approach via simulation studies. In Section 5, we study two real honeypot attack data. In the final section, we conclude our results and present a discussion towards some limitations and future work.

Long short-term memory (LSTM)
The LSTM network, a special kind of RNNs, introduced in [9] has been received remarkable attention because of the ability to learn long term patterns in sequential data and tremendously improve the prediction accuracy compared to other deep learning models. Therefore, this deep learning model has been successfully applied in many areas [8].  Figure 1 illustrates a LSTM block. Compared to regular RNN, the LSTM employs a different approach for the activation. Specifically, given time series data y 1 , . . . , y t where y i = (y 1,i , . . . , y n,i ) , i = 1, . . . , t, the activation h t of LSTM at step t is computed based on four pieces: gate input (a.k.a, information gate), the forget gate, the output gate, and the cell gate. The information gate input at step t is where σ is the sigmoid function, and W's are the parameter matrices of the network learned during the training. The forget gate input and the output gate input are computed as The cell gate input is computed as where represents the element-wise Hadamard product, the C t−1 is the cell state information from the previous step, and k t is where tanh is the hyperbolic tangent function. Finally, the activation at step t is computed as where h(t) ∈ R n is the final cell output.

GRU and mLSTM
As the variants of LSTM, GRU [2] and Mulitplicative LSTM (mLSTM) [15] have also received much attention recently due to their extraordinary prediction abilities. The GRU uses less trainable parameters than LSTM, by having only two gates, the reset gate (r t ) and the update gate (z t ). Its computational process is as follows Compared to the LSTM, the mLSTM simply replaces h t−1 with an intermediate state, m t , in order to calculate i t , f t , o t , and k t , where m t is For more and detailed discussions on the deep learning models, please refer to [9,8] 3 Deep learning framework for multivariate cyber risks Let y t = (y 1,t , . . . , y n,t ) be a vector of attack time series. In practice, the dimension of n can be very large, and particularly, the dependence among cyber risks can be very nonlinear. This motivates us to use the following approach for analyzing multivariate cyber risks, which utilizes the advantage of high accurate point estimates of deep learning and the advantage of EVT for modeling the high quantiles. Specifically, the proposed model is as follows, where t = ( 1,t , · · · , n,t ) is an unobservable zero mean white noise vector process, and f is a function by mapping the vector (y t−1 , . . . , y t−p ) to the mean vector of y t . Note that a special case is the VAR, where A i are coefficient matrices, i = 1, . . . , p, and µ 0 is an n-dimensional constant vector. We propose to model the mean of y t byŷ Stage 1: Capturing the multivariate dependence via deep learning In multivariate time series modeling and prediction, the most challenging part is to find a suitable f for capturing the complex dependence pattern. Deep learning as an emerging tool for multivariate time series modeling and forecasting has an extraordinary forecasting performance [6,7]. Therefore, we develop the following procedure for modeling the mean function f via deep leaning. Given time series y 1 , . . . , y t−1 , we define p as the number of input data points in each of them, and p also indicates the number of time steps to look back (i.e., lags). The prediction process can be formally described as:ŷ where Y t−1,p = (y t−1 , . . . , y t−p ) ∈ R n×p is the input matrix. We apply the following objective function to measure the error generated by the model, where λ is the tuning parameter, || · || 2 2 represents the squared L 2 norm, and W represents all the training parameters in the learning process. For instance, if the network is LSTM, then The optimization is defined as which can be solved by using the gradient descent method [14]. To select the best deep learning model, we perform a manual grid search by setting the hyper-parameters of the networks as follows: • Network models: where e t = (e 1,t , . . . , e n,t ) . The second stage is to model the residuals by some statistical distribution. Since the heavy tails are often observed in the attack data, we propose to model the high quantiles of the residuals by the EVT approach. This is in principle in line with the two-stage pseudo-maximum-likelihood approach in [19].
Recall that a popular EVT method is known as the peaks over threshold approach [20,4]. Specifically, given a sequence of i.i.d. observations X 1 , . . . , X n , the excesses X i − µ with respect to some suitably large threshold µ can be modeled by, under certain mild conditions, the Generalized Pareto Distribution (GPD). The survival function of the GPD is: , and ξ and σ are respectively called the shape and scale parameters. Split the data set into mini-batch of size r; 3: while j <= b do 10: for each batch from the training data do Compute J by performing forward propagation using λ;

12:
Update W using gradient descent based on J with the learning rate γ; OUTPUT: Deep learning model M * with the corresponding parameters W * , and fitted valueŝ y t , t = 1, . . . , T .

Stage 3: Prediction and and evaluations
After we determine the deep learning model and residual distribution, then we can use the forward propagation to predict the meanŷ The q-quantile of y i,t+1 can be predicted aŝ   Compute residuals e t = y t −ŷ t , t = 1, . . . , j;
The prediction performance for the mean prediction are evaluated based on MSE and mean absolute percentage error (MAPE) [10].
In order to assess the prediction accuracy of quantiles, we propose using the Value-at-Risk (VaR) [20] metric because it is directly related to the quantities of interest. Recall that for a random variable X t , the VaR at level α for some 0 < α < 1 is defined as VaR α (t) = inf {l : P (X t ≤ l) ≥ α}. An observed value that is greater than the predicted VaR α (t) is called a violation, indicating inaccurate prediction. In order to evaluate the prediction accuracy of the VaR values, we adopt the following two popular tests [3]: (i) the unconditional coverage test, denoted by LR uc , which evaluates whether or not the fraction of violations is significantly different from the model's violations; (ii) the conditional coverage test, denoted by LR cc , which is a joint likelihood ratio test for the independence of violations and unconditional coverage.

Simulation study
In this section, we perform a simulation study on assessing the performance of proposed approach. To mimicking the attack data with heavy tails, we randomly generate two data sets with size 5000 × 5 from the following models. In each experiment, the simulated data is split into three parts: the first 3500 is used for training, and the following 500 observations is used for the validation. The testing data is the last 1000 observations.
• VAR with heavy tail. The VAR model is set as follows.
The prediction is performed via Algorithm 2, where the threshold is set to .9-quantile of residuals. The prediction performance is reported in Table 1. For the comparison purpose, we also report the prediction performance of VAR, where the lag p is selected via AIC criterion. The VAR model is considered as a benchmark model in the sequel discussion. It is seen from For the quantiles prediction, we plot the violations for the proposed model and benchmark model in Figure 3   .99 level, it is seen that the proposed approach underestimates the number of violations for time series 4 and 5, which result in relatively smaller p-values of LRuc and LRcc tests in Figures 3(b) and 3(c). For the benchmark model, we observe that it underestimates four time series at level .99 from Figure 3(d).
We conclude that the proposed approach has an overall satisfactory prediction performance for high quantiles from level .95 to .98.
To simulate the dependence among Z i,t , we generate 5-dimension multivariate uniform distribution, where the dependence structure is specified via the R-vine copula. Specifically, the R-vine tree matrix is as follows A 5000×5 matrix from the above R-vine structure is generated, which is used for generating the dependent Z i,t , i = 1, . . . , 5. The simulated time series plots are displayed in Figure 4, and the extreme values can also be observed.
We again employ Algorithm 1 to train and validate the deep learning model, and the residuals are used to fit the GPD distribution in Eq. (3.2). It is discovered that the threshold of 90th percentile of residuals can produce satisfactory fitting performance for all the time series. Therefore, in the prediction procedure of Algorithm 2, the threshold is fixed at 90th percentile of residuals for each time series. Table 2 presents the performances For the high quantiles, the number of violations are displayed in Figures 5(a) and 5(d). It is seen that the proposed model has a very satisfactory prediction performance compared to that of benchmark model. The large p-values of LRuc and LRcc tests in Figures 5(b) and 5(c) compared to those in 5(e) and 5(f) also confirm that the proposed approach has an accurate prediction performance for the high quantiles.
To conclude, the proposed model has satisfactory prediction performances for both the point prediction and high quantile prediction. It can significantly outperform the benchmark model.

Empirical study
In this section, we study two real attack data which are collected by the honeypot instrument [24].

Honeypot data I with 9 dimensions
This honeypot data is publicly available on the web [12], which was collected via Amazon Web Service (AWS) virtual honeypots across the world. The dataset has 9 honeypot hosts (EU, Oregon, Singapore, SA, Tokyo, Norcal1, Norcal2, US-East, Sydney), and the attacks were recorded between 03/03/2013 to 09/08/2013, which includes 451,581 events. The recorded attack data includes attack time, targeted host, attackers' IP addresses and origin countries. This dataset was discussed in [18] based on the daily aggregation, where the long range dependence model was discovered and incorporated into the modeling process.
Daily aggregation. The daily aggregated data has 188 observations, we leave the last 38 observations as the prediction evaluation as that in [18] where the predicted MSE is reported. We employ Algorithm 1 to build the deep learning model where the training data has 100 observations, and validation data has 50 observations. The prediction performance based on Algorithm 2 is reported in Table 3. The proposed deep learning has a satisfactory prediction performance. In particular, it outperforms the approach in [18] for all the attack time series except for Singapore based on MSEs (0.0145 vs 0.0133). Since the daily aggregated time series only has a total of 188 observations for each host, the high quantile predictions are not performed.
Hourly aggregation. In practice, the network defender is often interested in assessing the hourly attacks [28,29,22]. In this section, we study the performance of proposed model based on the hourly aggregated data which has 4512 observations for each host. The hourly aggregated time series are plotted in Figure 6. It is observed that there are extreme attacks for all the hosts. Particularly, there exist tremendous numbers of attacks during some short periods for Oregon, Singapore, and Tokyo hosts. The data is split into three parts: training with 3500 observations, validation with 500 observations, and the last 512 observations are used for the prediction evaluation. The Algorithm 1 is employed on the training and validation datasets to select the best deep learning model, and Algorithm 2 is used for the prediction.  The prediction performance is reported in Table 4. For the comparison purpose, the prediction performance of benchmark VAR model is also reported. For the proposed approach, we can observe that it can significantly outperform the benchmark model in terms of MAPE metric. For the MSE metric, it is seen that the proposed approach can also outperform the bencmark model for all the time series except for SA (23.7633 vs 14.6417). For Tokyo time series, it is seen that the MSEs are very large for both models. This is mainly because of one extreme observation which can not be predicted well. For the high quantile prediction, we study the proposed model at the aggregate level. Specifically, we take the log-transformation on both the observations and aggregated fitted values on the training and validation dataset, and calculate the residuals, y i , and t = 1, . . . , 4000. The GPD distribution is then fitted to {e t } t=1,...,4000 with the threshold of .9-quantile of {e t } t=1,...,4000 . The estimated parameters are ξ = 0.3710 with standard error 0.0827, and σ = 0.0981 with standard error 0.0100, which are significant. The QQ-and PP-plots are displayed in Figure 7. It is seen that PP-plot is satisfactory. However, the tails of QQ-plot is off the diagonal line, which is mainly caused by the extreme large values in the attack data. Although there are several points off the diagonal line, we still use the fitted GPD distribution to predict the high quantiles as it does not affect the main conclusion. Algorithm 2 is employed on the attack data at the aggregated level. The assessment of predictedŝ t , t = 4001, . . . , 4512 is shown in Table 5. It is seen that the prediction performance is very satisfactory. The predicted numbers of violations are quite close to the true numbers of violations. They pass all the tests based on LRuc and LRcc tests except that at .98 level of LRcc test is relatively small (0.0241). We conclude that the proposed approach has an overall satisfactory prediction performance.

Honeypot data II with 69 dimensions
In this section, we study the other honeypot data which was used in [22,28]. The dataset was collected by a low-interaction honeypot during 4 November 2010 to 21 December 2010 with a total number of 1123 hours, which has 69 consecutive IP addresses. Each TCP flow initiated by a remote computer and an unsuccessful TCP handshake are deemed as attacks [28]. The data is split into three parts: the training with 750 observations, the validation with 150 observations, and the rest 223 observations are used for the prediction evaluation. The data is log-transformed to reduce the skewness. We employ Algorithm 1 for training and validating the deep learning model, and Algorithm 2 for the prediction. The performance of point predictions are shown in Table 6. It is observed that the proposed deep learning model has a very satisfactory prediction performance. Particularly, the mean (0.091) and median (0.088) of MAPEs by the proposed approach are much smaller than the mean (0.180) and median (0.162) of MAPEs by the benchmark model. For MSEs, the proposed approach also significantly outperforms the benchmark model in terms of mean (0.316 vs 0.718) and median (0.253 vs 0.716). The boxplots are further shown in Figure  8. It is also seen that the proposed approach significantly outperforms the benchmark model.
For the high quantiles, we also study the prediction performance at the aggregate level as it corresponds to the network level in practice where the defense can take place. Specifically, we calculate the residuals, where s t = t i=1 y i ,ŝ t = t i=1ŷ i , and y s are the log-transformed attack data, t = 1, . . . , 900. The GPD distribution is then fitted to {e t } t=1,...,900 with the threshold to be .   displayed in Figure 9. It is seen that PP-plot is satisfactory. The tails of QQ-plot is slightly off the diagonal line. This is again caused by the extreme large values in the attack data. But the overall fitting is fairly well. For the high quantile prediction, Algorithm 2 is employed to the testing data, where the threshold is set to be the .9-quantile of residuals. During the prediction process, we retrain the deep learning model three times by using Algorithm 1 to further improve the prediction accuracy. The p-values of VaR tests are presented in Table 7. We first observe that the proposed model slightly underestimates the number of violations, particularly at the level .95 (11 vs 17), which is mainly due to the extreme large attacks. The other levels have large p-values which indicates that the proposed model can predict them well.
To conclude, the proposed model has overall satisfactory fitting and prediction performances.

Conclusion and discussion
In this work, we propose a novel approach to predicting the multivariate cyber risks. The key idea is to provide the accurate point prediction by training the deep learning network, while employ the extreme value theory for modeling and predicting the high quantiles. One particular advantage of proposed approach is that it can easily handle high dimensional cyber risks thanks to the deep learning architectures, and further, the deep learning model can be retrained to improve the prediction accuracy if needed. The simulation and empirical studies confirm the feasibility of proposed approach and satisfactory fitting and prediction performances. There are some limitations for the current study. Since we model the multivariate dependence via deep learning architectures, the dependence is treated as a 'black-box'. This causes the loss of interpreting the dependence relationship among cyber risks compared to the vine copula approach [22], which is the common issue of deep learning approach. The real attack datasets studied in the current work were collected by honeypots. The other attack data (e.g. collected by telescope or other instruments) may exhibit different phenomena, which should be carefully analyzed before employing the approach developed in our work.
The current work can be extended in several directions. For example, the deep learning model(s) may be further developed for more accurate fitting and prediction [5]. Second, for the complex multivariate cyber risks, the GPD distribution may not be enough for capturing the tail behavior. Then, the non-stationary extreme value distribution may be exploited for the possible modeling [4,20].