Introduction

User modeling is an important task for web applications dealing with large traffic flows. They can be used for a variety of applications such as to predict future situations or classify current states. Furthermore, user modeling can improve detection or mitigation of Distributed Denial of Service (DDoS) attack [1,2,3], improve the quality of service (QoS) [4], individuate click fraud detection and optimize traffic management. In peer-to-peer (P2P) overlay networks, IP models can also be used for optimizing request routing [5]. Those techniques are used by severs for deciding how to manage the actual traffic. In this context, also outlier detection methods are often used if only one class is known. If, for example, an Intrusion Prevention System wants to mitigate DDoS attacks, it usually has only seen the normal traffic class before and it has to detect the outlier class by its different behavior. In this paper we deal with the management of DDos because nowadays it has become a major threat in the internet. Those attacks are done by using a large scaled networks of infected PCs (bots or zombies) that combine their bandwidth and computational power in order to overload a publicly available service and denial it for legal users. Due to the open structure of the internet, all public servers are vulnerable to DDoS attacks. The bots are usually acquired automatically by hackers who use software tools to scan through the network, detecting vulnerabilities and exploiting the target machine. Furthermore, there is also a strong need to mitigate DDoS attacks near the target, which seems to be the only solution to the problem in the current internet infrastructure. The aim of such a protection system is to limit their destabilizing effect on the server through identifying malicious requests. There are multiple strategies with dealing with DDoS attacks. The most effective ones are the near-target filtering solutions. They estimates normal user behavior based on IP packet header information. Then, during an attack the access of outliers is denied. One parameter that all methods have in common is the source IP address of the users. It is the main discriminant for DDoS traffic classification. However, the methods of storing IP addresses and estimating their density in the huge IP address space, are different. In this paper, we present a novel approach based on system identification techniques and, in particular, on the Hammerstein models.

Non-Linear Analytic Prediction of IP Addresses

Data driven identification of mathematical models of physical systems (i.e. nonlinear) starts with representing the systems as a black box. In other terms, while we may have access to the inputs and outputs, the internal mechanisms are totally unknown to us. Once a model type is chosen to represent the system, its parameters are estimated through an optimization algorithm so that eventually the model mimics at a certain level of fidelity the inner mechanism of the nonlinear system or process using its inputs and outputs. These approach is, for instance, widely used in the related big data analytics area (e.g., [6,7,8,9,10,11,12,13,14])

In this work, we consider a particular sub-class of nonlinear predictors: the Linear-in-theparameters (LIP) predictors. LIP predictors are characterized by a linear dependence of the predictor output on the predictor coefficients. Such predictors are inherently stable, and that they can converge to a globally minimum solution (in contrast to other types of nonlinear filters whose cost function may exhibit many local minima) avoiding the undesired possibility of getting stuck in a local minimum. Let us consider a causal, time-invariant, finite-memory,continuous nonlinear predictor as described in (1).

𝑠 ˆ(𝑛) = 𝑓 [𝑠(𝑛 − 1), . . . , 𝑠(𝑛 − 𝑁 )](1)

where 𝑓 [•] is a continuous function, 𝑠(𝑛) is the input signal and 𝑠 ˆ(𝑛) is the predicted sample. We can expand 𝑓 [•] with a series of basis functions 𝑓 𝑖 (𝑛), as shown in (2).

𝑠 ˆ(𝑛) = ∞ ∑︁ 𝑖=1 ℎ(𝑖)𝑓 𝑖 [𝑠(𝑛 − 𝑖)](2)

where ℎ(𝑖) a re proper coefficients. To make (2) realizable we truncate the series to the first 𝑁 terms, thus we obtain

𝑠 ˆ(𝑛) = 𝑁 ∑︁ 𝑖=1 ℎ(𝑖)𝑓 𝑖 [𝑠(𝑛 − 𝑖)](3)

In the general case, a linear-in-the-parameters nonlinear predictor is described by the inputoutput relationship reported in (4).

𝑠 ˆ(𝑛) = 𝐻 ⃗ 𝑇 𝑋 ⃗ (𝑛)(4)

