Enrico Bernardi
email: en-rico.bernardi@unibo.it

Silvia Romagnoli
email: sil-via.romagnoli@unibo.it

Computing the Volume of an high-dimensional semi-unsupervised Hierarchical copula

Keywords: copula volume, Hierarchical copula, clustering methods, high dimensional problems, MatLab code copula volume, Hierarchical copula, clustering unsupervised method, high dimensional problem, MatLab code

We propose an algorithm for the computation of the volume of a multivariate copula function (and the probability distribution of the counting variable linked to this multidimensional copula function) which is very complex for large dimensions. As it is common practice for large dimensional problem, we restrict to positive orthant dependence and we construct a Hierarchical copula which describes the joint distribution of random variables accounting for dependence among them. This approach approximates a multivariate distribution function of heterogenous variables with a distribution of a fixed number of homogenous clusters, organized through a semi-unsupervised clustering method. These clusters, representing the second level sectors of Hierarchical copula function, are characterized by an into-sector dependence parameter determined by a method which is very similar to the Diversity Score method. The algorithm, implemented in Mat-Lab code, is particularly efficient allowing us to treat cases with a large number of variables, as can be seen in our scalability analysis. As an application, we study the problem of valuing the risk exposure of an insurance company, given the marginals i.e. the risks of each policy.

Copula functions represent a methodology which has recently become the most significant new tool to handle in a flexible way the co-movement between markets, risk factors and other relevant variables studied in finance.

