Belkacem Berdjane
email: berdjaneb@yahoo.fr

Serguei Pergamenshchikov
email: serge.pergamenchtchikov@univ-rouen.fr

Optimal consumption and investment for markets with random coefficients. *

Keywords: Black-Scholes market, Stochastic volatility, Optimal consumption and Investment, Hamilton-Jacobi-Bellman equation, Feynman -Kac formula, Fixed point solution Mathematical Subject Classification (2000) 91B28, 93E20

We consider an optimal investment and consumption problem for a Black-Scholes financial market with stochastic coefficients driven by a diffusion process. We assume that an agent makes consumption and investment decisions based on CRRA utility functions. The dynamical programming approach leads to an investigation of the Hamilton Jacobi Bellman (HJB) equation which is a highly non linear partial differential equation (PDE) of the second oder. By using the Feynman -Kac representation we prove uniqueness and smoothness of the solution. Moreover, we study the optimal convergence rate of the iterative numerical schemes for both the value function and the optimal portfolio. We show, that in this case, the optimal convergence rate is super geometrical, i.e. is more rapid than any geometrical one. We apply our results to a stochastic volatility financial market.

Introduction

One of the principal questions in mathematical finance is the optimal consumption-investment problem for continuous time market models. This paper deals with an investment problem aiming at optimal consumption during a fixed investment interval [0, T ] in addition to an optimal terminal wealth at maturity T . Such problems are of prime interest for the institutional investor, selling asset funds to their customers, who are entitled to certain payment during the duration of an investment contract and expect a high return at maturity. The classical approach to this problem goes back to Merton [START_REF] Merton | Optimal consumption and portfolio rules in a continuous time model[END_REF] and involves utility functions, more precisely, the expected utility serves as the functional which has to be optimized. By applying results from the stochastic control, explicit solutions have been obtained for financial markets with nonrandom coefficients (see, e.g. [START_REF] Karatzas | Methods of Mathematical finance[END_REF], [START_REF] Korn | Optimal portfolios[END_REF] and references therein). Since then, there has been a growing interest in consumption and investment problems and the Merton problem has been extended in many directions. One of the generalisations considers financial models with random coefficients such as stochastic volatility markets (see, e.g., [START_REF] Fouque | Derivatives in Financial Markets with Stochastic Volatility[END_REF]). In this paper for the CRRA (Constant Relative Risk Aversion) utility functions we consider the optimal consumption-investment problem for a Black-Scholes type model with coefficients depending on a diffusion process which is referred as the external stochastic factor. The pure investment problem for such models is considered in [START_REF] Zariphopoulou | A solution approach to valuation with unhedgeable risk[END_REF] and [START_REF] Pham | Smooth solutions to optimal investment models with stochastic volatilities and portofolio constraints[END_REF]. In these papers the authors use the dynamic programming approach and they show that the nonlinear HJB equation can be transformed in a quasilinear PDE. The similar approach has been used in [START_REF] Kraft | Poerfolio problems stopping at first hitting time with applications to default risk[END_REF] for optimal consumptioninvestment problems with the default risk for financial markets with non random coefficients. Furthermore, in [START_REF] Fleming | An optimal consumption model with stochastic volatility volatility[END_REF], by making use of the Girsanov measure transformation the authors study a pure optimal consumption problem for stochastic volatility financial markets. In [START_REF] Castaneda-Leyva | Optimal consumption investment problems in incomplete markets with stochastic coefficients[END_REF] and [START_REF] Hernández-Hernández | Robust utility maximization in stochastic factor model[END_REF] the authors use dual methods. Usually, the classical existence and uniqueness theorem for the HJB equation is shown by the linear PDE methods (see, for example, chapter VI. [START_REF] Hernández-Hernández | Robust utility maximization in stochastic factor model[END_REF] and appendix E in [START_REF] Fleming | Deterministic and stochastic optimal control[END_REF]). In this paper we use the approach proposed in [START_REF] Delong | Optimal investment and consumption in a black-scholes market with stochastc coefficients driven by a non diffusion process[END_REF] for the optimal consumption-investment problem for financial markets with random coefficients depending on pure jumps processes. Unfortunately, we can not apply directly this method for our case since in [START_REF] Delong | Optimal investment and consumption in a black-scholes market with stochastc coefficients driven by a non diffusion process[END_REF] the HJB equation is the integro-differential equation of the first oder. In our case it is a highly non linear PDE of the second oder. Similarly to [START_REF] Delong | Optimal investment and consumption in a black-scholes market with stochastc coefficients driven by a non diffusion process[END_REF] we study the HJB equation through the Feynman -Kac representation. We introduce a special metric space in which the Feynman -Kac mapping is contracted. Taking this into account we show the fixed-point theorem for this mapping and we show that the fixed-point solution is the classical unique solution for HJB equation in our case. Moreover, using the verification theorem we provide the explicite expressions for the optimal investment and consumption which depend on the HJB solution. Therefore, to calculate the optimal strategies one needs to study numerical schemes for HJB equation. To this end we find an explicite upper bound for the approximation accuracy. Then, we minimize this bound and we get the optimal convergence rate for both the value function and the optimal financial strategies. It turns out that in this case this rate is super geometrical.

The rest of the paper is organized as follows. In section 2 we introduce the financial market, we state the main conditions on the market parameters and we write the HJB equation. In section 3 we state the main results of the paper. Section 4 presents a stochastic volatility model as an example of applications of our results. In Section 5 we study the properties of the Feynman -Kac mapping. In Section 6 the corresponding verification theorem is stated. The proofs of the main results are given in Section 7. In Section 8 we consider a numerical example. In Appendix some auxiliary results are given.

Market model

Let (Ω, F T , (F t ) 0≤t≤T , P) be a standard filtered probability space with two standard independent (F t ) 0≤t≤T adapted Wiener processes (W t ) 0≤t≤T and (V t ) 0≤t≤T taking their values in R d and R m respectively, i.e.

W t = (W 1 t , . . . , W d t ) ′ and V t = (V 1 t , . . . , V m t ) ′ .
The prime ′ denotes the transposition. Our financial market consists of one riskless bond (S 0 (t)) 0≤t≤T and d risky stocks (S i (t)) 0≤t≤T governed by the following equations:

     dS 0 (t) = r(t, Y t ) S 0 (t)dt , dS i (t) = S i (t)µ i (t, Y t )dt + S i (t) d j=1 σ ij (t, Y t ) dW j t , (2.1)
with S 0 (0) = 1 and S i (0) = s i for 1 ≤ i ≤ d. In this model r(t, y) ∈ R + is the riskless interest rate, µ(t, y) = (µ 1 (t, y), . . . , µ d (t, y)) ′ is the vector of stock-appreciation rates and σ(t, y) = (σ ij (t, y)) 1≤i,j≤d is the matrix of stock-volatilities. For all y ∈ R m the coefficients r(•, y), µ(•, y) ∈ R d and σ(•, y) ∈ M d are nonrandom càdlàg functions. Here M d denotes the set of quadratic matrix of order d. Moreover, in just the same way as in [START_REF] Zariphopoulou | A solution approach to valuation with unhedgeable risk[END_REF] we assume, that the stochastic factor Y valued in R m has a dynamics governed by the following stochastic differential equation:

dY t = F (t, Y t ) dt + βdU t , (2.2)
where F is a [0, T ] × R m → R m nonrandom function and β is fixed positive parameter. The process U is the standard Brownian motion defined as

U t = ρV t + 1 -ρ 2 σ * W t , (2.3)
where 0 ≤ ρ ≤ 1 and σ * is a fixed m × d matrix for which σ * σ ′ * = I m . Here I m is the identity matrix of order m.

Moreover, we set K = [0, T ] × R m and we note, that for the model (2.1) the risk premium is the K → R d function defined as

θ(t, y) = σ -1 (t, y)(µ(t, y) -r(t, y) 1 d ) , (2.4)
where

1 d = (1, . . . , 1) ′ ∈ R d .
Similarly to [START_REF] Klüppelberg | Optimal consumption and investment with bounded downside risk for power utility functions[END_REF] we consider the fractional portfolio process

ϕ(t) = (ϕ 1 (t), . . . , ϕ d (t)) ′ ∈ R d ,
i.e. ϕ i (t) represent the fraction of the wealth process X t invested in the i-th stock at the time t. The fractions for consumptions we denote by c = (c t ) 0≤t≤T . In this case the wealth process satisfies the following stochastic equation

dX t = X t (r(t, Y t ) + π ′ t θ(t, Y t ) -c t )dt + X t π ′ t dW t , (2.5)
where π t = σ(t, Y t )ϕ t and the initial endowment X 0 = x. Now we describe the set of all admissible strategies. A portfolio control (financial strategy

) ϑ = (ϑ t ) t≥0 = ((π t , c t )) t≥0 is said to be admissible if it is (F t ) 0≤t≤T -progressively measurable with values in R d × [0, ∞), such that π T := T 0 |π t | 2 dt < ∞ and T 0 c t dt < ∞ a.s. (2.6)
and the equation (2.5) has a unique strong a.s. positive continuous solution (X ϑ t ) 0≤t≤T on [0 , T ]. We denote the set of admissible portfolios controls by V.

In this paper we consider an agent using the CRRA utility function x γ for 0 < γ < 1. The goal is to maximize the expected utilities from the consumption on the time interval [0, T ] and from the terminal wealth. Then for any x ∈ R d , y ∈ R m and ϑ ∈ V the value function of agent is

J(x, y, ϑ) := E x,y T 0 c γ t (X ϑ t ) γ dt + (X ϑ T ) γ ,
were E x,y is the conditional expectation for X ϑ 0 = x and Y 0 = y. Our goal is to maximize this function, i.e. sup ϑ∈V J(x, y, ϑ) .

(2.7)

Remark 2.1. Note, that (2.7) is the classical Merton optimization problem for the market (2.1), (see, e.g. [START_REF] Merton | Optimal consumption and portfolio rules in a continuous time model[END_REF], [START_REF] Delong | Optimal investment and consumption in a black-scholes market with stochastc coefficients driven by a non diffusion process[END_REF]). Pure investment cases of this problem are studied in [START_REF] Zariphopoulou | A solution approach to valuation with unhedgeable risk[END_REF] and [START_REF] Pham | Smooth solutions to optimal investment models with stochastic volatilities and portofolio constraints[END_REF].

Conditions

To list the conditions for the model (2.1) we set

α(t, y) = F (t, y) + β * σ * θ(t, y) , (2.8)
where β * = γ 1ρ 2 β/(1γ). Moreover, we denote by C i,k (K) the space of the functions f (t, y 1 , . . . , y m ) which are i times differentiable with respect to t ∈ [0, T ] and j times differentiable with respect to (y j ) 1≤j≤m .

We assume that the market parameters satisfy the following conditions:

A 1 ) The functions r : K → R + , µ : K → R d + and σ : K → M d belong to C 1,1 (K)
and have bounded derivatives. Moreover, for all (t, y) ∈ K the matrix σ(t, y) is non-degenerated and