where 𝐻 ⃗ 𝑇 is a row vector containing predictor coefficients and 𝑋 ⃗ (𝑛) is the corresponding column vector whose elements are nonlinear combinations and/or expansions of the input samples.

Linear Predictor

Linear prediction is a well known technique with a long history [15]. Given a time series 𝑋 ⃗ , linear prediction is the optimum approximation of sample 𝑥(𝑛) with a linear combination of the 𝑁 most recent samples. That means that the linear predictor is described as eq. ( 5).

𝑠 ˆ(𝑛) = 𝑁 ∑︁ 𝑖=1 ℎ 1 (𝑖)𝑠(𝑛 − 𝑖)(5)

or in matrix form as

𝑠 ˆ(𝑛) = 𝐻 ⃗ 𝑇 𝑋 ⃗ (𝑛)(6)

where the coefficient and input vectors are reported in ( 7) and (8).

𝐻 ⃗ 𝑇 = [︀ ℎ 1 (1) ℎ 1 (2) . . . ℎ 1 (𝑁 ) ]︀(7)𝑋 ⃗ 𝑇 = [︀ 𝑠(𝑛 − 1) 𝑠(𝑛 − 2) . . . 𝑠(𝑛 − 𝑁 ) ]︀(8)

Non-Linear Predictor based on Volterra Series

As well as Linear Prediction, Non Linear Prediction is the optimum approximation of sample 𝑥(𝑛) with a non linear combination of the 𝑁 most recent samples. Popular nonlinear predictors are based on Volterra series [16]. A Volterra predictor based on a Volterra series truncated to the second term is reported in (9).

𝑥 ˆ(𝑛) = 𝑁 1 ∑︁ 𝑖=1 ℎ 1 (𝑖)𝑥(𝑛 − 𝑖) + 𝑁 2 ∑︁ 𝑖=1 𝑁 2 ∑︁ 𝑗=𝑖 ℎ 2 (𝑖, 𝑗)𝑥(𝑛 − 𝑖)𝑥(𝑛 − 𝑗)(9)

where the symmetry of the Volterra kernel (the ℎ coefficients) is considered. In matrix terms, the Volterra predictor is represented in (10).

𝑠 ˆ(𝑛) = 𝐻 ⃗ 𝑇 𝑋 ⃗ (𝑛)(10)

where the coefficient and input vectors are reported in ( 12) and (12).

𝐻 ⃗ 𝑇 = [︂ ℎ 1 (1) ℎ 1 (2) . . . ℎ 1 (𝑁 ) ℎ 2 (1, 1) ℎ 2 (1, 2) . . . ℎ 2 (𝑁 2 , 𝑁 2 ) ]︂ (11) 𝑋 ⃗ 𝑇 = [︂ 𝑠(𝑛 − 1) 𝑠(𝑛 − 2) . . . 𝑠(𝑛 − 𝑁 1 ) 𝑠 2 (𝑛 − 1) 𝑠(𝑛 − 1)𝑠(𝑛 − 2) . . . 𝑠 2 (𝑛 − 𝑁 2 ) ]︂(12)

Non-Linear Predictor based on Functional Link Artificial Neural Networks (FLANN)

FLANN is a single layer neural network without hidden layer. The nonlinear relationships between input and output are captured through function expansion of the input signal exploiting suitable orthogonal polynomials. Many authors used for examples trigonometric, Legendre and Chebyshev polynomials. However, the most frequently used basis function used in FLANN for function expansion are trigonometric polynomials [17]. The FLANN predictor can be represented by eq.( 13).

𝑠 ˆ(𝑛) = 𝑁 ∑︁ 𝑖=1 ℎ 1 (𝑖)𝑠(𝑛 − 𝑖) + 𝑁 ∑︁ 𝑖=1 𝑁 ∑︁ 𝑗=1 ℎ 2 (𝑖, 𝑗) cos[𝑖𝜋𝑠(𝑛 − 𝑗)]+ 𝑁 ∑︁ 𝑖=1 𝑁 ∑︁ 𝑗=1 ℎ 2 (𝑖, 𝑗) sin[𝑖𝜋𝑠(𝑛 − 𝑗)](13)

Also in this case the Flann predictor can be represented using the matrix form reported in (10).