A n-dimesional copula is a function C from the n-dimensional hypercube to [0, 1] such that: i) C is grounded; ii) its one-dimensional margins are the identity function; and iii) C is n-increasing, (see [START_REF] Nelsen | Introduction to Copulas[END_REF][START_REF] Joe | Multivariate Models and Dependence Concepts[END_REF]. The last requirement is strictly linked to the notion of copula's volume; in effect it is equivalent to require the positivity of the copula's volume between two points. In [START_REF] Cherubini | Computing the Volume of n-Dimensional Copulas[END_REF] an algorithm is proposed to compute analytically such volume; since it must evaluate a set of 2 n terms, that algorithm becomes rapidly unusable for problems of dimension higher than 20. Our main goal here is to present an algorithm allowing to treat large dimensional problems. This algorithm is composed by two parts: the first one is the preparing phase while the second is the computational one. In the preparing phase we reduce the dimension of the problem. Our algorithm organizes via a semi-unsupervised method (K-means) the variables into a fixed number of groups which will be then considered as clusters of homogeneous variables. Such approximation will be more accurate as the number of clusters that we decide to set, increases. The dependence parameters in the groups, or the number of variables to consider will be determined by a procedure which is in the same spirit of the Diversity Score approach. After this preparation procedure, the computational phase aggregates the groups with a Hierarchical Archimedean copula function since we restrict to the positive orthant dependence, as it's usual for the high dimensional problem. Along the same lines, [START_REF] Savu | Hierarchical Archimedean copulas[END_REF] propose a class of Hierarchical Archimedean copulas to model high-dimensional joint distributions of asset returns organized in supervised clusters. For the applications of this problem we recall that the computation of the survival copula, frequent in finance because it corresponds to the put-call parity in multidimensional setting (Cherubini et all., 2004), is a particular case of copula volume's computation. Furthermore we observe that the volume is involved in the computation of the distribution function of the counting variable linked to the multidimensional copula. The paper is organized as follows: in section 2 we present and discuss the preparing phase or the clustering method and the homogenous approximation, in section 3 we analyze the computational phase from a theoretical point of view and finally in section 4 we propose an empirical application and a scalability analysis of our algorithm. Section 5 concludes. In the Appendix we present the MatLab code for our procedures.

The preparing phase

Here we explain the clustering procedure of the original sample and the homogenous approximation, used in order to reduce the complexity of the problem.

The clustering method

The primary objective of cluster analysis is to partition a given data of multidimensional vectors (patterns) into so-called homogeneous clusters such that patterns within a cluster are more similar to each other than patterns belonging to different clusters. The best known and most widely used algorithm to reach this goal is K-means. Lately neural networks, for example, competitive-learning networks, self-organizing feature maps and adaptive resonance theory networks have been also often used to cluster data. Each approach has its own pros and cons. In the application chosen in section 4 we opted for K-means but the method proposed may be implemented similarly grouping the data with other clustering methods. K-means clustering is a partitioning method, i.e. the function K-means partitions the observations in our data into K mutually exclusive clusters and returns a vector of indices pointing to the clusters which have been assigned to each observation. K-means treats each observation in our data as an object having a location in space. It finds a partition in which objects within each cluster are as close to each other as possible, and as far of objects in other clusters as possible. Each cluster in the partition is defined by its member objects and by its centroid, i.e. the point to which the sum of distances from all objects in that cluster is minimal. K-means computes cluster centroids differently for each distance measure, to minimize the sum with respect to the chosen measure (it is the mean of the points in each cluster for our choice of distance measure). K-means then uses an iterative algorithm that minimizes the sum of distances from each object to its cluster centroid, over all clusters. This algorithm moves objects between clusters until the sum cannot be decreased further. The result is a set of clusters that are as compact and well-separated as possible. K-means is a semi-unsupervised method since it requires to choose a distance metric and a number of clusters further discussed in the following subsections. Moreover in section 4, we propose a comparison between our algorithm implemented with K-means and with a completely unsupervised clustering method, i.e. with Kohonen self-organizing maps (SOM).

Metric selection

In order to mathematically identify clusters in a data set, it is usually necessary to first define a measure of similarity or proximity which will establish a rule for assigning patterns to the domain of a particular cluster centered. As it is to be expected, one of the most popular is the Euclidean distance which is one of the most common Minkowski distance metric. By using the Euclidean distance, K-means algorithm tends to detect hypersphericalshaped clusters. Since clusters can be of arbitrary shapes and size, the Minkowsky metrics may not be a good choice a priori. In [START_REF] Su | A modified version of the K-means algorithm with a a distance based on cluster symmetry[END_REF] On the other hand other researchers have focussed more on the performance of a metric in a randomized data set, than in the choice of a particular metric, arguing that in the case of absence of noise of any kind, virtually all metrics are equivalent. In Francois et all. (2005) it is suggested to choose the metric according to the shape of the noise that is assumed on the data. While it is known that the Euclidean metric is optimal in presence of white Gaussian noise, other type of noise require other metric such as the fractional norm which is proved to be a better dissimilarity estimator than the Euclidean norm when a highly colored noise is present. Giving the kind of data that we are clustering in this problem, we choose to work with the squared Euclidean distance as a distance measure. In section 4 we discuss some experiments on the performance of this choice for some randomized data.

Cluster Number Selection

The choice of the number of groups, given the dimensionality of the problem, plays an important role since it affects the between and within dependence parameters and then also the precision of the method. This is a critical problem which involves the consideration of the trade-off between the complexity of the problem in term of time of computing and the precision of the method. Several papers address the problem of cluster number selection by using a K-means approach that exploits local changes of internal validity indices to split or merge clusters. In [START_REF] Pelleg | X-means: Extending K-means with efficient estimation of the number of clusters[END_REF] this issue is addressed by developing a dynamic K-means variant known as X-means that performs bisecting splits of clusters where the resulting clustering solution is addressed by the Bayesian Information Criterion. Here the bisecting K-means algorithm splits the data samples into two disjoint clusters and if the split increases the overall fitness measured by internal validity indices, the cluster is split and the bisecting K-means continues recursively. If none of the bisecting steps of each existing cluster leads to an improved result, the algorithm stops and returns the current clusters as result. In Muhr and Granitzer (2009) a comparison is proposed between an X-means method with several internal multivariate validity indices that measure cluster quality based on the total scatter of between-cluster and within-cluster sum of squares, as the Calinski and Harabasz index (see [START_REF] Calinski | A dendrite method for cluster analysis[END_REF], the Hartingan index (see [START_REF] Hartigan | Statistical theory in clustering[END_REF] and the Krzanowski and Lai index (see [START_REF] Krzanowski | A criterion for determining the number of groups in a dataset using sum of squares clustering[END_REF]. To further emphasize how the implementation of the K-means technique applies to our problem, we propose to determine the best cardinality of groups given a subjective tolerance δ, defined as the minimum number of groups such that the maximum distance, on the chosen metric, of each data to the relative centroid is not greater than the tolerance. In this way, if the distance of a data is greater than the tolerance, one more group is made. Formally, if we work on the Euclidean metric, we can define the δ-optimal number of groups, n δ , as the argument of the following optimization problem with constrain

n δ = arg min u j -c i ≤δ,∀i,j (CER(i))
where u j are the data, c i are the centroids and the clustering error (CER) to minimize is

CER(i) = i w=1 n i j=1 (u j -c i ) 2
where n i is the number of data in the i-th group. We can observe that CER is very similar to the clustering error defined as sum of confusions in Wang et all. (2006). We can notice that the selection of n δ may be included in the K-means procedure with a simple if/else code-block.

The Homogenous approximation

In this subsection we discuss the homogeneous and heterogeneous group, so it is useful to first specify what is the meaning in this context that the elements inside the group are or are not equal in value.

The clustering method assembles some groups which must be considered homogenous to simplify the computation in case of high dimensional problems. To this aim we develop a procedure allowing us to consider that the variables are homogenous. In order to do that we must compensate the homogenous assumption with a change in the dependence parameter inside the group or in the size of the group. This problem is similar to that of To determine the dependence parameters we propose a method along the lines of the Diversity Score approach, which is inapplicable in our particular setting though. Embrechts et al. ( 2003) have shown that there must be a higher degree of dependence for variates linked at a lower level, than between those linked only at higher level to satisfy the conditions for a proper multidimensional distribution.

For this reason we can't transform the groups into homogenous and independent samples but we must consider only the homogenous transformations.

In [START_REF] Rouvinez | Diversity scoring for market value CDOs[END_REF] the "equivalent number of assets" is defined, which is an alternative method to the Diversity score one, where the independent request is dropped.

If we compensate with the dependence parameter, we have the following;

Definition 1 The Equivalent Dependence Parameter (EQDP) for the copula function C of group i, i = 1, ..., n is α i such that:

C(u i 1 , u i 2 , ..., u i m , α i ) = C(c i , ..., c i , α i )
where u i j , j = 1, ..., m are the heterogenous marginals of the group i, C(c i , ..., c i , α i ) is an m-variate copula with homogeneous marginals c i equal to the centroid of the i-th cluster and α i is the standard estimation of the dependence parameter within the heterogenous group for the same copula function.

On the other hand, if we decide to change the number of the variables to be considered into the homogenous group, we have the following; Definition 2 The Equivalent Number of Variables (ENV) for the copula function C of group i, i = 1, ..., n is the dimension m i of copula C after the homogeneous approximation: where u i j , j = 1, ..., m are the heterogenous marginals of group i, C(c i , ..., c i , α i ) is an m i -variate copula with homogeneous marginals c i equal to the centroid of the i-th cluster and α i is the standard estimation of the dependence parameter within the heterogenous group for the same copula function.

C(u i 1 , u i 2 , ..., u i m , α i ) = C(c i , ..., c i , α i ) 7

The Computational phase

We recall that the volume of copula function in n dimensions may be represented in a combinatorial form which is well suited for computational purposes [START_REF] Cherubini | Computing the Volume of n-Dimensional Copulas[END_REF]. Given a copula function C, the volume of the n-dimensional box S = [u, v] with u, v ∈ [0, 1] n , u ≤ v, may be represented as:

V C (S) = 2 n -1 i=0 (-1) k(i) C(c n (p n (i)))
where p n (i) denotes the n-dimensional vector of zeros and ones corresponding to the binary representation of i, k(i) is the number of elements equal to one in p n (i), and c n (p n (i)) is an n-dimensional vector such that c j = v j if p j = 0 and c j = u j if p j = 1, where u j , v j and c j denote the j-th element of the corresponding vector. We refer to the previous representation where we introduce a Hierarchical copula function. These copulas extend the class of Archimedean copulas to account for non-exchangeability [START_REF] Joe | Multivariate Models and Dependence Concepts[END_REF]

HC(u 1 , u 2 , ..., u n ) ≡ φ -1 s (φ s ♦φ -1 1 (φ 1 (u 1 )+φ 1 (u 2 ))+...+φ s ♦φ -1 n (φ n (u n-1 )+φ n (u n )))
where φ s represents the generator of the Archimedean copula representing the dependence between groups, φ j , j = 1, ..., n represents the dependence within groups and ♦ is the compound operator. We can apply the previous representation of the volume of copula function to the Hierarchical class. u, v ∈ [0, 1] n , u ≤ v, may be represented as:

V HC (S) = n i=0 (-1) i D c (i,k) j=1 HC(c(p i,k (j)))
where the coordinates defining the box S are obtained by the matrix ũ, ṽ ∈

[0, 1] k × [0, 1]

Remark 4

We observe that D(i, k) can be determined by assuming to have a number of groups equal to n, i.e. to the sum of groups'cardinalities, each one with cardinality equal to 1, so that one can take into account the order of the groups. Then, computing D(i, k) = n i , we must filter these o.c.d. with respect to the groups' cardinalities, reducing to the compatible ones, i.e. D c (i, k).

For example, if we have two groups each one with cardinality three, the number of o.c.ds of i = 1 is D(1, 2) = D c (1, 2) = 6 since all the o.c.ds are compatibles with the groups' cardinalities. For all i we assume to order, in some way, the o.c.c.ds according to some order index j. This way, the second o.c.c.d. for i = 1 and two groups, each one with cardinality three, corresponds (for the chosen ordering criteria which is the row's order here) to one. We observe that we can choose any ordering criteria since the final order of the o.c.c.ds is not relevant.

p 1,2 (2)
We can observe that the problem is simplified assuming that the groups are homogeneous.

Proposition 5 (Hierarchical Homogeneous copula volume) Given a Hierarchical copula function HHC with k homogeneous groups of dimension m s , s = 1, .., k respectively such that

k s=1 m s = n, the volume of the k- dimensional box S = [u, v] with u, v ∈ [0, 1] k , u ≤ v,
may be represented as:

V HHC (S) = n i=0 (-1) i Dc (i,k) j=1 HHC(c(p i,k (j)))
where Dc (i, k) denotes the number of compatible combinatorial distribution (c.c.ds in the following) 3 of the integer i into k groups, p i,k (j) is a k × max(m s , s = 1, ..., k) matrix of zeros and ones corresponding to the j-th c.c.d. of i ones into k groups of dimensions m s , s = 1, ..., k, c(p i,k (j)) is a k × max(m s , s = 1, ..., k) matrix such that c s,w = v s if p s,w = 0 and c s,w = u s if p s,w = 1, where u s , v s are the s-th element of the corresponding vector and c s,w denote the (s, w)-th element of the corresponding matrix and HHC(c(p i (j))) is the hierarchical copula computed for the j-th c.c.d. accounting for the permutation of the elements into the groups.

Remark 6

We can observe that D(i, k) can be computed by a number of permutations with repetitions, i.e.

D(i, k) = n + k -1 k -1
For example if we assume two groups each one with cardinality three, the number of c. If we order these c.c.ds, then the second c.c.d. for i = 1 distributed into two groups each one with cardinality three, corresponds (for the chosen ordering criteria given here by the column's order) to

ds of i = 1 is D(1, 2) = Dc (1, 2) =
p 1,2 (2) = 0 0 0 1 0 0
We can observe that the computation of the Survival copula function is a particular case of the volume computation;

Corollary 7 (Survival copula function in n dimension) In the same setting of Proposition 5, the survival of the Hierarchical copula function HHC with marginals u ∈ [0, 1] k , may be represented as:

V HHC (S) = n i=0 (-1) i Dc (i,k) j=1 HHC(c(p i,k (j)))
where the box S = [u, e] and where e = (1, 1, ..., 1) is the unit vector in R k .

The volume of a copula is involved into the computation of the probability distribution of the counting variable linked to copula function;

Definition 8 (Counting variable linked to n-dimensional copula function)

In the same setting of Proposition 5, the counting variable linked to the copula function HHC(c(p l,k (j)) is given by

r(j) = (p l,k (j))
which counts the number of elements equal to one in the matrix p l,k (j) for the j-th c.c.d.. Clearly we have r(j) = l, ∀j. We can recover the distribution of the counting variable linked to the ndimensional copula function. To do this we observe that a number l must correspond, through his probability function, to a number in [0, 1] by considering all the c.c.ds of l ones into k groups (we enumerate this c.c.d. with the index j). The j-th c.c.d., corresponding to the matrix p l,k (j), generates a pair of coordinates4 , i.e. the j-th c.c.d. generates the box S j for which we compute the volume. The sum of the volumes computed for all the coordinates generated by all the c.c.ds for the same number of ones, i.e. for the same l ones into k groups, represents the probability to count l ones, or better to have l events whose probabilities are the marginal centroids.

Corollary 9 (Probability distribution of the counting variable linked to the n-dimensional copula) In the same setting of Proposition 5, the probability distribution of the counting variable linked to the n-dimensional copula HHC, with k homogeneous groups of dimension m s , s = 1, .., k respectively such that k s=1 m s = n and with centroids u ∈ [0, 1] k , is the function

F r : [0, n] → [0, 1] such that: F r (i) = i l=0 Dc (l,k) j=1 V HHC (S j (p l,k (j)))
where V HHC (S j ) is the volume of the hierarchical homogeneous copula computed for the box S j = [u j , v j ] ∈ R k × R k , determined for the j-th c.c.d. of l ones into k groups such that, if p s,w = 0, u j s = 0 and v j s = u s , if p s,w = 1, u j s = u s and v j s = 1, where p l,k (j) is a k × max(m s , s = 1, ..., k) matrix of zeros and ones corresponding to the j-th c.c.d. of l ones into k groups of dimensions m s , s = 1, ..., k and finally accounting for the permutation of the elements into the groups.

Referring to the representation of copula's volume of Proposition 5, in the Appendix we lay out the computational phase detailing the MatLab code of our algorithm.

Copula's volume vs Probability distribution

We observe that we have found all the elements which may be used to compute the copula's volume but also the distribution of counting variable linked to the multidimensional copula function. In this case, each row of the matrix 0-1 (p 4,3 [START_REF] Cherubini | Bivariate Option Pricing with Copulas[END_REF] in the example in the Appendix), generates the coordinates of a volume corresponding to a probability and then, by aggregation, corresponding to an element of the probability distribution of the counting variable (see corollary 9). Finally we remark that the computation of the survival copula function is a particular case of volume computation (see Corollary 7).

Empirical Application

We consider here an insurance company, that for simplicity operates only in the life branch. The company evaluates the marginal probabilities of death of each insured subject but the total probability of default isn't simply assessable for an high number of policies. We face here the problem of determining the survival copula function corresponding to the probability of total survivals. We observe that for large dimensional problems like the one we are dealing with, it's necessary to reduce its complexity and we propose here the application of a semi-unsupervised clustering method. In this example we consider 80 policies. K-mean clustering method based on the Euclidean metric with tolerance δ = 0.00039, generates 5 groups with the following cardinalities: clusterV ector = 11 31 30 6 2 represented by the histogram in Figure 2. Figure 1 represents the marginal death probabilities of the sample and underlines the clusters created by Kmeans approach. These groups are assumed homogeneous with marginals equal to the centroids of the same group, which represents the vector u in the Corollary 5: The dependence parameters within and between the groups are assumed equal to the EQDP in the sense specified before 5 : α within = 0.1 2 3 0.4 0.8 α between = 0.01

In this applications we find a Survival copula function equal to 0.9121 as a sum of 80 terms with sign represented in Figure 4 6 . We can fix the number of groups and the dependence parameters determined before and calculate the Survival copula function for a sample where we take off an element with marginal equal to the last, following the order depicted in Figure 1. Working backwards till a cardinality of the sample equal to 40, we can observe in Figure 4 the behavior of the Survival copula 5 We assume to know the estimations of the real parameters within and between the groups which may be estimated by an econometric methodology, like as the maximum likelihood method, possibly simplifying the problem by considering the limiting distributions proposed in [START_REF] Schönbucher | Taken to the limit: simple and not-sosimple loan loss distributions[END_REF] for Archimedean copulas with granularity adjustment, and then compensated to account for the homogeneous approximation as suggested in Definition 1. 6 We must recall the representation of the Survival copula function proposed in Corollary 7. The terms represented in Figure 3 function with respect to the size of the sample. As expected, the Survival copula function decreases for larger samples.

In order to validate our procedure, we make a comparison with the algorithm CR proposed in Cherubini and Romagnoli (2009) that we consider as our Benchmark, i.e. the truth. We consider the maximum dimension which allow to calculate the truth, i.e. we consider a sample of 20 heterogeneous elements (which correspond to the first 20 elements of our sample), and compare it to the hierarchical homogeneous survival (HH Survival) calculated with our procedure for several number of groups. The clusters are made with K-means and with Kohonen self-organizing map (SOM) and the distance is computed on the Euclidean metric. In table 1 we can observe the performance of our K-means procedure which clearly increases with the number of the groups.

In table 2 and 3 we investigate the performance of our procedure in a randomized data set. We start with the same data set used in table 1 randomized by a white Gaussian noise and by a colored Gaussian noise, i.e. a white Gaussian noise affecting only few marginals. The clusters are made by K-means based on the Euclidean distance. We compare the CR Survival with the K-means HH Survival, computed on a randomized data by a Monte Carlo procedure. In a random content we compute the K-means HH Survival for a maximum number of groups equal to 5 since K-means can't preserve a greater number of clusters during the iterations, reducing it when someone become empty. We can observe that for an higher number of groups, the error decreases for the white randomization but increases for the colored one. In line with Francois et all. (2005), we can assert that our approach, and therefore the K-means approach based on the Euclidean metric, can be a good choice in the presence of white Gaussian noise but the performance is not so good if we assume a colored randomization of data.

In table 4 we compare the CR Survival with the Kohonen SOM HH Survival computed on the 20-th dimensional data set, allowing also for white and colored randomization. Kohonen SOM groups the data on a map with six neurons and, yet after 200 iterations, the grid tends to approximate the data distribution. In the random data experiments, we compute SOM HH Survival by a Monte Carlo procedure.

If we compare table 4 with the previous tables relative to K-means applications, we observe that SOM with a small number of nodes behave in a way that is similar to K-means while, as it's pointed in Kohonen An interesting question concerns the scalability of the algorithm, i.e. the maximum dimension to which our techniques can be extended. In K-means HH algorithm the scalability must be analyzed with respect to the number of groups and the total dimension of the problem. The MatLab algorithm was lunched on a Dell Dimension DXP061 workstation, Intel(R) 2CPU, 6400@213 GHz, 213 GHz, 2.00 GB RAM. Table 5 reports the computing time for the first 5 percent of the probability distribution of a variable of dimension n = 100 for different numbers of groups. We propose a procedure to compute the volume of a multivariate copula function in high dimensional problems. Our approach is articulated in two steps:

we cluster, by a semi-unsupervised method, the variables in a fixed number of homogeneous groups by compensating with the into-sectors and between-sectors dependence parameters to take account for the homogeneous assumption);

we compute compute a Hierarchical copula function describing the joint distribution of random variables and its volume.

The empirical example presented above explains the strong potentiality of this procedure, as can be observed in the scalability analysis reported, and points also to the feasible aspect of the algorithm strongly dependent on the memory stack resources of MatLab program.

A Appendix

A.1 The combinatorial problem

To compute the volume of copula function, we must solve an interesting combinatorial problem. In fact we are concerned with the computing of the c.c.ds of an integer into the groups already determined by the clustering procedure. We must consider the c.c.ds for all the integers from zero to the dimension of the problem.

1 plist = zeros(0,n); 2 for i = 0:m 3 temp = partitions(total sumi * c, ... % If the number of ones in the group is bigger than its 7 % cardinality, the row is set to zero and then erased.

A.3 Preparing the Copula computation

Giving the c.c.ds' matrix, we compute the copula function into the groups and so we create a matrix of dimension equal to the number of groups × the maximum dimension of the groups (this is the matrix p i,k (j) in Proposition 5) where we distribute i ones into k groups according to the j-th c.c.d.. For example if the second row of the c.c.ds' matrix for the integer 4 into 3 groups with the same dimension equal to 3, is the following:

Hierarchical Copula computation

We compute the copula function into the groups; we aggregate the marginal probabilities with an Archimedean copula (in the following MatLab code we compute the Clayton's copula function with parameter equal to the EQDP of the corresponding group). This must be done for each row of the c.c.d. matrix.

1 % This function computes the Clayton's hierarchical second level 2 % copula, and the product of binomial coefficients allows to 3 % take account for the position into the groups. 4 5 % alpha is the vector of the parameters inside the groups 6 % alpha 0 is the parameters among groups 7 % numones is the c.c.ds' matrix for the same number of 8 % deaths into the groups 9 % sizecluster is the vector of the groups' cardinalities 10 % u and v are the vectors of volume's coordinates ,alpha,alpha 0) 14 numclusters=size(numones,2);

  a new non-metric measure is proposed, based on the idea of symmetry which can be used to group a given data set into a set of clusters of different geometrical structures. In Wang et all. (2006) through observing the characteristic property of data clustering, a new data-dependent dissimilarity measure is designed, called density-sensitive distance metric, which has the property of elongating the distance among points in different high density region and then can reflect the characters of data clustering. The density-sensitive K-means algorithm can identify non-convex clustering structures, thus generalizing the application area of the original K-means algorithm. In the same spirit Xing et all. (2002) determines the metric as distance which solves a convex-optimization problem for a given data set.

5 Page

 5

  portfolio of collateral to a synthetic pool of identical assets, each defaulting independently, in the analysis of cash flow collateralized debt obligation. The Diversity Score, which is defined as the number of assets in the synthetic pool, was firstly proposed in[START_REF] Cifuentes | The binomial expansion method applied to CBO/CLO analysis[END_REF] and following in the Moody's alternative version in[START_REF] Cifuentes | Buying and selling credit risk: a perspective on credit-linked obligations[END_REF] and some applications of the model was discussed inSchorin and Weinreich (1999) and in Duffie and[START_REF] Gârleanu | Risk and valuation of collateralized debt obligations[END_REF].

Page 7

 7 of 27 URL: http://mc.manuscriptcentral.com/gcom E-mail: ijcm@informa.

Proposition 3 ( 8 Page 8

 388 Hierarchical copula volume) Given a Hierarchical copula function HC with k groups of dimension m s , s = 1, .., k respectively such that k s=1 m s = n, the volume of the n-dimensional box S = [u, v] with of 27 URL: http://mc.manuscriptcentral.com/gcom E-mail: ijcm@informa.

  row corresponds to the first group and the second to the second

2 1 .

 21 since all the c.ds are compatibles with the groups' cardinalities. In this framework these two c.c.ds correspond to the following situations: the integer 1 is placed in the first group, independently of the position (for convention in the first place);2. the integer 1 is placed in the second group, independently of the position (for convention in the first place).

11 Page

 11

Figure 1 :

 1 Figure 1: Marginal Death probabilities of the sample with evidence of Kmeans clusters.

Figure 2 :

 2 Figure 2: K-means group's cardinality histogram.

15 PageFigure 3 : 16 Page

 15316 Figure 3: The components with sign of the HH Survival copula function for a sample of dimension 80.

Figure 4 : 17 Page

 417 Figure 4: The HH Survival copula of a sample clusterized in 5 groups for an increasing dimension of the sample.

4 candidate set( 1 : 21 Page6

 4121 (end-1)), ... 5 max count(1:(end-1)),fixed count-i); plist = [plist;[temp,repmat(i,size(temp,1),1)]]; 7 endA.2 Testing compatibility with the groupsMoving beyond the purely combinatorial framework, we now select the c.ds consistent with the groups' cardinalities generated by K-means.1 % This function tests the compatibility of the c.c.ds in MPerm 2 % with the groups' cardinalities.

3 4

 3 function[SizeThreshold]= SizeThreshold(MPerm,p) 5 6

[

  following the row's order7 , the corresponding 0-1 matrix is:Then we must substitute the elements of the matrix p 4,3[START_REF] Cherubini | Bivariate Option Pricing with Copulas[END_REF] with the coordinates requested to compute our object, following the methods suggested by Proposition 5 and the following Corollaries 7 and 9. If we want to compute the copula's volume, we must substitute the element of p 4,3[START_REF] Cherubini | Bivariate Option Pricing with Copulas[END_REF] with the corresponding element of the vectors defining the box. For example if the box is S = [col(0.3 × I), col(0.5 × I)] where I is the identity matrix in R 3 , the matrix p 4,3 (2) becomes:c(p 4,3 (2)) =

11 12

 11 function[HCM,pfactor]= claytonLevel2(numones,u,v,sizecluster...

13

 13

15 16 % 20 21

 1620 Given a row corresponding to a c.c.d., this function creates 17 % a matrix with a number of rows equal to the number of groups 18 % and a number of columns equal to the max cardinality of 19 % groups (it's p ik(j)) Amat=mmv(numones,u,v,sizecluster); 22 for j=1:numclusters 23 for i=1:sizecluster(1,j) 24 25 % The function G Clayton computes the copula's generators 26 % for each element of matrix Amat and for the parameter of 27 % the corresponding group.

2829

  Cmat(j,i)=G Clayton(Amat(j,i),alpha(1,j)); 30 end 31 end 32 Cmat=[Cmat(:,:)]; 33 34 % Here we sum over the rows of matrix Cmat till the cardinality 35 % of the corresponding group.

  Finally we compute the copula function at the first level, i.e. we aggregate the second level copula functions with an Archimedean copula with parameter equal to the corresponding EQDP, accounting for the position of the element into the groups.

  ://mc.manuscriptcentral.com/gcom E-mail: ijcm@informa.comInternational Journal of Computer Mathematics

  max{ms,s=1,...,k} whose (s, w)-th element corresponds to the w-th position in the s-th group 1 , where D c (i, k) denotes the number of the ordered compatible combinatorial distributions (o.c.c.ds in the following) 2 of the integer i into k groups, p i,k (j) is a k × max(m s , s = 1, ..., k) matrix of zeros and ones corresponding to the j-th o.c.c.d. of i ones into k groups of dimensions m

s , s = 1, ..., k, c(p i,k (j)) is a k × max(m s , s = 1, ..., k) matrix such that c s,w = ṽs,w if p s,w = 0 and c s,w = ũs,w if p s,w = 1, where ũs,w , ṽs,w and c s,w denote the (s, w)-th element of the corresponding matrix and HC(c(p i,k (j))) is the Hierarchical copula computed for the j-th o.c.c.d..

Table 1 :

 1 K-means HH Survival copula vs CR Survival for several number of groups. Clayton copula with α within = 0.1 and α between = 0.01 computed on a heterogeneous sample of dimension 20.

	Number of groups K-means HH Survival Error
	2	0.9588	0.0039
	4	0.9565	0.0016
	5	0.9560	0.0011
	10	0.9553	0.0004

18 Page 18 of 27 URL

	International Journal of Computer Mathematics
	Number of groups K-means HH Survival/white random Error
	2	0.9587	0.0038
	4	0.9566	0.0017
	5	0.9563	0.0014
	F o
	r
	P
	e e r
	R e v i e w
		O n
		l y

: http://mc.manuscriptcentral.com/gcom E-mail: ijcm@informa.com

Table 2 :

 2 K-means HH Survival copula vs CR Survival for several number of groups. Clayton copula with α within = 0.1 and α between = 0.01 computed on a white randomized heterogeneous sample of dimension 20.

	Number of groups K-means HH Survival/colored random Error
	2	0.9274	0.0275
	4	0.9237	0.0312
	5	0.9229	0.032

Table 3 :

 3 K-means HH Survival copula vs CR Survival for several number of groups. Clayton copula with α within = 0.1 and α between = 0.01 computed on a colored randomized heterogeneous sample of dimension 20.

Table 4 :

 4 SOM HH Survival copula vs CR Survival in several scenarios. Clayton copula with α within = 0.1 and α between = 0.01.

	and Hondela

URL: http://mc.manuscriptcentral.com/gcom E-mail: ijcm@informa.com

Table 6 :

 6 CPU computation time in seconds. K-means HH Probability Distribution of the left queue (till 10 percent) of a variable of dimension n = 60 (2007), larger self-organizing maps rearrange data in a way that is fundamentally topological in character.

  Table 6 reports the computing time for the first 10 percent of the probability distribution of a variable of dimension n = 60 for different numbers of groups. We can observe that the number of groups is the critical variable; the CPU time increases exponentially with the number of group till 15, for the next divisor, i.e. 20, problems of memory stack become unmanageable.

	International Journal of Computer Mathematics	Page 20 of 27
	5 Conclusion
	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y
	20

URL: http://mc.manuscriptcentral.com/gcom E-mail: ijcm@informa.com

The coordinates of the box S are computed as u = col(ũ ), v = col(ṽ ), where the function col : [0, 1] max{ms,s=1,...,k} × [0, 1] k → [0, 1] n sets orderly all the columns of a matrix into the same column.

We refer to D(i, k) as the number of the ways in which one can distribute the integer i into k groups taking into account the order into the groups (o.c.ds for short). If we also take into account the cardinalities of the groups, we refer to D c (i, k) as the number of the ordered compatible (for groups' cardinalities) combinatorial distributions.

We refer to D(i, k) as the number of the ways in which one can distribute the integer i into k groups without taking into account the order of the groups (c.ds for short). If we take into account the groups' cardinalities, then we refer to Dc (i, k) as the number of the compatible (for cardinalities' groups) combinatorial distributions.

The rûle which explains how the coordinates are generated, is presented in the following Corollary.

Acknowledgements: The authors would like to thanks our colleagues and friends Umberto Cherubini and Sabrina Mulinacci for useful comments and support and two anonymous referees which allows to improve this work.


Helena Sofia Rodrigues
email: sofiarodrigues@esce.ipvc.pt

Viana Sofia

Do

M Monteiro

Teresa

Delfim F M ; Torres
email: delfim@ua.pt

M Teresa

T Monteiro

Alan Zinober
email: a.zinober@sheffield.ac.uk

Dengue disease, basic reproduction number and control

Keywords: Dengue, basic reproduction number, stability, Cape Verde, control Dengue, basic reproduction number, stability, Cape Verde, control. 2010 Mathematics Subject Classification: 92B05, 93C95, 93D20

Dengue is one of the major international public health concerns. Although progress is underway, developing a vaccine against the disease is challenging. Thus the main approach to fight the disease is vector control. A model for the transmission of Dengue disease is presented. It consists of eight mutually-exclusive compartments representing the human and vector dynamics. It also includes a control parameter (insecticide) in order to fight the mosquito. The model presents three possible equilibria: two disease-free equilibria (DFE) and another endemic equilibrium. It has been proved that a DFE is locally asymptotically stable, whenever a certain epidemiological threshold, known as the basic reproduction number, is less than one. We show that if we apply a minimum level of insecticide, it is possible to

Introduction

The first recognized Dengue epidemic occurred almost simultaneously in Asia, Africa and North America in the 1780s, shortly after the identification and naming of the disease in 1779. It has spread especially in the tropical and subtropical regions around the world, and nowadays is a disease widely found in urban and semi-urban areas.

According to the World Health Organization, 50 to 100 million Dengue infections occur yearly, including 500,000 Dengue Haemorrhagic Fever cases and 22000 deaths, mostly among children [27]. Growing awareness of global climate change has stimulated several assessments of its likely effects on vector-borne disease as well as on health outcomes. Some of these studies have indicated that countries with a mild climate, such as in the Mediterranean, are at risk due to future climate conditions that may be favourable to this kind of disease [START_REF] Hopp | Global-scale relationships between climate and the dengue fever vector, Aedes Aegypti[END_REF]. This risk may be aggravated further due to the volume of international tourism and trade that this region experiences [START_REF] Semenza | Climate change and infectious diseases in europe[END_REF][START_REF]Mosquito ameaça populações[END_REF]. Travellers play an essential role in the global epidemiology. They act as viremic travellers, carrying the disease into areas with mosquitos that can transmit the infection. This is particularly true with respect to the reality in the archipelago of Cape Verde, where it is believed that the insects responsible for the outbreak that occurred in 2009 came from Brazil, transported by means of the frequent air transport [START_REF]Mosquito transmissor da dengue chega a Cabo Verde[END_REF].

There are two forms of Dengue: Dengue Fever and Dengue Haemorrhagic Fever. The first one is characterized by a sudden fever without respiratory symptoms, accompanied by intense headaches and lasts between three and seven days. The second has the previous symptoms but additionally also nausea, vomiting and fainting due to low blood pressure, and can lead to death in two or three days [START_REF] Derouich | A model of dengue fever[END_REF]. There are four distinct, but closely related, viruses that cause Dengue. Recovery from infection by one virus provides lifelong immunity against that virus but confers only partial and transient protection against subsequent infection by the other three viruses. There is good evidence that sequential infection increases the risk of developing Dengue Haemorrhagic Fever (DHF).

Activities, such as triage and management, are critical in determining the clinical outcome of Dengue. A rapid and efficient front-line response not only reduces the number of unnecessary hospital admissions but also saves lives.

There is no vaccine to protect against Dengue. With four closely related viruses that can cause the disease, a vaccine would need to immunize against all four types to be effective. There is limited understanding of how the disease typically behaves and how the virus interacts with the immune system. Another challenge is that some studies show that some secondary Dengue infection can leave to DHF, and theoretically a vaccine could be a potentially cause of severe disease if a solid immunity is not established against the four serotypes. Research to develop a vaccine is ongoing and the incentives to study the mechanism of protective immunity are gaining more support, now that the number of outbreaks around the world is increasing [START_REF]Dengue: guidelines for diagnosis, treatment, prevention and control[END_REF].

The spread of Dengue is attributed to expanding geographic distribution of the four Dengue viruses and their mosquito vectors, the most important of which is the predominantly urban species Aedes aegypti ; some notes about it are presented in Section 2. A mathematical model of the interaction between human and mosquito populations is presented in Section 3. Section 4 is concerned with the basic reproduction number, the equilibria of the epidemiological model and their stability. In Section 5 previous results are applied to a case study. Finally, concluding notes are given in Section 6.

Biological notes on Aedes aegypti

The mosquito Aedes aegypti is a tropical and subtropical specie widely distributed around the world, mostly between latitudes 35 cle of a mosquito has four distinct stages: egg, larva, pupa and adult. In the case of Aedes aegypti the first three stages take place in or near water whilst air is the medium for the adult stage [START_REF] Otero | A stochastic spatial dynamical model for Aedes aegypti[END_REF]. The eggs of Aedes aegypti can resist desiccation and low temperatures for up to one year. Although the hatching of mature eggs may occur spontaneously at any time, this is greatly stimulated by flooding. The larva moults four times in a period of a few days which culminates in the pupal stage. The pupal stage lasts from one day to a few weeks, depending on the temperature. At the end of this stage the adult emerges from the pupal skin [START_REF] Christophers | Aedes aegypti, the Yellow fever Mosquito: Its Life History, Bionomics and Structure[END_REF]. Studies suggest that most female mosquitoes may spend their lifetime in or around the houses where they emerge as adults. This means that people, rather than mosquitoes, rapidly move the virus within and between communities.

The adult stage of the mosquito is considered to last an average of eleven days in the urban environment. Dengue is spread only by adult females, that require a blood meal for the development of eggs; male mosquitoes feed on nectar and other sources of sugar. In this process the female acquires the virus while feeding on the blood of an infected person. After virus incubation for eight to ten days, an infected mosquito is capable, during probing and blood feeding, of transmitting the virus for the rest of its life. Aedes aegypti is one of the most efficient vectors for arboviruses because it is highly anthropophilic, frequently bites several times before complete oogenesis, and thrives in close proximity to humans [START_REF]Dengue: guidelines for diagnosis, treatment, prevention and control[END_REF].

It is very difficult to control or eliminate the Aedes aegypti mosquito because it adapts to the environment and becomes resistant to natural phenomena, e.g., droughts, or human interventions and control measures. Vector control is a key for combating mosquito-borne diseases and the major tool available for tackling the transmission of Dengue, a disease for which there is no vaccine nor prophylaxis. There are two primary methods for primary prevention: larval control and adult mosquito control, depending on the intended target [START_REF]Bioecologia do Aedes Aegypti. Biológico[END_REF]. In urban areas Aedes mosquitos breed on water collections in artificial containers such as cans, plastic cups, used tires, broken bottles, flower pots, etc. Proper solid waste disposal and improved water storage practices, including covering containers to prevent access by egg-laying fe-

The mathematical model

The mathematical model is based on [START_REF] Dumont | On a temporal model for the Chikungunya disease: modeling, theory and numerics[END_REF][START_REF] Dumont | Vector Control for the Chikungunya Disease[END_REF]. The novelty in this paper is the presence of the control parameter related to adult mosquito insecticide.

The notation used in our mathematical model includes four epidemiological states for humans: It is assumed that the total human population (N h ) is constant, so,

S h (t)
N h = S h + E h + I h + R h .
There are also four other state variables related to the female mosquitoes (the male mosquitoes are not considered in this study because they do not bite humans and consequently they do not influence the dynamics of the disease):

A m (t) aquatic phase (that includes the egg, larva and pupa stages) S m (t) susceptible (mosquitoes that are able to contract the disease) E m (t) exposed (mosquitoes that are infected but are not yet able to transmit to humans) I m (t) infected (mosquitoes capable of transmitting the disease to humans)

In order to analyze the effects of campaigns to fight the mosquito, there is also a control variable: c(t) level of insecticide campaigns Some assumptions in this model: • there is no immigration of infected individuals into the human population;

• the population is homogeneous, which means that every individual of a compartment is homogenously mixed with the other individuals;

• the coefficient of transmission of the disease is fixed and does not vary seasonally;

• both human and mosquitoes are assumed to be born susceptible; there is no natural protection;

• for the mosquito there is no resistant phase, due to its short lifetime.

The parameters of the model are:

N h total
                   dS h dt (t) = µ h N h -(Bβ mh I m N h + µ h )S h dE h dt (t) = Bβ mh Im N h S h -(ν h + µ h )E h dI h dt (t) = ν h E h -(η h + µ h )I h dR h dt (t) = η h I h -µ h R h (1)
and vector population

                   dA m dt (t) = µ b (1 - A m K )(S m + E m + I m ) -(η A + µ A )A m dS m dt (t) = -(Bβ hm I h N h + µ m )S m + η A A m -cS m dE m dt (t) = Bβ hm I h N h S m -(µ m + η m )E m -cE m dI m dt (t) = η m E m -µ m I m -cI m (2)
with the initial conditions

S h (0) = S h0 , E h (0) = E h0 , I h (0) = I h0 , R h (0) = R h0 , A m (0) = A m0 , S m (0) = S m0 , E m (0) = E m0 , I m (0) = I m0 . (3)
Notice that the equation related to the aquatic phase does not have the control variable c, because the adulticide does not produce effects in this stage of the life of the mosquito. To combat the larval phase it would be necessary to use larvicide. This treatment should be long-lasting and have World Health Organization clearance for use in drinking water. As we want to study only a short period of time, this type of treatment has not been considered here.

With the condition

S h + E h + I h + R h = N h , one can, in the example given, use R h = N h -S h -E h -I h and
Ω = {(S h , E h , I h , A m , S m , E m , I m ) ∈ R 7 + : S h + E h + I h ≤ N h , A m ≤ kN h , S m + E m + I m ≤ mN h }
be the region of biological interest, that is positively invariant under the flow induced by the differential system (1)-( 2) (see Appendix A).

Theorem 1. Let Ω be defined as above. Consider also

M = -(c(η A + µ A ) + µ A µ m + η A (-µ b + µ m )).
The system (1)-( 2) admits at most two disease free equilibrium points:

• if M ≤ 0, there is a Disease-Free Equilibrium (DFE), called Trivial Equilibrium, E * 1 = (N h , 0, 0, 0, 0, 0, 0);

• if M > 0, there is a Biologically Realistic Disease-Free Equilibrium (BRDFE), E * 2 = N h , 0, 0, kN h M η A µ b , kN h M µ b µm , 0, 0 . Proof. See Appendix B.
Remark 1. The condition M > 0 is equivalent, by algebraic manipulation, to the condition

(η A + µ A )(µ m + c) µ b η A > 1
which corresponds to the basic offspring number for mosquitos. Thus, if M < 0, then the mosquito population will collapse and the only equilibrium for the whole system is the trivial equilibrium. If M ≥ 0, then the mosquito population is sustainable.

It is necessary to determine the basic reproduction number of the disease, R 0 . This number is very important from the epidemiologistic point of view. It represents the expected number of secondary cases produced in a completed susceptible population, by a typical infected individual during its entire period of infectiousness [START_REF] Heffernan | Perpective on the basic reproductive ratio[END_REF][START_REF] Hethcote | The mathematics of infectious diseases[END_REF]. Following [START_REF] Van Den Driessche | Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission[END_REF], we prove:

(2) is R 2 0 = B 2 S h0 S m0 β hm β mh η m ν h N 2 h (η h + µ h )(c + µ m )(c + η m + µ m )(µ h + ν h ) .
The equilibrium point BRDFE is locally asymptotically stable if R 0 < 1 and

unstable if R 0 > 1.
Proof. See Appendix C.

Remark 2. In our model we have two different populations (human and vectors), so the expected basic reproduction number reflects the infection humanvector and also vector-human, that is,

R 2 0 = R hm × R mh . The term Bβ hm S m0 N h
represents the product between the transmission probability of the disease from humans to vectors and the number of susceptible mosquitos per human; 1 η h +µ h is related to the human's viremic period; and ηm c+ηm+µm represents the proportion of mosquitoes that survive the incubation period. Analogously, the term Bβ mh S h0 N h is related to the transmission probability of the disease between mosquitos and human, in a susceptible population; 1 (c+µm) represents the lifespan of an adult mosquito; and ν h (µ h +ν h ) is the proportion of humans that survive the incubation period.

When R 0 < 1, each infected individual produces, on average, less than one new infected individual, and therefore it is predictable that the infection will be cleared from the population. If R 0 > 1, the disease is able to invade the susceptible population.

* 3 = (S * h , E * h , I * h , A * m , S * m , E * m , I * m )
, where,

S * h = N h - (µ h + ν h )(µ h + η h ) µ h ν h I * h , E * h = µ h + η h ν h I * h , I * h = ξ χ , ξ = N h µ h -B 2 kβ hm β mh ν h η m M + µ b µ 2 m (η m + µ m )(µ h + ν h )(µ h + η h ) + c 2 µ b (η h + µ h )(µ h + ν h )(c + η m + 3µ m ) + cµ b µ m (µ h + ν h )(µ h (3µ m + 2) + η h (2η m + 3µ m )) , χ = Bβ hm (η h + µ h ) -µ b µ h (c + µ m )(c + η m + µ m ) -Bkβ mh η m M (µ h + ν h ), A * m = M η A µ b kN h , S * m = kNh 2 M µ b (cN h + BI * h β hm + N h µ m ) , E * m = µ m + c η m I * m , I * m = BI * h kN h β hm η m M µ b (c + µ m )(c + η m + µ m )(cN h + BI * h β hm + N h µ m ) . Proof. See Appendix B.
From a biological point of view, it is desirable that humans and mosquitoes coexist without the disease reaching a level of endemicity. We claim that proper use of the control c, can result in the basic reproduction number remaining below unity and, therefore, making BRDFE stable.

In order to make effective use of achievable insecticide control, and simultaneously to explain more easily to the competent authorities its effectiveness, we assume that c is constant. The goal is to find c such that R 2 0 < 1. [START_REF] Thomé | Optimal control of Aedes aegypti mosquitoes by the sterile insect technique and insecticide[END_REF][START_REF] Yang | Assessing the effects of temperature on dengue transmission[END_REF]. The simulations were carried out using the following values:

F

N h = 480000, B = 1, β mh = 0.375, β hm = 0.375, µ h = 1/(71 × 365), η h = 1/3, µ m = 1/11, µ b = 6, µ A = 1/4, η A = 0.08, η m = 1/11, ν h = 1/4, , m = 6, k = 3, K = kN h .
The initial conditions for the problem were:

S h0 = N h -E h -I h , E h0 = 216, I h0 = 434, R h0 = 0, A m0 = kN h , S m0 = mN h .
Considering nonexistence of control, the basic reproduction number for this outbreak in Cape Verde is approximately R 0 = 4.85, which is in agreement to other studies of Dengue in other countries [START_REF] Nishiura | Mathematical and Statistical Analyses of the spread of Dengue[END_REF]. The control c affects the basic reproduction number, and our aim is to find a control that puts R 0 less than one. Proposition 1. Let us consider the parameters listed above and consider c as a constant. Then R 0 < 1 if and only if c > 0.0837. The computational investigations were carried out using c = 0.084, which means that the insecticide is continuously applied during a 84 day period.

The software used was Scilab [START_REF] Campbell | Modeling and simulation in Scilab/Scicos[END_REF]. It is an open source, cross-platform numerical computational package and a high-level, numerically oriented programming language. For our problem we used the routine ode to solve the set of differential equations. By default, ode uses the lsoda solver of package ODEPACK. It automatically selects between the nonstiff predictor- corrector Adams method and the stiff Backward Differentiation Formula (BDF) method. It uses the nonstiff method initially and dynamically monitors data in order to decide which method to use. The graphics were also obtained using this software, using the command plot.

Figures 1 and2 show the curves related to human population, with and without control, respectively. The number of infected persons, even with small control, is much less than without any insecticide campaign. When the control is applied, the number of infected mosquitoes is close to zero. Note that the intention is not to eradicate the mosquitoes but instead the number of infected mosquitoes.

Conclusions and future work

In this paper a model based on two populations, humans and mosquitoes, with insecticide control has been presented. It has been shown that with a steady insecticide campaign it is possible to reduce the number of infected humans and mosquitoes, and can prevent an outbreak that could transform an epidemiological episode to an endemic disease. It has been proved algebraically that, if a constant minimum level of a control is applied (c = 0.084), it is possible to maintain the basic reproduction number below unity, guaranteeing the BRDFE. This value is corroborated in another numerical study [START_REF] Rodrigues | Insecticide control in a dengue epidemics model[END_REF].

In this work we considered a constant control. In future work, using a theoretical approach [START_REF] Rodrigues | Dynamics of dengue epidemics when using optimal control[END_REF], we intend to find the best function c(t), using optimal control theory. Instead of finding a constant control it will be possible to study other types of control, such as piecewise constant or even continuous but not constant. Numerical methods can be used to solve the model [START_REF] Simos | A new modified Runge-Kutta-Nystrom method with phase-lag of order infinity for the numerical solution of the Schrodinger equation and related problems[END_REF][START_REF] Simos | A modified Runge-Kutta method with phase-lag of order infinity for the numerical solution of the Schrodinger equation and related problems[END_REF][START_REF] Vigo-Aguiar | A parallel boundary value technique for singularly perturbed two-point boundary value problems[END_REF]. Additionally we could consider another strategy, a more practical one: due to logistics and health reasons, it may be more convenient to apply insecticide periodically and at some specific hours at night. The rapid increase in mosquito resistance to several chemical insecticides, and the damage caused by these to the environment, has resulted in the search for new control alternatives, such as the use of biological agents. Among the alternatives available, the use of Bacillus thuringiensis israelensis (Bti) has been adopted by several countries [START_REF]Diretrizes nacionais para prevenção e controle de epidemias de dengue[END_REF]. Laboratory testing shows that Bti has a high larvicide property and its mechanism of action is based on the production of an endotoxin protein that, when ingested by the larvae, causes death. So, in future work, we will also take into account the use of

-µm -ηm µ b µ b µ b 0 0 0 η A -Bβ hm I h N h -µm -c 0 0 0 0 0 0 Bβ hm I h N h -µm -ηm -c 0 0 0 0 0 0 ηm -µm -c 1 C C C C C C C C C C C C C A
, and F = (µ h N h , 0, 0, 0, 0, 0, 0) T . As M(X) has all off-diagonal entries nonnegative, M(X) is a Metzler matrix.

E * 2 = N h , 0, 0, kN hM η A µ b , kN hM µ b µm , 0, 0 .
This is biologically interesting only if M is greater than 0.

The third solution corresponds to a situation where humans and mosquitoes live together but the disease persists in the two populations. So the disease is not anymore an epidemic episode, but transforms into a endemic one. With some algebraic manipulations we obtain the following point:

E * 3 = (S * h , E * h , I * h , A * m , S * m , E * m , I * m ) where, S * h = N h - (µ h + ν h )(µ h + η h ) µ h ν h I * h , E * h = µ h + η h ν h I * h , I * h = N h µ h (-B 2 kβ hm β mh ν h η m M + µ b µ 2 m (η m + µ m )(µ h + ν h )(µ h + η h ) + c 2 µ b (η h + µ h )(µ h + ν h )(c + η m + 3µ m ) + cµ b µ m (µ h + ν h )(µ h (3µ m + 2) + η h (2η m + 3µ m )))/(Bβ hm (η h +µ h )(-µ b µ h (c+µ m )(c+η m +µ m )-Bkβ mh η m M)(µ h +ν h )); A * m = M η A µ b kN h , S * m = kNh 2 M µ b (cN h + BI * h β hm + N h µ m ) , E * m = µ m + c η m I * m , I * m = BI * h kN h β hm η m M µ b (c + µ m )(c + η m + µ m )(cN h + BI * h β hm + N h µ m ) As before, this equilibrium is only biologically interesting if M > 0.
With the Mathematica software we obtained a fourth solution. But some of the components are negative, which means that they do not belong to the Ω set.

C The basic reproduction number: proof of Theorem 2

The basic reproduction number is calculated in a disease free equilibrium. In this case we consider the most realistic one, BRDFE.

In [START_REF] Allen | The basic reproduction number in some discrete-time epidemic models[END_REF][START_REF] Van Den Driessche | Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission[END_REF] they considered the vector x T = (E h , I h , E m , I m ) which corresponds to the components related to the progression of the disease.

The subsystem used was:

This subsystem can be written as partitioned, dx dt = F (x) -V(x), where

x T = (E h , I h , E m , I m ), F (x) represents the components related to new cases of disease (in this situation in the exposed compartments) and V(x) represents the other components. Thus the subsystem (6) can be rewritten as

F (x) =       Bβ mh Im N h S h 0
Bβ hm I h N h S m 0       and V(x) =       (ν h + µ h )E h -ν h E h + (η h + µ h )I h (µ m + η m + c)E m -η m E m + (µ m + c)I m       .
Let us considerer the Jacobian matrices associated with F and V:

J F (x) =       0 0 0 Bβ mh S h N h 0 0 0 0 0 Bβ hm Sm N h 0 0 0 0 0 0       and J V(x) =       ν h + µ h 0 0 0 -ν h η h + µ h 0 0 0 0 µ m + η m + c 0 0 0 -η m µ m + c       .
According to [START_REF] Van Den Driessche | Reproduction numbers and sub-threshold endemic equilibria for compartmental models of disease transmission[END_REF] the basic reproduction number is R 0 = ρ(J F (x 0 ) J V -1 (x 0 ) ), where x 0 is a disease free equilibrium (BRDFE) and ρ(A) defines the spectral radius of a matrix A. Using Mathematica

  mosquitoes, are methods that are encouraged through community-based programmes [4]. Active monitoring and surveillance of the natural mosquito population should accompany control efforts to determine programme effectiveness.

•

  the total human population (N h ) is constant;

Page 8 of 20 URLTheorem 2 .

 202 : http://mc.manuscriptcentral.com/gcom E-mail: ijcm@informa.If M > 0, then the square of the basic reproduction number associated to (1)-

Theorem 3 .

 3 If M > 0 and R 0 > 1, then the system (1)-(2) also admits an

Fig. 1 :Fig. 2 :

 12 Fig. 1: Human compartments using control Fig. 2: Human compartments without control

Page 12 of 20 URLFig. 3 :Fig. 4 :

 2034 Fig. 3: Mosquito compartments using control Fig. 4: Mosquito compartments without control

B 2

 2 kβ hm β mh η m ν h M µ b (η h + µ h )(c + µ m ) 2 (c + η m + µ m )(µ h + ν h )and we obtain the value for the threshold parameter, with M > 0. Page 20 of 20 URL: http://mc.manuscriptcentral.com/gcom E-mail: ijcm@informa.com International Journal of Computer Mathematics

  consider an equivalent system for human population without considering the R h differential equation.

	4 Basic Reproduction number, equilibrium points and
	stability
	Let the set
	F o
	r
	P
	e e r
	R e v i e w
	O n
	l y

  country. A total of 17 224 cases, including six deaths, were reported from 18 of the 22 municipalities in Cape Verde by the end of 2009. The municipality of Praia, on Santiago island, notified the highest number of cases (13000 cases) followed by São Felipe in Fogo Island (3000 cases). We used the data for human population related to Cape Verde [4]. There have been no cases of Dengue in all the West African countries near Cape Verde islands since 2000 ([START_REF]Mosquito transmissor da dengue chega a Cabo Verde[END_REF]). They speak explicitly about mosquitoes coming from Brazil. Also the information that comes from the Ministry of Health in the capital of Cape Verde, Praia, confirms that the insects responsible for Dengue came most probably from Brazil, transported by means of air transport that perform frequent connections between Cape Verde and Brazil, as reported by José Rosa in the Radio of Cape Verde. With respect to Aedes aegypti, we have thus considered data from Brazil

	o
	r
	P
	e e r
	R e v i e w
	O n
	l y

  S h -(ν h + µ h )E h dI h dt (t) = ν h E h -(η h + µ h )I h dEm dt (t) = Bβ hm I h N h S m -(µ m + η m )E m -cE m dIm dt (t) = η m E m -µ m I m -cI m

	Page 19 of 20	International Journal of Computer Mathematics
		   	dE h dt (t) = Bβ mh	Im N h
		
		
		
		
		
		 
		F o
		r
		P
			e e r
			R e v i e w
				O n
				l y

URL: http://mc.manuscriptcentral.com/gcom E-mail: ijcm@informa.com

Basic Reproduction number, equilibrium points and stability

URL: http://mc.manuscriptcentral.com/gcom E-mail: ijcm@informa.comInternational Journal of Computer Mathematics

Dengue in Cape Verde

Conclusions and future work 12

Conclusions and future work 13

Acknowledgements

Work partially supported by Portuguese Foundation for Science and Technology (FCT) through the PhD Grant SFRH/BD/33384/2008 (Rodrigues) and the R&D units Algoritmi (Monteiro) and CIDMA (Torres). The authors are very grateful to three referees, for valuable remarks and comments, which has contributed significantly to the quality of the paper.

A Dynamical properties of Metzler systems

System (1)-( 2) can be rewritten in the following way:

where

Using the fact that F ≥ 0, the system (4) is positively invariant in R 7 + [START_REF] Abate | Box invariance in biologically-inspired dynamical systems[END_REF], which means that any trajectory of the system starting from an initial state in the positive orthant R 7 + remains forever in R 7 + .

B Equilibrium points

The equilibrium points are reached when the following equations hold:

Using the Mathematica software to solve the system (5), we obtained four solutions.

The first one is known as the Trivial Equilibrium, since the mosquitoes do not exist, so there is no disease:

In the second one, mosquitoes and humans interact, but there is only one outbreak of the disease, i.e., over time the disease goes away without it being necessary to kill all the mosquitoes. We have called this equilibrium point a Biologically Realistic Disease-Free Equilibrium (BRDFE), since it is a more reasonable situation to find in nature than the previous one:

which is equivalent to