sup (t,y)∈K |r(t, y)| + |µ(t, y)| + |σ -1 (t, y)| < ∞ . A 2 ) The K → R m function F (•, •) belongs to C 1,1 (K)
and its partial derivatives are bounded. This condition provides the existance of a unique strong solution for the equation (2.2). On the interval [t, T ] we denote this solution by Y t,y = (Y t,y v ) t≤v≤T with Y t,y t = y.

A 3 ) There exist [0, T ] × R → R continuously differentiable functions (a i,j (•, •)) 1≤i,j≤m such that for any 0 ≤ t ≤ T and for any

x = (x 1 , . . . , x m ) ′ from R m the ith component [α(•, •)] i of the vector α(t, x) has the form [α(t, x)] i = m l=1 a i,l (t, x j ) (2.9)
and, moreover,

α * = max i,j max t,z max |a i,j (t, z)| , ∂a i,j (t, z) ∂t , ∂a i,j (t, z) ∂z < ∞ .
In the sequel, for any vector x from R n we denote the ith component by [x] i . Remark 2.2. In Section 4 we check the condition A 3 ) for a two asset stochastic volatility financial market.

The HJB equation

Now we introduce the HJB equation for the problem (2.7). To this end, for any differentiable K → R function f we denote by D y f (t, y) its gradient with respect to y ∈ R m , i.e. Let now q ∈ R m+1 and M ∈ M m+1 be fixed parameters of the following form

D y f (t, y) = ∂ ∂y 1 f (t, y), . . . , ∂ ∂y m f (t, y) ′ . Let now z(t, ς) be any [0, T ] × R × R m → R two times differentiable function. Here ς = (x, y) ′ , x ∈ R and y ∈ R m . We set D ς z(t, ς) = z x (t, ς) , D y z(t, ς) ′ and D ς,ς z(t, ς) =    ∂ 2 z(t,
q = (q 1 , q) ′ and M = µ ; µ ′ µ ; M 0 , (2.10)
where q 1 , µ ∈ R, q, µ ∈ R m and M 0 ∈ M m . For these parameters with q 1 > 0 we define the Hamilton function as

H(t, ς, q, M) = x r(t, y)q 1 + (F (t, y)) ′ q + 1 γ 1 γ q 1 γ 1 -1 + |θ(t, y)q 1 + β 1 -ρ 2 σ ′ * µ| 2 2|µ| + β 2 2 trM 0 , (2.11)
where γ 1 = (1γ) -1 . In this case the HJB equation is

z t (t, ς) + H(t, ς, D ς z(t, ς), D ς,ς z(t, ς)) = 0 , z(T, ς) = x γ .
(2.12)

To study this equation we make use of the distortion power transformation introduced in [START_REF] Zariphopoulou | A solution approach to valuation with unhedgeable risk[END_REF], i.e. for a positive function h we represent z(t, y) as

z(t, ς) = x γ h ε (t, y) , ε = 1 -γ 1 -ρ 2 γ . (2.13)
It is easy to deduce that the function h satisfies the following quasilinear PDE:

               h t (t, y) +Q(t, y)h(t, y) + (α(t, y)) ′ D y h(t, y) + β 2 2 trD y,y h(t, y) + 1 q * 1 h(t, y) q * -1 = 0 ; h(T, y) = 1 , (2.14)
where

Q(t, y) = γ(1 -ρ 2 γ) 1 -γ r(t, y) + |θ(t, y)| 2 2 (1 -γ) . (2.15)
Note that, by the condition A 1 ) this function and its y-gradient are bounded. Therefore, we can set

Q * = sup (t,y)∈K Q(t, y) and D * = sup (t,y)∈K |D y Q(t, y)| . (2.16)
Our goal is to study the equation (2.14). By making use of the probabilistic representation for the linear PDE (the Feynman-Kac formula) we show in Proposition 5.4, that the solution of this equation is the fixed-point solution for a special mapping of the integral type which will be introduced in the next section.

Main results

First, to study the equation (2.14) we introduce a special functional space. Let X be the set of

K → [1, ∞) functions from C 0,1 (K) such that f 0,1 = sup (t,y)∈K (|f (t, y)| + |D x f (t, y)|) ≤ r * , (3.1)
where

r * = r 0 + m H * q * (2 √ T + T ) + D * T e (α * +Q * )T , (3.2)
r 0 = e Q * T + (e Q * T -1)/Q * q * .
The parameter H * will be defined below in Section 3.1. To define a metrics in X we set

κ = Q * + ζ + 1 + mL * and L * = (1 + r * q * + T D * )e α * T , (3.3)
where the positive parameter ζ will be specified later, q * = (1-ρ 2 γ) -1 . Now for any f and g from X we define the metrics as

̺ * (f, g) = f -g * , (3.4)
where

f * = sup (t,y)∈K e -κ(T -t) (|f (t, y)| + |D x f (t, y)|) .
Moreover, to write the Feynman-Kac representation for the equation (2.14) we need to use the stochastic process (η v ) 0≤v≤T governed by the following stochastic differential equation

dη t = α(t, η t )dt + β dU t . (3.5)
It should be noted that the condition A 3 ) provides the existence and the uniqueness of a strong solution on any time interval [t, T ] and initial condition y from R m . We denote this solution by (η t,y v ) t≤v≤T . Using this process we define the X → X Feynman-Kac mapping L:

L f (t, y) = E G(t, T, y) + 1 q * T t H f (t, s, y) ds , (3.6)
where G(t, s, y) = exp s t Q(u, η t,y u )du and

H f (t, s, y) = E f (s, η t,y s ) 1-q * G(t, s, y) . (3.7)
To solve the equation (2.14) we need to find the fixed-point solution for the mapping L in X , i.e.

L h = h . (3.8)
To this end we construct the following iterated scheme. We set h 0 ≡ 1 and

h n (t, y) = L h n-1 (t, y) for n ≥ 1 . (3.9)
First, we study the behavior of the deviation

∆ n = h(t, y) -h n (t, y) .
In the following theorem we show that the contraction coefficient for the operator (3.6) is given by

λ = 1 + mL * ζ + 1 + mL * (3.10)
and an appropriate choice of ζ gives the super-geometrical convergence rate for the sequence (h n ) n≥1 .

Theorem 3.1. Under the conditions A 1 )-A 3 ) the equation (3.8) has a unique solution h in X such that for any n ≥ 1 and ζ > 0 sup

(t,y)∈K |∆ n (t, y)| + D y ∆ n (t, y) ≤ B * λ n , (3.11)
where

B * = e κT (1 + r * )/(1 -λ) and κ is given in (3.3).
The proof of this theorem is given in Section 7. Now we can minimize the upper bound (3.11) over ζ > 0. Indeed, setting ζ = ζ/(1 + mL * ) and T = (1 + mL * )T , we obtain

B * λ n = C * exp{g n ( ζ)} ,
where

C * = (1 + r * )e (Q * +1+mL * )T and g n (x) = x T -ln x -(n -1) ln(1 + x) .

Now we minimize this function over

x > 0, i.e. min x>0 g n (x) = x * n T -ln x * n -(n -1) ln(1 + x * n ) ,
where

x * n = ( T -n) 2 + 4 T + n -T 2 T . Therefore, for ζ = ζ * n = (1 + mL * )x * n
we obtain the optimal upper bound (3.11).

Corollary 3.2. Under the conditions A 1 )-A 3 ) the equation (3.8) has a unique solution h in X such that for any n ≥ 1 sup

(t,y)∈K |∆ n (t, y)| + D y ∆ n (t, y) ≤ U * n , (3.12)
where

U * n = C * exp{g * n }.
Remark 3.1. One can check directly that for some δ > 0

U * n = O(n -δn ) as n → ∞ .
This means that the convergence rate is more rapid than any geometrical one, i.e. it is super geometrical.

Theorem 3.3. Assume, that the conditions A 1 )-A 3 ) hold. Then the optimal value of J(x, y, ϑ) for optimization problem (2.7) is given by

max ϑ∈V J(x, y, ϑ) = J(x, y, ϑ * ) = x γ (h(0, y)) ε ,
where h(t, y) is the unique solution of the equation (2.14). Moreover, for all 0 ≤ t ≤ T an optimal financial strategy ϑ

* = (π * , c * ) is of the form        π * t = π * (t, Y t ) = θ(t, Y t ) 1 -γ + ε 1 -ρ 2 βσ * D y h(t, Y t ) (1 -γ) h(t, Y t ) ; c * t = c * (t, Y t ) = (h(t, Y t )) -q * .
(3.13)

The optimal wealth process (X * t ) 0≤t≤T satisfies the following stochastic equation

dX * t = a * (t, Y t )X * t dt + X * t (b * (t, Y t )) ′ dW t , X * 0 = x , (3.14)
where

a * (t, y) = |θ(t, y)| 2 1 -γ + ε 1 -ρ 2 β (1 -γ)h(t, y) (σ * D y h(t, y)) ′ θ(t, y) + r(t, y) -(h(t, y)) -q * ; b * (t, y) = θ(t, y) 1 -γ + ε 1 -ρ 2 β (1 -γ)h(t, y) σ * D y h(t, y) .
The proof of this theorem is given in Section 7.

Remark 3.2. Note that the optimal strategy (3.13) coincides with the well-known Merton strategy in the case ρ = 1, i.e. when the process Y is independent of the brownian motion W generating the financial market (2.1). Note also that, in the case 0 ≤ ρ < 1, the optimal investment strategy in (3.13) depends on the external factor Y through the derivative of the function h. The first term in this expression is the well-known Merton strategy and the second term is the impact of the external factor.

To calculate the optimal strategy in (3.13) we use the sequence (h n ) n≥1 , i.e. we set

π * n (t, y) = θ(t, y) 1 -γ + ε 1 -ρ 2 β (1 -γ)h n (t, y) σ * D y h n (t, y)
and c * n (t, y) = (h n (t, y)) -q * . Theorem 3.1-Theorem 3.3 imply the following result Theorem 3.4. Assume, that the conditions A 1 )-A 3 ) hold. Then for any n ≥ 1 sup

(t,y)∈K |π * (t, y) -π * n (t, y)| + |c * (t, y) -c * n (t, y)| ≤ B * 1 U * n ,
where

B * 1 = γ 1 βε 1 -ρ 2 |σ * |r * + q * . Remark 3.3.
Note that in this paper we use only the power utility function x γ with 0 < γ < 1. It seems that the question about the upper bound (3.11) is open for general utility functions of the CRRA type. The method to obtain the bound (3.11) is based on the explicite form of the utility function. The heart of this method is to show that the operator (3.6) is contracted in a suitable functional space. To this end one needs to make use of the power uitility.

Formula for H *

To write the upper bound for the partial derivatives of the function H f (t, s, y) (see Lemma A.2) we define for any q > 0

φ q = √ 2 e T φ q and ι q = qα * (β -2 + T ) , (3.15)
where φ q = q Q * + β 2 α * /2 + (2m 2 + 1)q 2 β -2 α 2 * . Using these parameters we set

H * = 2 max 2φ 1 (ι 1 + α * ) + Ψ * (φ 2 ) 1/2 , 1 β √ 2π e ι 2 , (3.16)
where