𝑠 ˆ(𝑛) = 𝐻 ⃗ 𝑇 𝑋 ⃗ (𝑛)(14)

where the coefficient and input vectors of FLANN predictors are reported in (22) and (23).

𝐻 ⃗ 𝑇 = ⎡ ⎣ ℎ 1 (1) ℎ 1 (2) . . . ℎ 1 (𝑁 ) ℎ 2 (1, 1) ℎ 2 (1, 2) . . . ℎ 2 (𝑁 2 , 𝑁 2 ) ℎ 3 (1, 1) ℎ 3 (1, 2) . . . ℎ 3 (𝑁 2 , 𝑁 2 ) ⎤ ⎦ (15) 𝑋 ⃗ 𝑇 = ⎡ ⎣ 𝑠(𝑛 − 1) 𝑠(𝑛 − 2) . . . 𝑠(𝑛 − 𝑁 ) cos[𝜋𝑠(𝑛 − 1)] cos[𝜋𝑠(𝑛 − 2)] . . . . . . cos[𝑁 2 𝜋𝑠(𝑛 − 𝑁 2 )] sin[𝜋𝑠(𝑛 − 1)] sin[𝜋𝑠(𝑛 − 2)] . . . . . . sin[𝑁 2 𝜋𝑥(𝑠 − 𝑁 2 )] ⎤ ⎦(16)

Non-Linear Predictors based on Hammerstein Models

Previous research [18] shown that many real nonlinear systems, spanning from electromechanical systems to audio systems, can be modeled using a static non-linearity. These terms capture the system nonlinearities, in series with a linear function, which capture the system dynamics as shown in Figure 1. Indeed, the front-end of the so called Hammerstein Model is formed by a nonlinear function whose input is the system input. Of course the type of non-linearity depends on the actual physical system to be modeled. The output of the nonlinear function is hidden and is fed as input of the linear function. In the following, we assume that the non-linearity is a finite polynomial expansion, and the linear dynamic is realized with a Finite Impulse Response (FIR) filter. Furthermore, in contrast with [18], we assume a mean error analysis and we postpone the analysis in the robust framework in future work. In other word,

x (n) 2 x (n) 3 x (n) m . . . + p m FIR x(n) y(n) p 2 p 3𝑧(𝑛) = 𝑝(2)𝑥 2 (𝑛) + 𝑝(3)𝑥 3 (𝑛) + . . . 𝑝(𝑚)𝑥 𝑚 (𝑛) = 𝑀 ∑︁ 𝑖=2 𝑝(𝑖)𝑥 𝑖 (𝑛)(17)

On the other hand, the output of the FIR filter is:

𝑦(𝑛) = ℎ 0 (1)𝑧(𝑛 − 1) + ℎ 0 (2)𝑧(𝑛 − 2) + . . . + ℎ 0 (𝑁 )𝑧(𝑛 − 𝑁 ) = = 𝑁 ∑︁ 𝑗=1 ℎ 0 (𝑗)𝑧(𝑛 − 𝑗)(18)

Substituting (17) in (20) we have:

𝑦(𝑛) = 𝑁 ∑︁ 𝑖=1 ℎ 0 (𝑖)𝑧(𝑛 − 𝑖) = 𝑁 ∑︁ 𝑗=1 ℎ 0 (𝑗) 𝑀 ∑︁ 𝑖=2 𝑝(𝑖)𝑥 𝑖 (𝑛 − 𝑗) = 𝑀 ∑︁ 𝑖=2 𝑁 ∑︁ 𝑗=1 ℎ 0 (𝑗)𝑝(𝑖)𝑥 𝑖 (𝑛 − 𝑗)(19)

Setting 𝑐(𝑖, 𝑗) = ℎ 0 (𝑗)𝑝(𝑖) we write

𝑦(𝑛) = 𝑀 ∑︁ 𝑖=2 𝑁 ∑︁ 𝑗=1 𝑐(𝑖, 𝑗)𝑥 𝑖 (𝑛 − 𝑗)(20)

This equation can be written in matrix form as

𝑠 ˆ(𝑛) = 𝐻 ⃗ 𝑇 𝑋 ⃗ (𝑛)(21)

