{'fraction_non_alphanumeric': 0.03602865916069601, 'fraction_numerical': 0.011719549641760491, 'mean_word_length': 4.607173601147776, 'pattern_counts': {'":': 0, '<': 12, '<?xml version=': 0, '>': 10, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 11, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
Statistical Modelling of Computer Network Traffic Event Times

This paper introduces a statistical model for the arrival times of connection events in a computer network. Edges between nodes in a network can be interpreted and modelled as point processes where events in the process indicate information being sent along that edge. A model of normal behaviour can be constructed for each edge in the network by identifying key network user features such as seasonality and self-exciting behaviour, where events typically arise in bursts at particular times of day. When monitoring the network in real time, unusual patterns of activity could indicate the presence of a malicious actor. Four different models for self-exciting behaviour are introduced and compared using data collected from the Imperial College and Los Alamos National Laboratory computer networks.


Introduction
Statistical anomaly detection tools (Lazarevic et al. 2003, Neil et al. 2013 have an important role to play in the next generation of cyber-security defences, complementing more traditional "signature-based" techniques which rely on packet inspection to match known indicators of malicious content held in a database. In contrast, anomaly detection techniques harvest data on normal behaviour in a computer network, and monitor for any significant deviations in observed traffic from probability models built on those historic data. An advantage of anomaly detection techniques is that attack vectors which have not previously been observed may still be detected, and for this reason any departures from normality are of potential interest. We construct a model of a computer's normal behaviour in an enterprise computer network using only the observed event times of computer-computer connections. Lambert et al. (2001) modelled timing patterns in network data using a "dynamic Poisson timing model" with time partitioned into seasonal periods and an exponentially weighted moving average intensity parameter. The approach considered in this article provides a more general counting process framework that captures the self-exciting nature of computer network flow data, where event data can be seen to occur in bursts. This is partly due to human activity, with users having sessions of activity connecting from a particular client to a particular server, and partly due to the data collection process where a stream of packets are aggregated into grouped summaries of packet flows with some arbitrariness in how the records are divided.
Let y 1 , y 2 , . . . be the observed event times in seconds of a counting process Y (t) of con-nections between two computer nodes. The conditional intensity of Y (t) will be assumed to have a general form where µ(t) ≥ 0 is a function which captures seasonal variation in the intensity and For a fixed, non-increasing function ω(u) ≥ 0, ψ i (u) takes one of two forms: Hawkes process: ψ i (u) = ω(u) The process models of Hawkes (1971) and Wold (1948) fundamentally differ in the following way: Under the Hawkes model, all events occurring before time t potentially contribute to the intensity at t; in contrast, defining y * t = max i {y i : y i < t}, the Wold process intensity model simplifies to and so only depends on the time since the most recent event. Both provide potentially plausible models for the burstiness in computer network data, and so their performance will be compared.
The compensator for Y (t) is The structure of the remainder of the article is as follows: Section 2 describes the motivating network traffic flow data, obtained from the computer networks of Imperial College and Los Alamos National Laboratory. Section 3 then introduces a piecewise constant model for the inherent seasonalities present in computer network traffic. Section 4 presents four models for self-exciting behaviour using Hawkes and Wold processes. Section 5 describes how changepoint detection algorithms can be used to estimate the parameters of both the seasonal and self-exciting components of the models. The performance of these approaches to modelling computer network traffic are compared in Section 6. Section 7 presents an alternative model for discrete event time data, and demonstrates performance over an entire computer network.

Computer network data
Two sources of computer network traffic data are considered: network flow ("NetFlow") data from the Imperial College computer network, and authentication logs obtained from the enterprise network of Los Alamos Laboratory (LANL) . NetFlow records are a high level aggregation of the packets sent between two IP addresses over the same ports, under the same protocol in a short space of time, and summarise what might be interpreted as a single communication between the two addresses. NetFlow data are difficult to anonymise, and so cannot be publicly shared.
The authentication log data have been anonymised and are publicly available , and specify the times at which users performs authentication events on computers within the internal network. An event type is recorded that indicates which kind of authentication event occurred, such as a network log-on, interactive log-on, a workstation screen lock, and so on. Turcotte et al. (2017) and Price-Williams et al. (2017) noted that some authentication event types exhibit strong periodic patterns, indicating automated network activity. For this analysis, we only consider event types that are considered "interactive" and can be linked to a user being present at their computer. College IP address to a particular internet sever over a 4-week period. The right image plots the "interactive log-on" events, which were recorded in seconds, for a random user (user 689229) from the LANL network. It is apparent that the event times in both examples exhibit seasonal variations, with all events occurring within the typical hours of a working day, and clear two day breaks that correspond to weekends. More interestingly, particularly in the NetFlow data, the events are seen to occur in bursts, even when diurnal patterns have been taken into account.

Modelling seasonality
Seasonal patterns appear in computer network traffic data since human users are more likely to be active on weekdays during the day time than they are at night time or at the weekend.
In the counting process conditional intensity (1), these seasonal variations are represented by the function µ(t). For seasons of total duration S, models for µ(t) can be assumed to be repeated in blocks of time S, such that µ(t) ≡ µ(t mod S) for all t > 0. Assuming a weekly seasonality of S = 604,800 seconds, we therefore simply require a model for µ(t) on Within a single week there will be daily fluctuations in activity, with weekend days in particular likely to witness fewer network events than weekdays. However, besides these fluctuations in volume, it will be further assumed that the activity patterns each day will have the shape. Formally, we assume where S = S/7 is the duration of one day in seconds, x is the integer part of x,μ(t) is a probability density function on [0, S ] andμ 0 , . . . ,μ 6 are positive daily multipliers such that the seasonal component of the conditional intensity (1) on the ith day of the week at the time of day t seconds isμ(t)μ i .
Given some historical event data, the parametersμ 0 , . . . ,μ 6 can be simply estimated by the average number of events occurring on each week day, and so it simply remains to estimate the density functionμ(t) of event arrivals throughout the day. For computational tractability and simplicity,μ(t) will be assumed to be piecewise constant with an unknown number of changepoints m, denoted σ 1 , . . . , σ m , ordered such that 0 = σ 0 < σ 1 < . . . < σ m < σ m+1 = S ; letμ 0 ,μ 1 , . . . ,μ m be the corresponding densities in each changepoint segment, such that and m i=0μ i (σ i+1 − σ i ) = 1. Estimation of (4) will be considered later in Section 5.

Modelling self-exciting behaviour
In the general counting process conditional intensity (1), self-exciting behaviour is represented by the function ψ(u) and the scalar λ. Given an estimate for µ(t), modelling selfexciting behaviour can be simplified using the time-rescaled process where M −1 is the inverse of (2). For each event time which further simplifies to for the Hawkes model and The choice of the non-increasing excitation function ω controls the increase in a user's intensity following each event and the subsequent rate of decay. Two alternatives are con-sidered: where α, β > 0 in (8), and (9) is a step function with changepoints 0 ≡ τ 0 < τ 1 < . . . < τ and decreasing step heights λ 0 > . . . > λ −1 > λ ≡ 0. Substituting either of the two excitation functions (8) and (9) into the Hawkes and Wold models (6) and (7) gives rise to four different models for self-exciting behaviour. Cartoon illustrations of these models are shown in Figure 2.

Parameter estimation
Given the event times z 1 , z 2 , . . . , z n of the time-rescaled counting process (5) Numerical maximum likelihood estimation is straightforward for the exponential excitation function (8), since in this case there are just two unknown parameters under either the Hawkes or Wold models. More details are provided in Ozaki (1979), Laub et al. (2015).
In contrast, the non-increasing step function excitation model (9) (6), and the right column corresponds to the Wold model (7). The excitation function ω is either an exponential decay (top row) or a step function (bottom row). changepoints and their locations can be successfully estimated using a changepoint algorithm, using the Bayesian information criterion (BIC) (Schwarz 1978) to prevent over-fitting the number of parameters. This will be considered in detail in Section 5.
Under the full Hawkes model, even applying changepoint approaches is too complex.
Instead, a simpler version of the Hawkes-step function model will be considered, where the number of steps is fixed to be 1. In this case, simple numerical optimisation of (10) can be deployed.

Changepoint detection
Estimation of the changepoints for both the seasonal density (4) and for the piecewise constant excitation function (9) under the Wold model can be efficiently performed using pruned exact linear time (PELT) method of Killick et al. (2012). For an observed sequence of random variables x 1:n , the algorithm finds changepoints 0 ≡ τ 0 < τ 1 < . . . < τ m < τ m+1 ≡ n which minimise the Bayesian Information Criterion (Schwarz 1978 where L is an estimated likelihood function and α > 0 notionally represents the number of additional free parameters introduced to the model by adding a changepoint. When x 1:n are the event times from an inhomogeneous Poisson process, and α = 2, since introducing a new changepoint adds two parameters: the location of the changepoint, and the intensity within the new segment. To estimate the seasonal density, historical event times from Y (t) are mapped onto [0, S ] and changepoint detection for (4) proceeds with PELT by treating these as a Poisson process.
For estimating the excitation step function (9) in the Wold model, the event times z 1:n from the transformed process Z(t) (5) must be further transformed into Poisson process event times before the same estimation procedure can be used. Let d i = z i+1 − z i , i = 1, . . . , n − 1 be the waiting time between successive events and let d (1) , . . . , d (n−1) be the corresponding order statistics. Defining δ i = i j=1 (n − j + 1)d (n−j) , then under the excitation model (9), δ 1 , . . . , δ n−1 are event times from an inhomogeneous Poisson process with conditional intensity at t given by Letτ 1:ˆ ,λ andλ 1:ˆ +1 be the maximum likelihood changepoints and intensities subject to the constraint thatλ 1:ˆ +1 are non-increasing. These can be calculated recursively, To ensure the BIC estimate of the intensity (11) is still non-increasing, the PELT algorithm is constrained here to only allow changepoints which are a subsequence from the full, constrained maximum likelihood solutionτ 1:ˆ .

Modelling NetFlow event times
The conditional intensity function (1), assuming either the Hawkes (6) or Wold (7) models, with exponential (8) or piecewise constant (9) excitation functions are now fit to the NetFlow event times from Figure 1, representing the times an Imperial College computer contacted an internet server over a 28-day period. Under each model, the parameters are estimated using the first two weeks of data, and then goodness-of-fit is compared on the second two weeks of data.
To examine the effect of modelling seasonality, it is first assumed that ψ i (u) = 0 in (1), and that the seasonal densityμ(t) is given either by (4) or else it is constant and equal to 1/S. The maximum likelihood estimate for the parameter λ is given simply by n/T , the average number of events observed per unit of time.
If the model for the conditional intensity (1) is correct, then the time rescaling theorem (Brown et al. 2002) states that Λ(y 1 ), . . . , Λ(y n ) are the event times of a homogeneous Poisson process with unit rate, where Λ(t) is the compensator function (3). The inter-arrival times of this time-transformed process, Λ(y i ) − Λ(y i−1 ), i = 1, . . . , n with y 0 ≡ 0, should therefore be exponentially distributed with rate 1. Under the null hypothesis of the intensity model being correct, the upper tail p-value is given by Under the intensity model, p i ∼ Uniform (0, 1), while large (small) p-values suggest smaller (larger) than expected waiting times between event. A KolmogorovSmirnov (KS) test (Massey Jr 1951) can be used to assess the goodness of fit of the event times to the model by comparing the empirical cumulative distribution function (CDF) of p 1:n to the CDF of a Uniform(0, 1) distribution. Table 1 shows the KS test statistic for each of the proposed models for the conditional intensity (1). For a visual assessment, Q-Q plots of the p-values p 1:n are presented in Figure 3.  The seasonal baseline model provides only a limited improvement over the assumption of a constant intensity, suggesting that capturing seasonal variability alone is not sufficient for realistic modelling of the arrival times of NetFlow events; the shape of the Q-Q curve for the seasonal model in Figure 3 implies the p-values are too large, meaning the majority of the events arrive more quickly than expected. Incorporating any of the four models for capturing self-excitation causes a significant jump in performance. The Wold step function model narrowly outperforms the Hawkes exponential method, possibly because the number of parameters is not fixed and is estimated from the data using a changepoint detection algorithm, providing greater flexibility.

Modelling computer network data without seasonality
In some situations it may not be feasible to usefully model the seasonal component of a user's conditional intensity (1); for example, if a user has irregular work patterns. Furthermore, it is possible for the self-exciting component to largely compensate for a poorly specified seasonal model, such as µ(t) ≡ 1, as the self-exciting terms will naturally increase the intensity function more during seasonal periods where underlying activity is higher.

Discrete event times
In some real world cyber-security data sets, event times may be recorded quite crudely.
For example, the authentication logs from the Los Alamos National Laboratory enterprise network introduced in Section 2 have event times recorded only to the nearest second. Particularly when the event times are bursty, this can lead to several events appearing to be observed at the same moment in time, and so a Poisson process model is no longer appropriate. One solution could be to bin all events that appear in the same second and treat them as one combined event in any further analysis. A problem with this solution, particularly in anomaly detection, is that a bursty sequence of events may actually represent anomalous behaviour and binning the events within one second may weaken or remove the anomalous signal. An alternative approach considered here is to explicitly model the arrivals as a discrete time process.
Following on from Section 6, attention is restricted here to the best performing model from that section, which was the Wold model with self-excitation described by a step function, with no seasonal component. This model is particularly straightforward to convert into discrete time. In continuous time, under the Wold model for Y (t) with step function self-excitation, the conditional intensity function simplifies to where λ 1 > λ 2 > . . . > λ > 0. This translates to piecewise exponentially distributed waiting times between events. To recast this model in discrete time, first the changepoints τ i can be constrained to be integer-valued; second, since the geometric distribution is the discrete analogue for the exponential distribution (both share the lack of memory property), it will be assumed that on entering segment S i = {τ i , . . . , τ i+1 − 1} without having observed an event, the next event occurs in S i with probability and if so, the event time follows a truncated geometric distribution on S i . This formulation allows an unbounded number of events to occur at each time point. The parameters of (13) are estimated using the same methodology as in Section 5, with the likelihood function adapted to the revised model (14).

Analysis of user authentication event data
As mentioned above, the authentication data from the Los Alamos National Laboratory (LANL) computer network, Section 2, are only recorded to the nearest second with many coincidental values, and so modelling must proceed in discrete time. For a comparison with the proposed approach, a baseline homogeneous model is also implemented with a discrete hazard function which is constant over time.
For each user in the network, the first 28 days of data are used to estimate the parameters of both models. The performances are then tested on the remaining 62 days of data. The analysis is restricted to the 3119 users in the LANL network with at least 200 events in both sections of data. For each user, the Kolmogorov-Smirnov (KS) test statistic is calculated to measure how well the two discrete models fit the test data.

Conclusion
A general framework for specifying the conditional intensity for normal, user-driven computer network behaviour has been presented, capturing key user behavioural features such as seasonality and self-exciting burstiness. Such models have important applications in anomaly detection in statistical cyber-security, (Neil et al. 2013, Turcotte 2013, potentially decreasing the number of false positive detections through more realistic models of normal network and user behaviour.
A Wold model for self-excitement, using step functions to model decay in excitation, was shown to perform very well across a large sample of users, and negated the need for modelling any seasonal variations in the data, such as diurnal activity patterns. The Hawkes exponential and Hawkes step function models could be used to detect anomalies along each edge of a computer network by monitoring the peak magnitude of the estimated conditional intensity; if this intensity rose above a nominated threshold value, a flag could be recorded indicating a possible attack with unusually high activity.

Supplementary Material
Supplementary materials avaliable online at https://github.com/Matt0312/SToCND contain python code to implement the methods introduced in the article.