Ψ * = T α * + (T D * + T α * (α * + 1 + β 2 /2) + β -2 α * ) 2 and ι = β 2 T (max (ι 1 , ι 2 ) + α * ) 2 .

Two-asset stochastic volatility model

In this section we consider an important example of the model (2.1) used in financial markets with the stochastic volatilities (see, e.g., [START_REF] Pham | Smooth solutions to optimal investment models with stochastic volatilities and portofolio constraints[END_REF]). This model is defined as follows:

dS 0 (t) = r S 0 (t)dt and dS 1 (t) = µ 1 S 1 (t)dt + S 1 (t)σ 1 (t, Y 1 t ) dW 1 t + b 1 dW 2 t ; dS 2 (t) = µ 2 S 2 (t)dt + S 2 (t)σ 2 (t, Y 2 t ) dW 1 t + b 2 dW 2 t ,
where b 1 = b 2 , the process

Y t = (Y 1 t , Y 2 t )
′ is gouverned by the following stochastic differential equations

dY 1 t = F 1 (t, Y 1 t ) dt + βdU 1 t , dY 2 t = F 2 (t, Y 2 t ) dt + βdU 2 t .
Here U i t = ρV i t + 1ρ 2 W i t for some 0 ≤ ρ ≤ 1, i.e. the parameter σ * = 1 in (2.3). Assume, that the [0, T ] × R → R functions F i are bounded and have bounded derivatives, i.e.

F max = max 1≤i≤2 sup (t,y) max |F i (t, y)| , ∂F i (t, y) ∂t , ∂F i (t, y) ∂y < ∞ .
We assume also, that the [0, T ] × R → R functions σ i are such, that Obviously, that in this case the matrix σ(t, y 1 , y 2 ) is non-degenerated and

σ max = max 1≤i≤2 sup (t,
σ -1 (t, y 1 , y 2 ) =      b * 1 σ 1 (t, y 1 ) ; - b * 1 σ 2 (t, y 2 ) - b * 2 σ 1 (t, y 1 ) ; b * 2 σ 2 (t, y 2 )      , where b * 1 = b 2 (b 2 -b 1 ) -1 and b * 2 = (b 2 -b 1 ) -1 . Therefore, the components of the function (2.4) are [θ(t, y 1 , y 2 )] i = b 1,i σ 1 (t, y 1 ) + b 2,i σ 2 (t, y 2 ) , 1 ≤ i ≤ 2 , where b 1,1 = b * 1 (µ 1 -r), b 2,1 = -b * 1 (µ 2 -r), b 1,2 = -b * 2 (µ 1 -r) and b 2,2 = b * 2 (µ 2 -r).
From this we obtain, that in this case function (2.8) has the form (2.9), i.e.

[α(t,

y 1 , y 2 )] i = F i (t, y i ) + β * b 1,i σ 1 (t, y 1 ) + b 2,i σ 2 (t, y 2 ) , 1 ≤ i ≤ 2 .
It is easy to see, that for this model the conditions A 1 )-A 3 ) hold with

α * ≤ F max + b * 1 + σ max σ min and b * = β * max 1≤i,j≤2 (| b 1,i |) .
5 Properties of the mapping L Proposition 5.1. Assume, that the conditions A 1 )-A 3 ) hold. Then L f ∈ X for any f ∈ X .

Proof. Obviously, that for any f ∈ X the mapping L f ≥ 1. Moreover, setting

f s = f (s, η t,y s ) , (5.1)
we represent L f (t, y) as

L f (t, y) = E G(t, T, y) + 1 q * T t E f s 1-q * G(t, s, y)ds . (5.2)
Therefore, taking into account that f s ≥ 1 and q * ≥ 1 we get

L f (t, y) ≤ e Q * (T -t) + T t 1 q * e Q * (s-t) ds ≤ r 0 , (5.3)
where the upper bound r 0 is defined in (3.2). Moreover, Lemmas A.4-A.5 yield

∂ ∂y i L f (t, y) = E ∂ ∂y i G(t, T, y) + 1 q * T t ∂ ∂y i H f (t, s, y) ds .
Therefore, applying here (5.3), (A.16) and Lemmas A.2 we get

L f (t, y) + D y L f (t, y) ≤ r * .
Now we have to show that the function L f (t, y) is continuously differentiable with respect to t for any f ∈ X . Indeed, to this end we consider for any f from X the equation (2.14), i.e.

               u t (t, y) +Q(t, y)u(t, y) + (α(t, y)) ′ D y u(t, y) + β 2 2 trD y,y u(t, y) + 1 q * 1 f (t, y) q * -1 = 0 ; u(T, y) = 1 .
(5.4)

Setting here u(t, y) = u(Tt, y) we obtain a uniformly parabolic equation for u with initial condition u(0, y) = 1. Moreover, due to the conidtions A 1 ) and A 3 ) the functions Q and α have bounded partial derivatives with respect to y. Therefore, for any f from X Theorem 5.1 from [START_REF] Ladyženskaja | Linear and quasilinear equations of parabolic type[END_REF] (p. 320) with 0 < l < 1 provides the existence of the unique solution of (5.4) belonging to C 1,2 (K). Applying the Itô formula to the process u(s, η t,y s ) e

s t Q(v,η t,y v )dv t≤s≤T
and taking into account the equation (5.4) we get

u(t, y) = L f (t, y) . (5.5)
Therefore, the function L f (t, y) ∈ C 1,2 (K), i.e. L f ∈ X for any f ∈ X . Hence Proposition 5.1.

Proposition 5.2. Under the conditions A 1 )-A 3 ) the mapping L is a contraction in the metric space (X , ̺ * ), i.e. for any f , g from X

̺ * (L f , L g ) ≤ λ̺ * (f, g) , (5.6)
where the parameter 0 < λ < 1 is given in (3.10).

Proof. First note that, for any f and g from X and for any

y ∈ R m |L f (t, y) -L g (t, y)| ≤ 1 q * E T t G(t, s, y) f s 1-q * -( g s ) 1-q * ds ≤ T t E G(t, s, y) f s -g s ds .
We recall that f s = f (s, η t,y s ) and g s = g(s, η t,y s ). Taking into account here that G(t, s, y) ≤ e Q * (s-t) we obtain

|L f (t, y) -L g (t, y)| ≤ T t e Q * (s-t) E| f s -g s |ds .
Taking into account in the last inequality, that

| f s -g s | ≤ e κ(T -s) ̺ * (f, g) a.s. , (5.7)
we get

e -κ(T -t) L f (t, y) -L g (t, y) ≤ 1 κ -Q * ̺ * (f, g) . (5.8)
Moreover, by virtue of Lemma A.4 we obtain, that for any 1

≤ i ≤ m ∂ ∂y i L f (t, y) = E ∂ ∂y i G(t, T, y) + 1 -q * q * T t E ( f s ) -q * f 0 (s) ′ υ i (s)G(t, s, y)ds + 1 q * T t E ( f s ) 1-q * ∂ ∂y i G(t, s, y)ds , (5.9)
where f 0 (s) = f 0 (s, η t,y s ), f 0 (s, z) = D z f (s, z) and υ i (s) = ∂ η t,y s /∂ y i . Therefore, for any f and g from X ∂L f (t, y) ∂y i -∂L g (t, y)

∂y i = 1 -q * q * T t E (̟ 1 (s)) ′ υ i (s)G(t, s, y)ds + 1 q * T t E ̟ 2 (s) ∂ ∂y i G(t, s, y)ds ,
where

̟ 1 (s) = f 0 (s) ( f s ) q * - g 0 (s) ( g s ) q * and ̟ 2 (s) = ( f s ) 1-q * -( g s ) 1-q * .
Note, that similarly to (5.7)

| f 0 (s) -g 0 (s)| ≤ e κ(T -s) ̺ * (f, g) a.s. .

Using this, one can check directly that

|̟ 1 (s)| ≤ (1 + r * q * )e κ(T -s) ̺ * (f, g)
and

|̟ 2 (s)| ≤ (q * -1)e κ(T -s) ̺ * (f, g) .
Therefore, these upper bounds and Lemma A.3 imply

∂L f (t, y) ∂y i - ∂L g (t, y) ∂y i ≤ e κ(T -t) L * κ -Q * ̺ * (f, g) ,
where L * is given in (3.3). Therefore, sup

(t,y)∈K e -κ(T -t) D y L f (t, y) -D y L g (t, y) ≤ mL * κ -Q * ̺ * (f, g) .
Taking into account the definition of κ in (3.3), we obtain the inequality (5.6). Hence Proposition 5.2.

Proposition 5.3. Under the conditions A 1 )-A 3 ) the equation L h = h has a uniqe solution in X .

Proof. Indeed, using the contraction of the operator L in X and the definition of the sequence (h n ) n≥1 in (3.9) we get, that for any n ≥ 1

̺ * (h n , h n-1 ) ≤ λ n-1 ̺ * (h 1 , h 0 ) , (5.10)
i.e. the sequnce (3.9) is fundamental in (X , ̺ * ). Therefore, due to Proposition A.1 this sequence has a limit in X , i.e. there exits a function h from X for which

lim n→∞ ̺ * (h, h n ) = 0 .
Moreover, taking into account that h n = L h n-1 we obtain, that for any n ≥ 1

̺ * (h, L h ) ≤ ̺ * (h, h n ) + ̺ * (L h n-1 , L h ) ≤ ̺ * (h, h n ) + λ̺ * (h, h n-1 ) .
The last expression tends to zero as n → ∞. Therefore ̺ * (h, L h ) = 0, i.e. h = L h . Propostion 5.2 implies immediately that this solution is unique.

u = L h = h .
Therefore, the function h satisfies the equation (2.14). Moreover, this solution is unique since h is the unique solution of the equation (3.8).

Remark 5.1. The representation (5.5) for the solution of the equation (5.4) is called the probabilistic representation or the Feynman -Kac formula (see, e. g., [START_REF] Yu | Two-Scale Stochastic Systems[END_REF], p. 194).

Verification theorem

In this section we state the verification theorem from [START_REF] Klüppelberg | Optimal consumption and investment with bounded downside risk for power utility functions[END_REF]. Consider on the interval [0, T ] the stochastic control process given by the Ndimensional Itô process

dς ϑ t = a(t, ς ϑ t , ϑ t ) dt + b(t, ς ϑ t , ϑ t ) dW t , t ≥ 0 , ς ϑ 0 = x ∈ R N , (6.1)
where (W ) 0≤t≤T is a standard k -dimensional Brownian motion. We assume that the control process ϑ takes values in some set Θ. Moreover, we assume that the coefficients a and b satisfy the following conditions

• for all t ∈ [0, T ] the functions a(t, •, •) and b(t, •, •) are continuous on R N × Θ;

• for every deterministic vector υ ∈ Θ the stochastic differential equation dς ϑ t = a(t, ς ϑ t , ϑ) dt + b(t, ς ϑ t , ϑ) dW t has a unique strong solution. Now we introduce admissible control processes for the equation (6.1). We set F t = σ{W u , 0 ≤ u ≤ t} for any 0 < t ≤ T . Definition 6.1. A stochastic control process ϑ = (ϑ t ) 0≤t≤T is called admissible on [0, T ] with respect to equation (6.1) if it is (F t ) 0≤t≤T -progressively measurable with values in Θ, and equation (6.1) has a unique strong a.s. continuous solution (ς ϑ t ) 0≤t≤T such that