where

𝐻 ⃗ 𝑇 = ⎡ ⎣ 𝑐(2,𝑋 ⃗ 𝑇 = ⎡ ⎣ 𝑠 2 (𝑛 − 2) 𝑠 2 (𝑛 − 3) . . . 𝑠 2 (𝑛 − 𝑁 ) 𝑠 3 (𝑛 − 2) 𝑠 3 (𝑛 − 3) . . . 𝑠 3 (𝑛 − 𝑁 ) 𝑠 𝑀 (𝑛 − 2) 𝑠 𝑀 (𝑛 − 3) . . . 𝑠 𝑀 (𝑛 − 1) 𝑠 𝑀 (𝑛 − 3) . . . 𝑠 𝑀 (𝑛 − 𝑁 ) ⎤ ⎦ (23)

Predictor Parameters Estimation

So far we saw that all the predictors can be expressed, at time instant 𝑛, as

𝑠 ˆ(𝑛) = 𝐻 ⃗ 𝑇 𝑋 ⃗ (𝑛)(24)

with different definitions of the input, 𝑋 ⃗ (𝑛), end parameters vectors 𝐻 ⃗ 𝑇 . There are two well known possibilities for estimating the optimal parameter vector.

Block-based Approach

The Minimum Mean Square estimation is based on the minimization of the mathematical expectation of the squared prediction error 𝑒(𝑛) = 𝑠(𝑛) − 𝑠 ˆ(𝑛)

𝐸[𝑒 2 ] = 𝐸[(𝑠(𝑛) − 𝑠 ˆ(𝑛)) 2 ] = 𝐸[(𝑠(𝑛) − 𝐻 ⃗ 𝑇 𝑋 ⃗ (𝑛)) 2 ](25)

The minimization of ( 25) is obtain by setting to zero the Laplacian of the mathematical expectation of the squared prediction error:

∇ 𝐻 𝐸[𝑒 2 ] = 𝐸[∇ 𝐻 𝑒 2 ] = 𝐸[2𝑒(𝑛)∇ 𝐻 𝑒] = 0(26)

which leads to the well known unique solution

𝐻 ⃗ 𝑜𝑝𝑡 = 𝑅 ⃗ −1 𝑥𝑥 𝑅 ⃗ 𝑠𝑥(27)

where

𝑅 ⃗ 𝑥𝑥 (𝑛) = 𝐸[𝑋 ⃗ (𝑛)𝑋 ⃗ 𝑇 (𝑛)](28)

is the statistical auto-correlation matrix of the input vector 𝑋 ⃗ (𝑛) and

𝑅 ⃗ 𝑠𝑥 (𝑛) = 𝐸[𝑠(𝑛)𝑋 ⃗ (𝑛)](29)

is the statistical cross-correlation vector between the signal 𝑠(𝑛) and the input vector 𝑋 ⃗ (𝑛).

The mathematical expectations of the auto and cross correlation are estimated using

𝑅 ⃗ 𝑥𝑥 (𝑛) = ∑︀ 𝑛 𝑘=1 𝑋 ⃗ (𝑛)𝑋 ⃗ 𝑇 (𝑛) 𝑛(30)

is the statistical auto-correlation matrix of the input vector 𝑋 ⃗ (𝑛) and

𝑅 ⃗ 𝑠𝑥 (𝑛) = ∑︀ 𝑛 𝑘=1 𝑠(𝑘)(𝑛)𝑋 ⃗ (𝑛) 𝑛 (31)

Adaptive Approach

Let us consider a general second order terms of a Volterra predictor

𝑦(𝑛) = 𝑁 −1 ∑︁ 𝑘=0 𝑁 −1 ∑︁ 𝑟=0 ℎ 2 (𝑘, 𝑟)𝑥(𝑛 − 𝑘)𝑥(𝑛 − 𝑟)(32)

It can be generalized for higher order term as

𝑁 ∑︁ 𝑘 1 =1 • • • 𝑁 ∑︁ 𝑘𝑝=1 𝑐 𝑘 1 • • • 𝑐 𝑘𝑝 𝐻 𝑝 {︀ 𝑥 𝑘 1 (𝑛), • • • 𝑥 𝑘𝑝 (𝑛) }︀ (33)

where

𝑁 ∑︁ 𝑘=1 𝑐 𝑘 𝑥 𝑘 (𝑛). (34)

For the sake of simplicity and without loss of generality, we consider a Volterra predictor based on a Volterra series truncated to the second term

𝑟 ˆ(𝑛) = 𝑁 1 ∑︁ 𝑖=1 ℎ 1 (𝑖)𝑟(𝑛 − 𝑖) + 𝑁 2 ∑︁ 𝑖=1 𝑁 2 ∑︁ 𝑗=𝑖 ℎ 2 (𝑖, 𝑗)𝑟(𝑛 − 𝑖)𝑟(𝑛 − 𝑗)(35)

By defining

𝐻 𝑇 (𝑛) = |ℎ 1 (1), • • • , ℎ 1 (𝑁 1 ) , ℎ 2 (1, 1), • • • , ℎ 2 (𝑁 2 , 𝑁 2 )|(36)

and

𝑋 𝑇 (𝑛) = ⃒ ⃒ 𝑟(𝑛 − 1), • • • , 𝑟 (𝑛 − 𝑁 1 ) , 𝑟 2 (𝑛 − 1) 𝑟(𝑛 − 1)𝑟(𝑛 − 2), • • • , 𝑟 2 (𝑛 − 𝑁 2 ) | (37)

Eq (35) can be rewritten as follows

𝑟 ˆ(𝑛) = 𝐻 𝑇 (𝑛)𝑋(𝑛).(38)

In order to estimate the best parameters 𝐻, we consider the following loss function

𝐽 𝑛 (𝐻) = 𝑛 ∑︁ 𝑘=0 𝜆 𝑛−𝑘 [︀ 𝑟 ˆ(𝑘) − 𝐻 𝑇 (𝑛)𝑋(𝑘) ]︀ 2(39)

where 𝜆 𝑛−𝑘 weights the relative importance of each squared error. In order to find the 𝐻 that minimizes the convex function (39) it is enough to impose its gradient to zero, i.e.,

∇ 𝐻 𝐽 𝑛 (𝐻) = 0 (40) That is equivalent to 𝑅 𝑋𝑋 (𝑛)𝐻(𝑛) = 𝑅 𝑟𝑋 (𝑛)(41)

where

𝑅 𝑋𝑋 (𝑛) = ∑︀ 𝑛 𝑘=0 𝜆 𝑛−𝑘 𝑋(𝑘)𝑋 𝑇 (𝑘) 𝑅 𝑟𝑋 (𝑛) = ∑︀ 𝑛 𝑘=0 𝜆 𝑛−𝑘 𝑟(𝑘)𝑋(𝑘)(42)

It follows that the best 𝐻 can be computed by

𝐻(𝑛) = 𝑅 −1 𝑋𝑋 (𝑛)𝑅 𝑟𝑋 (𝑛)(43)

Since

𝑅 𝑋𝑋 (𝑛) = 𝜆𝑅 𝑋𝑋 (𝑛 − 1) + 𝑋(𝑛)𝑋 𝑇 (𝑛)(44)

it follows that

𝑅 −1 𝑋𝑋 (𝑛) = 1 𝜆 [︀ 𝑅 −1 𝑋𝑋 (𝑛 − 1) − 𝑘(𝑛)𝑋 𝑇 (𝑛)𝑅 −1 𝑋𝑋 (𝑛 − 1) ]︀(45)

where 𝑘(𝑛) is equal to

𝑘(𝑛) = 𝑅 −1 𝑋𝑋 (𝑛 − 1)𝑋(𝑛) 𝜆 + 𝑋 𝑇 (𝑛)𝑅 −1 𝑋𝑋 (𝑛 − 1)𝑋(𝑛)(46)

Instead, matrix 𝑅 𝑟𝑋 (𝑛) in ( 43) can be written as

𝑅 𝑟𝑋 (𝑛) = 𝜆𝑅 𝑟𝑋 (𝑛 − 1) + 𝑟(𝑛)𝑋(𝑛)(47)

Thus, inserting Eq (47) and Eq (45)in Eq (43) and rearranging the terms, we obtain

𝐻(𝑛) = 𝐻(𝑛 − 1) + 𝑅 −1 𝑋𝑋 (𝑛)𝑋(𝑛)𝜖(𝑛)(48)

where

𝜖 = 𝑟 ˆ(𝑛) − 𝐻 𝑇 (𝑛 − 1)𝑋(𝑛)(49)

By recalling Eq. ( 46), we can write Eq. (48) as

𝐻(𝑛) = 𝐻(𝑛 − 1) + 𝜖(𝑛)𝑘(𝑛)(50)

By introducing, the new notation,

𝐹 (𝑛) = 𝑆 𝑇 (𝑛 − 1)𝑋(𝑛)(51)

The previous equations can be resumed by the following system

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ 𝐿(𝑛) = 𝑆(𝑛 − 1)𝐹 (𝑛) 𝛽(𝑛) = 𝜆 + 𝐹 𝑇 (𝑛)𝐹 (𝑛) 𝛼(𝑛) = 1 𝛽(𝑛)+ √ 𝜆𝛽(𝑛) 𝑆(𝑛) = 1 √ 𝜆 [︀ 𝑆(𝑛 − 1) − 𝛼(𝑛)𝐿(𝑛)𝐹 𝑇 (𝑛) ]︀ 𝜖(𝑛) = 𝑟 ˆ(𝑛 − 1) − 𝛼(𝑛)𝐿(𝑛)𝐹 𝑇 (𝑛) 𝜖(𝑛) = 𝐻(𝑛 − 1) + 𝐿(𝑛) 𝜖(𝑛) 𝛽(𝑛)(52)

It should be noted that by using Eq (52) the estimation adapts in each step in order to decrease the error. Thus, the system structure is somehow similar to the Kalman filter.

Finally, we define the estimation error as

𝑒(𝑛) = 𝑟(𝑛) − 𝐻 𝑇 (𝑛)𝑋(𝑛)(53)

It is worth noting that the computation of the predicted value from Eq. (38) requires 6𝑁 tot +2𝑁 2 tot operations, where 𝑁 tot = 𝑁 1 + 𝑁 2 (𝑁 2 + 1) /2.

Experiments

In order to prove the effectiveness of the proposed approach, in this section we present our experimental results for a real dataset. Specifically, we consider the requests made to the 1998 World Cup Web site between April 30, 1998 and July 26, 1998 In the dataset, 87 days are reported. We use the first one in order to initialise the estimator in (35) and we use the others as test set by using a rolling horizon method (as in [18]). Particularly, for each day 𝑡 we compute the estimation by using all the IP observations in the previous days [0, 𝑡 − 1]. The results are reported in Figure 2. The increase of errors in June is due to the sudden increment of different IP accessing the website due to the start of the competition (see [19]). It should be noted that the estimation error decrease exponentially despite dealing with several millions of IPs.

Since the computation of the optimal coefficient 𝐻(𝑛) may require some time, we measure the percentage of available data that our approach needs in order to provide good results. Particularly, in this experiment we consider the average estimation error done by the model at time 𝑡 by considering a subset of the IPs observed in interval [0, 𝑡 − 1]. The experimental results on real data cubes are depicted in Figure 3. It should be noted that the Hammerstain model outperform the results by the Volterra's kernel as well as the clustering techniques. In more detail, the clustering techniques are the one less performing. This is due to the nature of the clustering techniques that exploit the geometric information of the data more than their time dependency. We highlight that despite the calculation of 𝐻(𝑛) is computational intensive, this does not effect the real time applicability of the method. In fact, the access decision is taken by considering the estimator 𝑟 ˆ(𝑛) that is computed once per day. Thus the computation of 𝐻(𝑛) does not need to be fast.

Conclusions and Future Work

In this paper, we presented a new way to deal with cyber attack by using Hammerstein models. Experimental results clearly confirm the effectiveness of the proposed techniques for a real data set, outperforming other well-known techniques. Future work will have two objectives. First, we want to consider the problem in a stochastic optimization settings, as for example in [20]. Second, we want to test the approach on other case studies, by also exploiting knowledge management methodologies (e.g., [21,22]).