T 0 |a(t, ς ϑ t , ϑ t )| + |b(t, ς ϑ t , ϑ t )| 2 dt < ∞ a.s.. (6.2)
We denote by V the set of all admissible control processes with respect to the equation (6.1). Moreover, let f : [0, T ] × R m × Θ → [0, ∞) and h : R m → [0, ∞) be continuous utility functions. We define the cost function by

J(t, x, ϑ) = E t,x T t f (s, ς ϑ s , ϑ s ) ds + h(ς ϑ T ) , 0 ≤ t ≤ T ,
where E t,x is the expectation operator conditional on ς ϑ t = x. Our goal is to solve the optimization problem

J * (t, x) := sup ϑ∈V J(t, x, ϑ) . (6.3)
To this end we introduce the Hamilton function, i.e. for any 0 ≤ t ≤ T , ς, q ∈ R N and symmetric N × N matrix M we set

H(t, ς, q, M) := sup ϑ∈Θ H 0 (t, ς, q, M, ϑ) , (6.4)
where H 0 (t, ς, q, M, ϑ) := a ′ (t, ς, ϑ)q + 1 2 tr bb ′ (t, ς, ϑ)M + f (t, ς, ϑ) .

In order to find the solution to (6.3) we investigate the Hamilton-Jacobi-Bellman equation

   z t (t, ς) + H(t, ς, D ς z(t, ς), D ς,ς z(t, ς)) = 0 , t ∈ [0, T ] , z(T, ς) = h(ς) , ς ∈ R N . (6.5)
Here z t denotes the partial derivative of z with respect to t, D ς z(t, ς) the gradient vector with respect to ς in R N and D ς,ς z(t, ς) denotes the symmetric hessian matrix, that is the matrix of the second order partial derivatives with respect to ς.

We assume that the following conditions hold: H 1 ) The functions f and h are non negative.

H 2 ) There exists [0, T ] × R N → (0, ∞) function z(t, ς) from C 1,2 ([0, T ] × R N )
which satisfies the HJB equation (6.5).

H 3 ) There exists a measurable function ϑ

* : [0, T ] × R N → Θ such that for all 0 ≤ t ≤ T and ς ∈ R N H(t, ς, D ς z(t, ς), D ς,ς z(t, ς)) = H 0 (t, ς, D ς z(t, ς), D ς,ς z(t, ς), ϑ * (t, ς)) .
H 4 )There exists a unique strong solution to the Itô equation

dς * t = a * (t, ς * t ) dt + b * (t, ς * t ) dW t , t ≥ 0 , ς * 0 = x , (6.6
)

where a * (t, •) = a(t, •, ϑ * (t, •)) and b * (t, •) = b(t, •, ϑ * (t, •)).
Moreover, the optimal control process ϑ * t = ϑ * (t, ς * t ) for 0 ≤ t ≤ T belongs to V. H 5 ) There exists some δ > 1 such that for all 0 ≤ t ≤ T and ς

∈ R N sup τ ∈T t E t,ς (z(τ, ς * τ )) δ < ∞
where T t is the set of all stopping times in [t, T ].

Theorem 6.2. Assume that V = ∅ and H 1 ) -H 5 ) hold. Then for all t ∈ [0, T ] and for all x ∈ R N the solution of the Hamilton-Jacobi-Bellman equation (6.5) coincides with the optimal value of the cost function, i.e. z(t, x) = J * (t, x) = J * (t, x, ϑ * ) ,

where the optimal strategy ϑ * is defined in H 3 ) and H 4 ).

Proofs

7.1 Proof of Theorem 3.1

Proposition 5.3 implies the first part of this theorem. Moreover, from (5.10) it is easy to see, that for each n ≥ 1

̺ * (h, h n ) ≤ λ n 1 -λ ̺ * (h 1 , h 0 ) .
Thanks to Proposition 5.1 all the functions h n belong to X , i.e. by the definition of the space X

̺ * (h 1 , h 0 ) ≤ sup (t,y)∈K |h 1 (t, y) -1| + |D y h 1 (t, y)| ≤ 1 + r * .
Taking into account that sup

(t,y)∈K |∆ n (t, y)| + D y ∆ n (t, y) ≤ e κT ̺ * (h, h n ) ,
we obtain the inequality (3.11). Hence Theorem 3.1.

Proof of Theorem 3.3

We apply the Verification Theorem 6.2 to Problem 2.7 for the stochastic control differential equation (2.5). For fixed ϑ = (π, c), where π ∈ R d and c ∈ [0, ∞), the coefficients in model (6.1) are defined as

a(t, ς, ϑ) = x (r(t, y) + π ′ θ(t, y) -c), F (t, y) ′ b(t, ς, ϑ) = x π ′ ; 0 ′ m β 1 -ρ 2 σ * ; βρ I m , where ς = (x, y) ′ ∈ R N , N = m + 1, k = d+ m, 0 m = (0, . . . , 0) ′ ∈ R m , I m is the identity matrix of the order m. Note that bb ′ (t, ς, ϑ) = x 2 |π| 2 ; xβ 1 -ρ 2 π ′ σ ′ * xβ 1 -ρ 2 σ * π ; β 2 I m .
Therefore, according to the definition of H 0 in (6.4), for any q and M of the form (2.10)

H 0 (t, ς,q, M, ϑ) = x r(t, y)q 1 + (F (t, y)) ′ q + (x γ c γ -xcq 1 ) + 1 2 x 2 µ|π| 2 + xπ ′ θ(t, y)q 1 + β 1 -ρ 2 σ ′ * µ + β 2 2 trM 0 .
To check the conditions H 2 ) -H 4 ) we need to calculate the Hamilton function (6.4) for Problem 2.7 which is defined as

H(t, ς, q, M) = sup ϑ∈R d ×[0,∞)
H 0 (t, ς, q, M, ϑ) = H 0 (t, ς, q, M, ϑ 0 ) .

One can check directly, that for µ < 0 and q 1 > 0 we obtain the form (2.11), where the optimal function

ϑ 0 = ϑ 0 (t, ς, q, M) = (π * (t, ς, q, M), c * (t, ς, q, M)) and          π * = π * (t, ς, q, M) = θ(t, y)q 1 + β 1 -ρ 2 σ ′ * µ x|µ| ; c * = c * (t, ς, q, M) = 1 x γ q 1 γ 1 . (7.1)
Proposition 5.4 implies that the HJB equation (2.12) has a solution from C 1,2 (K) defined in (2.13). It should be noted that for this function q 1 = z x (t, x, y) = γx γ-1 h(t, y) > 0 and µ = z xx (t, x, y) = γ(γ -1)x γ-1 h(t, y) < 0. Therefore, in this case the condition H 3 ) holds with ϑ * (t, ς) = (π * (t, ς), c * (t, ς)) and

     π * (t, ς) = θ(t, y) 1 -γ + ε 1 -ρ 2 β σ * D y h(t, y) (1 -γ)h(t, y) ; c * (t, ς) = (h(t, y)) -q * , (7.2)
where the function h is solution of the equation (2.14). Now we check the conditions H 4 ) and H 5 ). Note that the equation (6.6) is identical to the equations (2.2) and (2.5). Due to the condition A 1 ) the coefficients a * (t, y) and b * (t, y) are bounded. Therefore, the equation (2.2) has a unique strong solution which for 0 ≤ t ≤ s ≤ T can be represented as

X * s = X * t e s t a * (v,Y v )dv E t,s ,
where

E t,s = exp s t (b * (v, Y v )) ′ dW v - 1 2 s t |b * (v, Y v )| 2 dv .
It is clear, that for the bounded function b * (v, y) the process (E t,s ) t≤s≤T is a quadratic intergrable martingale and by the Doob inequality

E sup t≤s≤T E 2 t,s ≤ 4 E E 2 t,T < ∞ .
This implies directly, that for any 0 ≤ t ≤ T , x > 0 and

y ∈ R m E sup t≤s≤T X * s |X * t = x , Y t = y < ∞ . Now we recall that z(t, ς) = x γ (h(t, y)) ε
where h(t, y) is some positive function bounded by r * . Thus

E t,ς z m (τ, X * τ , Y τ ) ≤ (r * ) mε E t,ς (X * τ ) γm . Therefore, taking m = 1/γ > 1 one gets E t,ς z m (τ, X * τ , Y τ ) ≤ (r * ) ε/γ E t,ς sup t≤s≤T X * s < ∞
which implies the conditions H 4 ) and H 5 ). Therefore, thanks to Theorem 6.2 we get Theorem 3.3.

Numerical example

In this section through Scilab we calculate the function h(t, y) using the sequence (3.9) with n = 14 iterations. The curve is obtained in the following stochastic volatility market settings: the market consists on one riskless asset (the bond) and a risky one (that means d = 1). Moreover, we set m = 1, T = 1, r(t, y) = 0.01(1 + 0.5 sin(yt)) , µ(t, y) = 0.02(1 + 0.5 sin(yt))

and σ(t, y) = 0.5 + sin 2 (yt). The parameters of the economic factor are F (t, y) = 0.1 sin(yt), β = 1 and ρ = 0.5. The utility function parameter is γ = 0, 75.

1.00

1.05 The second figure illustrates the super geometrical convergence rate for the functions (h n ) n≥0 . The curve represents the accuracy δ n calculated at each step by We observe the values of the accuracy: δ 5 ≃ 0.0001, δ 8 ≃ 10 -08 and δ 14 ≃ 10 -16 . Moreover δ n = 0 ∀n ≥ 15, i.e. numerically the limit function h is reached at the 14 th iteration.

δ n = sup (t,y) |h n (t, y) -h n-1 (t, y)| .

Appendix

A.1 Properties of the space (X , ̺ * ) Proposition A.1. (X , ̺ * ) is a complete metrical space.

Proof. Indeed, let (f n ) n≥1 be a fundamental sequence from (X , ̺ * ), i.e. lim m,n→∞

̺ * (f n , f m ) = 0 .
Taking into account that the norm • * defined in (3.4) is equivalent to the norm (3.1), we obtain that this sequence is fundamental in Banach space C 0,1 (K). Therefore, there exists a function f ∈ C 0,1 (K) such that lim

n→∞ ̺ * (f n , f ) = 0 .
The definition of the metrics ̺ * in (3.4) implies immediately that f ∈ X . Hence Proposition A.1.

A.2 Properties of the function H f

In this subsection we study the smoothness of H f with respect to y.

Lemma A.2. Assume that the condition A 3 ) holds. Then, for any

0 < t ≤ s ≤ T max 1≤i≤m sup y∈R m sup f ∈X ∂ ∂y i H f (t, s, y) ≤ H * 1 √ s -t + 1 , (A.1)
where H * is given in (3.16).

Proof. First, for any y ∈ R m we introduce the auxiliary Brownian motion on the interval [t, s] as

ξ v = y + β(U v -U t ) .
By making use of this process and the Girsanov theorem (see [START_REF] Liptser | Statistics of Random Process I. General Theory[END_REF] p. 254) we can represent the mapping H f as

H f (t, s, y) = E f 1 (s, ξ s ) e Φ(ξ) ,
where

f 1 = f 1-q * , Φ(ξ) = s t Q 1 (v, ξ v ) dv + β -2 m l=1 s t [α(v, ξ v )] l d[ξ v ] l (A.2)
and

Q 1 (v, y) = Q(v, y) - |α(v, y)| 2 2β 2 .
Now we fixe some 1 ≤ i ≤ m and we set

ς v = [ξ v ] i = y i + β ([U v ] i -[U t ] i ) . (A.3)
We recall, that [x] i is the ith component of the vector x ∈ R m . Taking the conditional expectation with respect to ς s we represent H f as

H f (t, s, y) = R H f (s, y, z) p(z, y i )dz , (A.4)
where

H f (s, y, z) = E f 1 (s, ξ s ) e Φ(ξ) |ς s = z , p(z, y i ) = 1 ν √ 2π e -(z-y i ) 2 2ν 2
and ν 2 = β 2 (ts) .

To calculate this conditional expectation note, that one can check directly that for any t < v 1 < . . . < v k < s and for any bounded

R k → R function G E G(ς v 1 , . . . , ς v k )|ς s = z = E G(B v 1 , . . . , B v k ) , (A.5)
where B v is the well-known Brownian Bridge on the interval [t, s] (see, for example, [START_REF] Karatzas | Brownian Motion and Stochastic Calcul[END_REF], p. 359), i.e.

B v = y i + v -t s -t (z -y i ) + ς 0 v - v -t s -t ς 0 s (A.6)
and

ς 0 v = β ([U v ] i -[U t ] i ).
Moreover, using the representation (2.9) we can represent the stochastic integral with respect to the ith component in (A.2) as

s t [α(v, ξ v )] i d[ξ v ] i = s t a i (v, ξ v ) dς v + s t a i,i (v, ς v ) dς v , (A.7)
where

a i (v, x) = m l=1,l =i a i,l (v, x l ) for x = (x 1 , . . . , x m ) ′ .

By putting

A i (v, z, y i ) = z y i a i,i (v, u)du
and using the Ito formula, we can rewrite the last stochastic integral in (A.7) as

s t a i,i (v, ς v ) dς v = A i (s, ς s , y i ) - s t A i (v, ς v , y i )dv , where A i (v, z, y i ) = z y i ∂ ∂v a i,i (v, u) du + β 2 2 ∂ ∂z a i,i (v, z) .
(ξ) = A i (s, ς s , y i ) + Φ 1,i (ξ, y i ) , where Φ 1,i (ξ, y i ) = β -2 s t αi (v, ξ v ) ′ dξ v + β -2 s t a i (v, ξ v )dς v + s t Q 2,i (v, ξ v , y i )dv and Q 2,i (v, x, y i ) = Q 1 (v, x) -A i (v, x i , y i ). Therefore, by the property (A.5) H f (s, y, z) = e A i (s,z,y i ) E f 1 (s, ξ i s ) e Φ 1,i ( ξ i ,y i ) , (A.9)
where

ξ i v = [ξ v ] 1 , . . . , [ξ v ] i-1 , B v , [ξ v ] i+1 . . . , [ξ v ] m ′ . Taking into ac- count that a i (v, ξ i v ) = a i (v, ξ v ) and s t αi (v, ξ i v ) ′ d ξ i v = s t αi (v, ξ i v ) ′ dξ v we get Φ 1,i ( ξ i , y i ) = β -2 s t αi (v, ξ i v ) ′ dξ v + β -2 s t a i (v, ξ v )dB v + s t Q 2,i (v, ξ i v , y i )dv .
Through the definition of B v in (A.6) we rewrite this equality as

Φ 1,i ( ξ i , y i ) = (z -y i -ς 0 s ) β 2 (s -t) s t a i (v, ξ v )dv + s t Q 2,i (v, ξ i v , y i )dv + s t g i (v, ξ i v ) ′ dU v , (A.10)
where the function g i (v, x) ∈ R m is defined as

[g i (v, x)] j = β -1 [α(v, x)] j for j = i ; a i (v, x) for j = i .
Due to the condition A 3 ) this function is bounded, i.e.

|g i (v, x)| ≤ 2β -1 α * . (A.11)
Let us show now, that for any q > 0 Ee qΦ 1,i ( ξ i ,y i ) ≤ φ q e ι q |z-y i | , (A.12)

where φ q and ι q are given in (3.15). First note, that by the condition

A 3 ) | A i (v, z, y i )| ≤ α * |z -y i | + α * β 2 /2
and, therefore,

Q 2,i (v, x) = Q 1 (v, x) -A i (v, z, y i ) ≤ Q * + α * |z -y i | + α * β 2 /2 . Putting ζ i v = g i (v, ξ i v )
and using the coefficients defined in (3.15), we obtain, that for any 0 < t < s < T

Φ 1,i ( ξ i , y i ) ≤ s t ζ i v ′ dU v + α 1 |ς 0 s | + α 2 |z -y i | + α 3 . (A.13)
Obviously, that for each 1 ≤ j ≤ m the process ([ζ i v ] j ) t≤v≤s is adapted to the filtration (B j v ) t≤v≤s , where

B j v = σ ([U u ] j ) t≤u≤v , ( Ǔ j u ) t≤u≤s .
Therefore, for any q > 0 the process

̺ j v = e q v t [ζ i v ] j d[U v ] j -q 2 2 v t [ζ i u ] 2 j du
is martingale, i.e. E ̺ j s = 1. This and (A.11) yield

E e q s t [ζ i v ] j d[U v ] j = E ̺ j s e q 2 2 s t [ζ i v ] 2 j dv ≤ e q 2 β -2 T α 2 * .
Now we use the following multidimensional version of the Hölder inequality. For any integrated variables (η j ) 1≤j≤m

E m j=1 |η j | ≤ m j=1 E |η j | r j 1/r j ,
where (r j ) 1≤j≤m are positive numbers such, that m j=1 r -1 j = 1. This inequality implies

E e q s t g i (v, ξ i v ) dU v = E m j=1 e q s t [ζ i v ] j d[U v ] j ≤ m j=1 Ee qm s t [ζ i v ] j d[U v ] j 1/m ≤ e q 2 β -2 m 2 T α 2 * .
Now from (A.13) for any q > 0 we can estimate the expectation in (A.12) in the following way

Ee qΦ 1,i ( ξ i ) ≤ e q α 2 |z-y i |+q α 3 E e q s t (ζ i v ) ′ dU v +q α 1 |ς 0 s | ≤ e q α 2 |z-y i |+q α 3 E e 2q s t (ζ i v ) ′ dU v 1/2 Ee 2q α 1 |ς 0 s | 1/2 , where α 1 = α * β -2 , α 2 = α * β -2 + T and α 3 = T Q * + β 2 α * /2 .
Taking into account here, that for any a ∈ R

Ee a|ς 0 s | ≤ 2e a 2 β 2 (s-t)/2 ,
we come to the upper bound (A.12). Obviously, that in (A.9) the function

|A i (s, z)| ≤ α * |z -y i |. Therefore, sup f ∈X H f (s, y, z) ≤ φ 1 e (ι 1 +α * )|z-y i | . (A.14)
Moreover, from (A.9) we get

∂ ∂y i H f (s, y, z) = -a i,i (s, y i ) H f (s, y, z) + e A i (s,z,y i ) H 1 f,i (s, y, z) , (A.15)
where A.3 Properties of the process (3.5)

H 1 f,i (s, y, z) = E f 1 (s, ξ i s ) e Φ 1,i ( ξ i ,y i ) ∂ ∂y i Φ 1,i ( ξ i , y i ) . Therefore, sup f ∈X H 1 f,i (s, y, z) ≤ E e 2Φ
In this subsection we study the properties of the process η = (η t,y s ) t≤s≤T Lemma A.3. Under the conditions A 1 )-A 2 ) the process (η t,y s ) t≤s≤T is almost sure continuously differentiable with respect to y ∈ R m for any t ≤ s ≤ T , i.e. for any 1 ≤ i ≤ m there almost sure the derivative υ i (s) = ∂ η t,y s /∂ y i such that One can check directly that the processes υ i (s) satisfies the following differential equations d ds υ i (s) = A s υ i (s) , υ i (t) = e i ,

where A s = α 0 (s, η t,y s ), e i = (0, . . . , 0, 1, 0, , . . . , 0) ′ (only ith component is equal to 1). Now by applying here the Gronwall-Bellman inequality we obtain the upper bounds for the derivatives υ i (s). Hence Lemma A.3.

A.4 Properties of the function G

Now we study the partial derivatives of the function G(t, s, y) defined in (3.6). To this end we need the following general result. This Lemma follows immediately from the Lebesgue dominated convergence theorem.

Lemma A.

Proposition 5 . 4 .

 54 Under the conditions A 1 )-A 3 ) the equation (2.14) has a unique solution which is the solution h of the fixed-point equationL h = h.Proof. Choosing in (5.4) the function f = h and taking into account the representation (5.5) and the equation (3.8) we obtain, that the solution of the equation(5.4)

Figure 1 :

 1 Figure 1: The function h(t,y)

Figure 2 :

 2 Figure 2: The accuracy δ n at each iteration

  Now for any vectorx = (x 1 , . . . , x m ) ′ we set xi = x 1 , . . . , x i-1 , 0, x i+1 , . . . , x m ′ . (A.8)Using these notations one can represent the function (A.2) as Φ

  s)| ≤ e α * T a.s.. Proof. First we introduce the matrix of the first partial derivatives of the function α(v, y) asα 0 (t, z) = ∂[α(t, z)] k ∂z l , 1 ≤ k, i ≤ m .

Lemma A. 4 .

 4 Let F = F (y, ω) be a R × Ω → R random bounded function such that for some nonrandom constant c * d dy F (y, ω) ≤ c * a.s. . Then d dy E F (y, ω) = E d dy F (y, ω) .

  Moreover, through (A.6) and (A.10) we can calculate directly, thatΨ i (t, s, y i ) =where for any t ≤ v ≤ t and any vectorx = (x 1 , . . . , x m ) ′ |Ψ 1,i (v, x)| ≤ D * + α * α * + 1 + β 2 /2 ., where the upper bound Ψ * is given in(3.16). Therefore, φ 1 + Ψ * (φ 2 ) 1/2 and ι * = max (ι 1 , ι 2 /2) + α * . Now from (A.4) we obtain Taking into account here, that for any a > 0 R e a|z-y i | p(z, y i )dz ≤ 2e we obtain the upper bound (A.1). Hence Lemma A.2.

	Ψ 1,i (v, x, y i ) =	∂Q(v, x) ∂x i	-	m j=1	[α(v, x)] j	∂a j,i (v, x i ) ∂z a 2 ν 2 2	-	∂v ∂a i,i (v, x i )
	and	R	-|z -y i | e a|z-y i | β 2 2 ∂ 2 a i,i (v, x i ) ∂z 2 ν 2 p(z, y i )dz ≤ 2a e + ∂a i,i (v, y i ) ∂v	. a 2 ν 2 2	+	√ ν √ 2 π	,
	From the condition A 3 ) we can estimate this function
	max 1≤i≤m	sup (v,x)∈K
	Therefore,
	E Ψ 2 i (t, s) = E	t	s	m j=1 , j =i	∂a j,i (v, B v ) ∂z	2 (s -v) 2 (s -t) 2 dv
			+ E
	sup f ∈X f,Now from (A.14) and (A.15) it follows that H 1
		ι 2 2 |z-y i | ∂a j,i (v, B v ) ∂z (s -v) E Ψ 2 i (t, s) . d[U v ] j ∂ H f (s, y, z) s p(z, y i )dz ∂y i (s -t) (s -v) j=1 , j =i m t sup f ∈X ∂ ∂y i H ∂H f (t, s, y) ∂y i = R + s t (s -t) Ψ 1,i (v, ξ i v , y i ) -a i (v, ξ v ) β 2 (s -t) + R H f (s, y, z) (z -y i ) ν 2 p(z, y i )dz .	dv ,

1

,i ( ξ i ) E Ψ 2 i (t, s) , where Ψ i (t, s, y i ) = ∂Φ 1,i ( ξ i , y i )/∂y i . Taking into account here the bound (A.12) we obtain

sup f ∈X H 1 f,i (s, y, z) ≤ (φ 2 ) 1/2 e s t (sv) (st) Ψ 1,i (v, ξ i v ) -a i (v, ξ v ) β 2 (st) dv 2 . This implies directly E Ψ 2 i (t, s) ≤ Ψ 2 * i (s, y, z) ≤ Ψ * (φ 2 ) 1/2 e ι 2 2 |z-y i | . f (s, y, z) ≤ H * e ι * |z-y i | ,

where H * = α *

  [START_REF] Fouque | Derivatives in Financial Markets with Stochastic Volatility[END_REF]. Under the conditions A 1 )-A 2 ) there exist the partial derivatives (∂G(t, s, y)/∂y i ) 1≤i≤m such that

	max 1≤i≤m ≤ D and sup y∈R m ∂G(t, s, y) ∂y i ∂ ∂y i E G(t, s, y) = E ∂ ∂y i G(t, s, y) .
	Proof. Lemma A.3 implies immediately, that
	∂G(t, s, y)
	∂y

* T e (α * +Q * )T (A.16) i = G(t, s, y)G i (t, s, y) where G i (t, s, y) = s t Q 0 (u, η t,y u ) ′ υ i (u)du, Q 0 (u, z) = D z Q(u, z)

and υ i (u) = ∂ η t,y u /∂ y i . Now Lemma A.3 and Lemma A.4 imply directly this lemma.

The paper is supported by the RFBR-Grant 09-01-00172-a.


Jérôme Franc ¸ois Cynthia Wagner
email: cynthia.wagner@uni.lu

Radu State
email: radu.state@uni.lu

Thomas Engel
email: thomas.engel@uni.lu

DANAK: Finding the odd!

With the growth of network connectivity and network sizes, the interest in traffic classification respectively attack and anomaly detection in network monitoring and security related activities have become very strong. In this paper, a new tool called DANAK has been developed for the detection of anomalies in Netflow records by referring to spatial and temporal information aggregation in combination with Machine Learning techniques. Spatially aggregated Netflow records are fed in a new designed kernel function in order to analyze Netflow records on context and quantitative information. To strengthen the analysis of large volumes of Netflow records, support vector machines are applied. The proposed method has been validated by extensive experimentation on real data sets, including numerous attack strategies of different roots.

I. INTRODUCTION

The need for network security and monitoring activities has become very strong with the growing sizes of networks and the increasing requests for permanent network connectivity. The demand for new and easy-to-use network analysis methods, abstracting complex processes for counter-fighting attacks or detecting network anomalies, is still present. Interesting challenges in networking monitoring are, for example, the automation of security monitoring, like the interaction of different monitoring modules, or to reduce the necessity of interaction between network operator and monitoring framework. Another issue is the large amount of available data that is monitored, where it is nearly impossible to store and evaluate all data. One kind of available information on network borders are Netflow [START_REF] Claise | Cisco Systems Netflow Services Export Version 9[END_REF] records, which can be exported by mostly all routers today, but storing and analyzing these large quantities of Netflow records instantly is a real problem.

A question that arises here is, if it is really necessary to evaluate all Netflow records or if abstractions of Netflow records are sufficient to provide the same outcomes. In recent research, it has often been referred to condensed forms of packets or sampling methods, where only a selected portion of Netflow records were used for evaluation, but a more novel approach is to use an aggregation-based technique.

In this paper, a new monitoring framework called DANAK (Detecting Anomalies in Netflow records by spatial Aggregation and Kernel methods) is described, which aims to detect anomalies and attacks in Netflow records. The analysis of all individual Netflow records is nearly impossible, therefore a spatial and temporal Netflow aggregation based technique is used. For the analysis of the aggregated Netflow records, it has been referred to kernel functions, a subdomain of Machine Learning Techniques.

A new kernel function is described, which captures topological changes and traffic evolution in aggregated Netflow records and by this, deduce if network traffic is anomalous or normal. To strengthen the results of the implemented approach, in addition to the kernel function, a Phase Space Analysis function has been implemented and the different profiles are evaluated by referring to a classification algorithm in order to gain a better understanding of the traffic evolution.

The paper is structured as follows: In section II, the modules of the monitoring framework are presented. The section describes the main components of DANAK and presents theoretical background information about the spatial aggregation of Netflow records, the kernel function and the Phase Space Analysis method. Section III describes the attack injection module as well as the different evaluated attacks and their outcomes. Section IV presents relevant research work in this domain and section V describes future work and presents the conclusions about DANAK.

II. THE ARCHITECTURE

In the following section presents the most relevant theoretical components including background information and the implemented modules of the monitoring framework DANAK. An overview of DANAK is given in Figure 1. The routers, with Netflow record exporting functionality, capture packets on the network and forward the records to the Netflow Collector module, which stores them. After the capture of Netflow records, DANAK comes into play.

The first module of DANAK, the S-aggregation-module, performs the task of spatial -temporal aggregation of Netflow records and is responsible for the generation of traffic profiles. These profiles represent traffic volumes, aggregated by IPaddress prefixes within a time window, having a tree-like representation. Therefore, they are of use to track traffic changes in the network.

The second module, PS-Embedding and Kernel Function, is the main component of DANAK. This module performs a Phase Space Analysis, the processing and the evaluation of the profiles by applying a new kernel calculus model. This method refines profiles by leveraging historical knowledge. The results obtained from the kernel function are applied to a Machine Learning module, the Classifier. This module groups the previously obtained kernel values into anomalies or conventional traffic.

A. Spatial-Temporal Aggregation of Netflow Records

Spatial and temporal aggregation of IP related information was first presented by a tool called Aguri in [START_REF] Cho | An aggregation-based traffic profiler[END_REF], [START_REF] Kaizaki | Detection of Denial of Service attacks using AGURI[END_REF]. Aguri is a near real-time tool, which uses full packet captures as input and then spatially aggregates network data into traffic profiles. A clear advantage of this approach is that due to aggregation, overviews on subnet basis can be presented instead of single IP flow basis. In [START_REF] Cho | An aggregation-based traffic profiler[END_REF], the authors are able to detect different attacks, like flooding or denial of service attacks. In this paper, a similar method to Aguri has been implemented. A major differences to Aguri is, DANAK needs Netflow records and Aguri full packet captures as input format. In Figure 1, this module is denoted by S-Aggregation. The logged Netflow records are spatially aggregated into tree-like structures by using a user-defined aggregation-threshold α over a time interval of η seconds. DANAK outputs traffic profiles in tree-like form by spatially assembling subnet, host and quantitative traffic information by extraction from the captured Netflow records. Netflow records are aggregated into a tree representation, having a Patricia tree as primal structure [START_REF] Morrison | PATRICIA--Practical Algorithm To Retrieve Infromation Coded in Alphanumeric[END_REF], with a bounded size for scalability reason. Each node in the tree represents a subnet. When a record has to be inserted, it is compared with other nodes using a pre-order traversal. Netflow records are compared based on their longest common IP prefix. If a match is found, the volume part (called vol i and expressed in bytes or packets) is updated, otherwise a new branching point with a new node is created. At the end of the aggregation task, only nodes having a volume vol i > α are kept. The output generated by DANAK are traffic profiles for source and destination information. For a detailed description of the aggregation task, we refer the reader to [START_REF] Cho | An aggregation-based traffic profiler[END_REF].

The traffic profiles contain a summary of the contextual and quantitative evolution of network traffic for a small time interval and serve as input information to the kernel function. An extract for a destination traffic profile is represented in Figure 2. This traffic profile is generated in DANAK with an aggregation threshold of α = 5%. The used time interval is η = 5 seconds. The first line in Figure 2 represents the overall transmitted traffic information in bytes for a traffic profile. At the second line starts the generated traffic profile with information about the network traffic. Each line is a node and is divided as follows,

• The first column in a traffic profile line represents the IP-address subnet. • The second column holds the traffic load in bytes.

• The third column represents the percentage of traffic load in a node (including aggregated nodes). • The last column holds the percentage of traffic load of the current subtree.

B. The Kernel Function for Traffic Profiles

For the evaluation of the generated traffic profiles, it has been referred to a commonly used method in Machine Learning, the kernel functions [START_REF] Vapnik | Statistical Learning Theory[END_REF], [START_REF] Schoelkopf | Estimating the Support of a High-Dimensional Distribution[END_REF]. A kernel function can be described as a similarity function for for complex input data and the distances can directly be derived without exhaustive calculations. Vapnik [START_REF] Vapnik | Statistical Learning Theory[END_REF] defines a kernel function K as a similarity mapping K : X × X → [0, ∞[, where X is an input space and K(x, y) = i φ i (x)φ i (y) = φ(x) • φ(y) a similarity score, with φ i (x) being a feature function over a sample x. The more, a kernel function respects symmetry [K(x, y) = K(y, x)] and positive-definiteness.

In this paper, a new kernel function is introduced for the evaluation of the generated traffic profiles. It calculates the similarity between two input profiles, T n and T m , by counting the amount of common nodes and not by performing an exhaustive computation over the entire space [START_REF] Vapnik | Statistical Learning Theory[END_REF]. A schematic figure of the kernel function process can be seen in Figure 3.

For the traffic profile evaluation, the kernel function needs the source and destination IP-addresses and the traffic information as input metrics. The first kernel function metric uses source and destination IP-address information (first column in Figure 4) respecting the CIDR1 scheme [START_REF] Rekhter | An Architecture for IP Address Allocation with CIDR[END_REF]. In the kernel function for profiles, an IP-subnets can be described as, IP = (pref ix i , pref ixlength i ), where pref ix i is the IP network part and pref ixlength i the length of the host identifier part (shown in Figure 4).

! < '() = && ! < ;'#&& ! > ;'# = && ! > '() = && ?@#A&! < *B! < '() C&! < ;'# D& $<;&! > *B! > '() C&! > ;'# D&
The second metric for the kernel function handles the traffic volume information (bytes) for a node, expressed in percent. This percentage represents the traffic volume in a node including other aggregated nodes. This parameter is called vol i (third column in Figure 4). To illustrate the selection of the relevant kernel function parameters, an example is shown in Figure 4. A spatially-aggregated profile can be described as

T i =< T src i , T dst i > with T src i
for a source profile and T dst i for a destination profile. Each profile type is defined as a set of nodes T i = {n 1 , ..., n j }, where a node n x is defined as n x = ( prefix x , prefixlength x , vol x ).

To complete the definition of spatial-temporal aggregation, the temporal aggregation can be defined as the timeline of traffic profiles and can be described as The output of the kernel function about these two profiles K(T n , T m ) gives a similarity score being the sum over all aggregation tree nodes and can be defined as,

K(T n , T m ) = 1 2 i∈T src n ,j∈T src m s(i, j) × v(i, j) + 1 2 i∈T dst n ,j∈T dst m s(i, j) × v(i, j) (1)
with s measuring the similarity between subnets and v catching the similarity between volumes. The similarity score s is the measure for the similitude of traffic profiles. For this, the similarity function compares traffic profiles and calculates the similarity score based on their longest common IP prefix. The higher the similarity score, the more similar the traffic profiles are. The similarity function part s for source and destination traffic profiles can be defined as follows:

s(i, j) =                                  2 pref ixlength j 2 pref ixlength i if pref ix i prefix of pref ix j 2 pref ixlength i 2 pref ixlength j if pref ix j prefix of pref ix i 0 otherwise (2)
The matching function part v is known as a classical Gaussian kernel where σ scales the width of this kernel function. The matching function part v is computed for source and destination traffic profiles and performs the volume calculi. The matching function v can be defined as:

v(i, j) = exp - | vol i -vol j | 2 σ 2 (3)
A good choice of σ is therefore necessary for the results. The optimization of σ can be realized by a grid search [START_REF] Wagner | Game Theory driven monitoring of spatial-aggregated IP-Flow records[END_REF] or by applying the Kruskal algorithm [START_REF]On the shortest spanning subtree of a graph and the traveling salesman problem[END_REF].

C. The Phase Space Method with delayed coordinates

In this section, the Phase Space Analysis with delayed coordinates is briefly presented. Traffic profiles summarize network activity for a period of η seconds. In the case of an attack, it generally operates over more than one time period (η seconds). The same is valid for benign traffic. The principle of the Phase Space Analysis with delayed coordinates allows to reconstruct missing dimensions by using previous and delayed information (function values) as additional coordinates [START_REF] Takens | Detecting Strange Attractors in Turbulence[END_REF], [START_REF] Mandelbrot | The Fractal Geometry of Nature[END_REF], [START_REF] Zalewski | Strange Attractors and TCP/IP Sequence Number Analysis[END_REF]. By this, a historical profile knowledge (including missing dimensions) can be explored.

#
In [START_REF] Zalewski | Strange Attractors and TCP/IP Sequence Number Analysis[END_REF], Phase Space Analysis is explained as a ndimensional space that fully describes a n-variable system, in this case, the generated traffic profiles. A three-dimensional representation of one-dimensional data is generated by calculating the first-order differences between traffic profiles. By this, the function dynamics become more suitable and useful for further processing. The obtained coordinates are also known as function attractors [START_REF] Mandelbrot | The Fractal Geometry of Nature[END_REF], which can be used for further processing in the classifier module.

The Phase Space Analysis with delayed coordinates can be defined as follow: Consider an input set s, and x, y, z as the point coordinates, then the equation for a sample n can be defined as,

x[n] = s[n -2] -s[n -3] (4) y[n] = s[n -1] -s[n -2] (5) z[n] = s[n] -s[n -1] (6)
In this paper, s[n] = K(T n , T n-1 ) and the first-order difference equations are applied to previous, delayed and actual traffic profiles kernel function values in order to obtain an accurate and extended view on the evolution of the network topology and traffic. This allows to catch the dynamism of changes unlike direct historical comparisons of one sample with prior ones. The obtained values are used as input for the evaluation with the classification method.

III. EXPERIMENTAL RESULTS

This section highlights the major outcomes by using DANAK for the classification of traffic profiles. A major problem for the validation of the approach was to find a labeled dataset. A lot of previous works used the Lincoln data set2 , but nowadays it is considered outdated. Therefore in this work, it has been referred to a real data set provided by a large network operator in Luxembourg and the dataset is assumed to be free from malicious network traffic. This assumption is based on the fact that a secondary semi-automated traffic screening has been realized with aid of a network operator specific solution [START_REF]Peakflow SP: Traffic anomaly detection[END_REF].

A technical description of an average sized data set of a duration of 15 minutes is given in Table I.

The parameters for the aggregation task of the Netflow records have been set to a minimal value of α = 5% for the aggregation threshold, in order to have a finest profile granularity and the traffic profile time window set to η = 5 seconds. The generated data set included 180 traffic profiles, where one profile represents a monitoring of 5 seconds.

A. Injection attacks

For the validation of the monitoring framework DANAK, the provided network operator data set has been extended by injecting synthetic attack traces into real traces using the tool Flame [START_REF] Braickhoff | FLAME: a flow-level anomaly modeling engine[END_REF]. Flame is a realistic Netflow-based attack injection tool. The input for the flow reader of Flame are the cleaned Netflow record files from the network operator. In the flow generator and flow deleter, Netflow records can be generated or deleted, depending of the attack typ. The main idea of Flame is to generate traffic by injecting different attacks, but also to affect normal traffic scenarios, as for example in DDoS, where traffic is reduced to the overload. An overview of the Flame architecture is shown in Figure 5.

The injected attacks have been limited to time periods of 5 minutes. This duration is necessary to evaluate the attack detection with various background traffic. In [START_REF] Braickhoff | FLAME: a flow-level anomaly modeling engine[END_REF], different attack models derived from real-world scenarios have been described. In this paper, some of these attack models have been applied to evaluate the performance of DANAK. Additionally to this performance measure, some existing attack models of Flame have been extended in order to create new attack models. The list below shows the attacks that have been generated for the validation of DANAK.

• Nachi scan3 : The Nachi worm first tests the hosts' reachability by ICMP scanning. The attack model considers this, such that one host of the original data set sends single ICMP packets of 92 bytes to destination IP addresses respecting the following properties: a shift of 40 to 45 IP-addresses between two consecutive scans, a positive/negative IP-address shift of 400 each 200 scans and a shift between 45 to 110 IP-addresses each 800 • SSH scan + TCP flood: TCP scans aim to probe SSH servers by trying to log in. Each flow is built of 1 to 4 packets. The inter-arrival time oscillates between 1 to 50 milliseconds. The destination IP-addresses are scanned in a sequential way until 400 scans are approximately executed. Then a shift of 200 to 400 addresses has been observed. In order to simulate a real scenario, 5% of the attacks are considered successful, triggering the victims TCP flood attack.

• DDoS UDP flood: One host receives UDP packets which have been sent from various ports and from multiple IP addresses having different ports. This attack operates by bursts of 40 flows and has a break for 60 to 120 milliseconds.

• DDoS TCP flood: This DoS attack is directed against web servers running on TCP port 80 with 3 packets and 128 bytes. Bursts can be observed 10 packets before a break of 60 to 120 milliseconds. • DDoS UDP flood + traffic deletion: The attack is similar to the DDoS UDP flooding attack, but each additional flow originated by the victim has a deletion probability of 0.2 due to the victim overload.

B. Interpretation of the Kernel Function

In Figure 6, the results of the kernel function for n traffic profiles are graphically represented. In Figure 6, the z- 6) is represented, where a violent DDoS UDP Flood attack has been injected into the data set. The attack has been injected between traffic profiles 60 and 120, which represents a 5 minute attack injection period. It can be observed that the anomalies can clearly be detected by the kernel function and be visualized in the graph. In Figure 7, a realistic DDoS Flood attack has been injected into the data set. Also in this case, the attack occurs between traffic profile 60 and 120, but this time the attack cannot be easily detected by kernel function, as shown in the graph. This result really shows that in this case, the kernel function is not sufficient to identify the anomalies and that it cannot be seen by human expertise per se.

!"!!#$!!% &"!!#$&'% ("!!#$&'% )"!!#$&'% *"!!#$&'% '"!!#$&'% +"!!#$&'% ,"!!#$&'% !% (!% *!% +!% -!% &!!% &(!% &*!% &+!% &-!% !"#$"

C. Interpretation of the Phase Space Embedding Analysis

The following section discusses the utility of using Phase Space Analysis in the detection of anomalies in Netflow records as it was shown that the kernel function (via the z coordinates) by itself, is not strong enough to identify as additional coordinates as illustrated in section II-C. In this paper, it is used to built a historical profile knowledge (including missing dimensions) which shows to be very effective. In Figure 8, a DDoS UDP Flooding attack is visualized without using Phase Space Embedding. The figure represents a simple history of the actual traffic profile T i (expressed as kernel value) compared to the three previous profiles T i-1 , T i-2 and T i-3 . In Figure 9, the same DDoS UDP attack is visualized, but this time by using Phase Space Analysis. The axis in the graph represent the x, y, z -coordinates obtained by applying the first order difference equations, presented in section II-C.

!"!!#$!!% &"!!#$&'% ("!!#$&'% )"!!#$&'% *"!!#$&'% '"!!#$&'% +"!!#$&'% ,"!!#$&'% !% (!% *!% +!% -!% &!!% &(!% &*!% &+!% &-!% !"#$"
By comparing the two figures (Figure 8 and 9), it can been seen that the anomalies can be visualized by applying Phase Space Analysis. In Figure 9, anomalies are now are regrouped and disjoint compared to the case in Figure 8. By combining kernel function and Phase Space Analysis different types of attacks can be visually detected, but to strengthen the validation of the presented approach, a classification algorithm is used in order to identify anomalous and benign network activities. To evaluate the Phase Space Analysis with the kernel function values obtained by the traffic profile processing, the SimpleCART algorithm has been used. The CART algorithm [START_REF] Breiman | Classification and Regression Trees[END_REF] is a supervised learning algorithm and stands for Classification and Regression Trees. It relies on a binary treebuilding algorithm and aims to identify the optimal splitting parameter for the binary splitting. The SimpleCART of WEKA implements a minimal cost pruning method [START_REF] Breiman | Classification and Regression Trees[END_REF]. To fine-tune the experiments, a cross-validation parameter of 5 folds has been set. Cross-validation is the equivalent of dividing the data set into n equally sized sub-data sets (in this case n = 5). The decision tree is generated from subsets where one of the subsets is left out. The performance of these others ones is assessed. This is repeated for all subsets and the average is taken. Each injected attack has been evaluated by using the SimpleCART algorithm. Table II summarizes the classifier results for each attack by the True-Positive and the False-Positive Rate. Finally, it is important to note that, it has been assumed that the data set is free of malicious traffic, but there is no ground truth. Hence, some False Positives may reveal as real attacks in the original data set.

It can be observed that this classifier can accurately detect the different attacks which have been injected into the data set. This highlights that the use of kernel functions for traffic profile evaluation is valid, even if weak attacks could not be detected visually in the graph (see Figure 7). The same conclusions hold for the Phase Space Analysis. The True Positive rates vary around 88% to 94% for the attacks. The False Positive rates are acceptable (2 -20%), even if for one attack, the rate is higher than 30%. The previous paper by Wagner et al. [START_REF] Wagner | Machine Learning Approach for IP-Flow Record Anomaly Detection[END_REF] shows that False Positive rates can be lowered. However in [START_REF] Wagner | Machine Learning Approach for IP-Flow Record Anomaly Detection[END_REF], Netflow records were used as input (without aggregation) and compared in time windows, which is equivalent to compute pairwise distances of approximately 22 000 flows per window. This paper shows the benefit of combining Phase Space Analysis with kernel function to Netflow records, which is only possible if Netflow records are aggregated, as done in DANAK. By inspecting the data set more deeply, it can be observed that the size of traffic profiles varies between 23 to 52 nodes per tree. In the worst case this means to compute 52 2 kernel functions for the DANAK method which is very low in comparison of 22000 2 for the simple Netflow record method. Obviously, the aggregation process needs some additional optimization, as proposed in [START_REF] Cho | An aggregation-based traffic profiler[END_REF], [START_REF] Kaizaki | Detection of Denial of Service attacks using AGURI[END_REF], especially for the bounding of the size of a tree during its construction phase. By introducing tree-optimization, trees can be constructed on the fly. In other words, aggregation is a good tradeoff for improving the scalability of the presented approach by maintaining a high True Positive rate with only a low cost in terms of False Positives.

IV. RELATED WORK

The main interest for using Netflow records in the monitoring task is the compact format and the fact that most commercially available routers support Netflow record export. In the domain of Network monitoring and security, a lot of research has been done yet and a lot of different techniques been presented [START_REF] Bahl | Towards highly reliable Enterprise Network Services via Inference of Multi-level Dependencies[END_REF], [START_REF] Lakhina | Mining Anomalies Using Traffic Feature Distributions[END_REF], [START_REF] Zang | Traffic monitor deployment in IP networks[END_REF], either by introducing dependency models or using statistics. A more recent approach is to refer to rule-based methods, as for example in [START_REF] Karagiannis | BLINC: Multilevel Traffic Classification in the Dark[END_REF], [START_REF] Duffield | Rulebased Anomaly Detection on IP Flows[END_REF]. In [START_REF] Duffield | Rulebased Anomaly Detection on IP Flows[END_REF], the authors use packet based techniques to generated signatures, whereas in [START_REF] Karagiannis | BLINC: Multilevel Traffic Classification in the Dark[END_REF], the authors try to classify traffic by applications.

Netflow records contain relevant information about traffic on the network and the evaluation insights to the operators in case of suspicious events, but require strong storage requirements. In Luxembourg, peak Rates of 60 000 flows /second have been regularly observed at ISP networks supporting high link loads. Therefore, different authors have described the Netflow captures by referring to Netflow sampling [START_REF] Estan | Building a better NetFlow[END_REF], [START_REF] Paredes-Oliva | Portscan Detection with Sampled NetFlow[END_REF], but finding good sampling rates remains challenging and precious information is not considerd. Another interesting approach was presented in [START_REF] Giura | NetStore: An efficient storage infrastructure for network forensics and monitoring[END_REF], where the authors have implemented a new column-oriented storage infrastructure. In [START_REF] Giura | NetStore: An efficient storage infrastructure for network forensics and monitoring[END_REF], the authors have been inspired by spatial and temporal aggregation for their tool. Spatial-temporal aggregation was first presented in [START_REF] Cho | An aggregation-based traffic profiler[END_REF], [START_REF] Kaizaki | Detection of Denial of Service attacks using AGURI[END_REF]. A major advantage of this technique is that key information is preserved and less relevant information aggregated into so called traffic profiles.

Phase Space Analysis has first been in [START_REF] Takens | Detecting Strange Attractors in Turbulence[END_REF] for chaotic dynamical systems to reconstruct attractors from observations by a generic function. This analysis technique is commonly used for analyzing dynamic systems, like nonlinear systems or indeterministic chaos theories [START_REF] Mandelbrot | The Fractal Geometry of Nature[END_REF], [START_REF] Takens | Detecting Strange Attractors in Turbulence[END_REF]. The method of using Phase Space Analysis with delayed coordinates in computer science has first been presented by [START_REF] Zalewski | Strange Attractors and TCP/IP Sequence Number Analysis[END_REF]. In [START_REF] Zalewski | Strange Attractors and TCP/IP Sequence Number Analysis[END_REF], Phase Space Analysis has been used for measuring TCP/IP sequence number generation quality, used for the estimation of attack feasibility and for the Pseudo Random Number Generator function behaviour analysis.

Machine Learning Techniques are accurate approaches for evaluating Netflows or IP flow related data, lass presented in the works for Flow Mining [START_REF] Lee | Mining in a data-flow environment: experience in network intrusion detection[END_REF], [START_REF] Caracas | Mining semantic relations using NetFlow[END_REF]. A sub-domain of Machine Learning Techniques are Kernel methods, which have first been presented in computer science in the 1990's [START_REF] Vapnik | Statistical Learning Theory[END_REF]. An advantage of using kernel methods is that complex data can be evaluated without complex calculus techniques. Besides the kernel functions in computer security, classification techniques are commonly used nowadays for the evaluation of complex data [START_REF] Caracas | Mining semantic relations using NetFlow[END_REF], [START_REF] Kahn | A new intrusion detection system using support vector machines and hierarchical clustering[END_REF], [START_REF] Mcgregor | Flow Clustering using machine learning techniques[END_REF].

V. CONCLUSION

In this paper, a new approach for the evaluation of Netflow records has been introduced by referring to spatial -temporal aggregation. Various attacks on a real data set, obtained from a large network operator in Luxembourg have been analyzed in this new monitoring tool for validating approach. The contribution of DANAK is twofold. First, a new tool based on spatial aggregation has been implemented for processing large quantities of Netflow records. Second, the evaluation of Netflow records in their quantitative and contextual context, where a new kernel function including a Phase Space Analysis for calculating the similarities between traffic profiles has been presented. It has been shown that scalable and efficient flow anomaly detection is possible and corresponding performance metrics well chosen. As future work is planned to optimize the spatial-temporal method for the traffic profile generation and to find kernel function parameters such that attacks can be observed more easily.

Fig. 1 .

 1 Fig. 1. General View of the DANAK Framework Architecture

  !"#$%&'&&(())*+,+& )-)-)-).)&&,/)00)&1/-2(3&.&())-))34&& &&&&25-)-)-).6&&,5)6(0&1,-)23&.&0+-0(34&& &&&&&&&()(-)-)-).*&&,,)**)&1,-))3&.&00-(034&& &&&&&&&&&&()(-(6*-5/-).0)&&+,/*6/&15-*53&.&((-2234&& &&&&&&&&&&&&&()(-(6*-+/-((,.60&&,5/5*0&1,-(63&.&,-(634&& &&&&&&&&&&()(-(+5-(0*-).(2&&,5/*,(&1,-(63&.&,-(634&& &&&&(//-)-)-)./&&++((+)&1+-)(3&.&5+-**34&& &&&&&&&(//-((,-(+5-).0)&&,,055/&1,-)03&.&,-)034&& &&&&&&&(/,-0(6-(60-).00&&,2)60*&1,-653&.&,-6534&& &&&&&&&(/,-0(6-(//-).0)&&+(005*&15-/+3&.&5-/+34&& &&&&&&&(/+-(*5-(0*-).(*&&+6+000&15-+)3&.&(+-0/34&& &&&&&&&&&&(/+-(*5-(//-).0(&&,22,*5&1,-/,3&.&,-/,34&& &&&&&&&&&&(/+-(*5-(5)-).0(&&,5()+/&1,-()3&.&,-()34&& &&&&&&&(/+-(*5-(20-).(*&&,,2,/6&1,-)*3&.&05-+*34&& &&&&&&&&&&(/+-(*5-(20-).0)&&502*5)&1,-+03&.&,-+034&& &&&&&&&&&&(/+-(*5-0)*-).0(&&,5(++6&1,-()3&.&,-()34&& &&&&&&&&&&(/+-(*5-0/)-).0(&&5)**+6&1,-,63&.&,-,634&& &&&&&&&&&&(/+-(*5-0/*-).0(&&,**5(+&1,-6,3&.&,-6,34&&

Fig. 2 .

 2 Fig. 2. Generated Tree of Aggregation Task for a Destination Profile

!

  "#$%&'()&*&+,,-./+&& /0/0/0/12&&+3435&&6202+7&1&2203879& &&2/0320/03&&222.55&&6802471&302+79& &&& &:& !"#$%&;'#&*&+,,-./+&& /0/0/0/1/&&85+8&&6/0247&1&3304579& &&&&&&/0/0/0/12&+4,25&620/5713505+79& & &:& !"#$%&'()&*&2+3.-32&& /0/0/0/1/&&5-.2&&6/0-47&1&330++79& &&&2/0/0/0/13&2.-8-620/-71./05379& & &:& !"#$%&;'#&*&85-,42+&& /0/0/0/1/&&4-+2&&6/0,,7&1&330++79& &2/0/0/0/13&&28.-8&62022714302479& &&&&&&&&&&&&&&&&&&:& K(T n ,T m ) = "

Fig. 3 .

 3 Fig. 3. Traffic Profiles in Tree Form as Input for the Kernel Function

Fig. 4 .

 4 Fig. 4. Parameter Selection for the Kernel Function

{< T src 1 , T dst 1 >

 1 , . . . , < T src M , T dst M >} with T src i being the aggregation tree for time window i for the source information and T dst i for destination information. Consider two traffic profiles denoted by T n =< T src n , T dst n > and T m =< T dst m , T dst m >.

Fig. 5 .

 5 Fig. 5. The Flame Architecture

Fig. 6 .

 6 Fig. 6. Kernel Function Evaluation of a Violent DDoS UDP Flood Attack

Fig. 7 .

 7 Fig. 7. Kernel Function Evaluation of a Weak DDoS UDP Flood Attack

Fig. 8 .

 8 Fig. 8. DDoS UDP Flood Attack without Phase Space Embedding Analysis

Fig. 9 .

 9 Fig. 9. DDoS UDP Flood Attack with Phase Space Embedding Analysis

CIDR: Classless Inter-Domain Routing, Standard scheme for IP-address allocation and IP packet routing. IP addresses are decomposed into network part and host identifier part. Example: 192.168.0.0/16 means 192.168.0.0 is the network part and /16 the routing prefix size.

http://www.ll.mit.edu/mission/communications/ist/corpora/ideval/data/ index.html, last accessed: march 2011

http://www.mcafee.com/threat-intelligence/malware/default.aspx?id= 101013, last accessed: march 2011

http://www.iss.net/security center/reference/vuln/Netbios Name Scan. htm, last accessed: march 2011

http://www.cs.waikato.ac.nz/ ∼ ml/weka/index.html, last accessed: march

ACKNOWLEDGMENT

We especially acknowledge Prof. T. Duhautpas from Restena Luxembourg for his counsel.