Nicolas Debarsy
email: nicolas.debarsy@cnrs.fr

James P Lesage
email: james.lesage@txstate.edu

Using convex combinations of spatial weights in spatial autoregressive models

Keywords: Markov Chain Monte Carlo estimation, Taylor series approximation, log-marginal likelihood, multiple weight matrices

Spatial regression models rely on simultaneous autoregressive processes that model spatial or cross-sectional dependence between cross-sectional observations using a weight matrix. A criticism of applied spatial regression methods is that reliance on geographic proximity of observations to form the weight matrix that specifies the structure of cross-sectional dependence might be unrealistic in some applied modeling situations. In cases where the structure of dependence or connectivity between (cross-sectional) observations arises from non-spatial relationships, spatial weights are theoretically unjustifiable. Some literature addresses the structure of dependence between observations by introducing geographic proximity as well as other types of non-spatial proximity, resulting in a model that utilizes multiple weight matrices. Each set of weights reflect a different type of dependence specified using linear combinations of observations defined by alternative characteristics. The multiple weight matrix approach results in a simultaneous autoregressive process that poses a number of challenges for parameter estimation and model interpretation.

We focus on literature that relies on a single connectivity matrix constructed from a convex combination of multiple matrices, each of which reflects a different type of dependence or interaction structure. The advantage of this approach is that the resulting simultaneous autoregressive process is amenable to conventional spatial regression estimation algorithms as well as methods developed for interpretation of estimates from these models. Estimates of the scalar parameters used to form the convex combination of the weight matrices can be used to produce an inference regarding the relative importance of each type of dependence.

Introduction

Spatial regression models typically rely on spatial proximity or Euclidean distance between observations to specify the structure of simultaneous dependence between observations. For example, neighboring regions that have common borders with each region are defined as those on which each observation is dependent. Alternatively, dependent observations for each region can be defined as a function of inverse distance to other observations, perhaps up to some cut-off distance. However, this is often viewed as restrictive and inappropriate in applications where spatial proximity or distance is theoretically unjustifiable as a basis for the underlying dependence structure of the simultaneous autoregressive process (Corrado and Fingleton, 2012). This criticism is most relevant when spatial regression models are used in broader contexts to model spatial phenomena involving interregional flows of goods or population, knowledge transmission between regions, or non-spatial phenomena such as student peer groups or social networks. In these scenarios, the spatial location of observations may not be theoretically appealing as a basis for the underlying structure of dependence between observations.

A simultaneous autoregressive process models variation in an N × 1 variable vector y using other observations in the vector y. An N ×N matrix W is used to define specific observations on which each observation is dependent.

The matrix W has zeros on the diagonal to prevent an observation from being dependent on itself, and is typically row-normalized to have row-sums of one. This type of normalization aids with interpretation of the model, allowing us to compare the magnitude of scalar dependence parameter estimates arising from models specified with alternative definitions of the matrix W . A simultaneous autoregressive process specifies a relationship between the vector of observations in y and the N × 1 vector W y, resulting from the product of the N × N matrix with the N × 1 vector y. The relationship takes the form: y = ρW y + ε, where the scalar parameter ρ defines the strength of dependence between the N × 1 vector y and N × 1 vector W y.

In the case of row-normalized W , the relationship is between observations in the vector y and a linear combination of observations on which each observation is dependent. The N × 1 vector ε is a stochastic disturbance term, typically assumed to be normally distributed with zero mean and a scalar variance, σ 2 . We focus on maximum likelihood-based estimation procedures where normality is important, but alternative estimation methods, such as the generalized method of moments do not require normally distributed disturbances. In situations where the dependence structure in the matrix W is defined based on spatial aspects of a sample of observations located in space, we can label this a spatial autoregressive process.

One can also define a simultaneous autoregressive process involving the model disturbances, for example, y = µι N +u, u = λW u+ε where the scalar µ = E(y), ι N is an N × 1 vector of ones, and u is a N × 1 set of disturbances. Here, the disturbances exhibit dependence on disturbances from other observations identified by the matrix W , with the strength of dependence defined by the scalar parameter λ. Generally, simultaneous autoregressive processes model outcomes (y) or disturbances (u) as related to a linear combination of other observations defined by the dependence structure. This could involve a scenario where each observation is possibly dependent on a different number of observations in the dependence set, with possibly different weights assigned to each of these sets of dependent observations.

Simultaneous autoregressive processes can be used to extend ordinary regression models, producing specifications such as: y = ρW y + Xβ + ε, which has been labeled a spatial autoregressive model (SAR), where X is an N × K matrix of explanatory variables and ε an N × 1 vector of normally distributed zero mean, constant variance disturbances. Of course, this type of model specification can be viewed more generally as a simultaneous autoregressive model when the application involves non-spatial observations. Another spatial regression specification is the spatial error model (SEM), y = Xβ + u, u = λW u + ε, where the N × 1 vector of disturbances u follow a spatial autoregressive process. Matrix-vector products such as W y, W u are referred to as spatial lags, and we can also use a matrix-matrix product W X to create spatial lags of the K explanatory variables of the model. In the case of non-spatial W matrices, we label these simultaneous lags, where W y, W u, W X simply reflect linear combinations of observations from the dependence set of observations. These reflect a linear combination of explanatory variables from dependent observations specified by the matrix W . These additional explanatory variables can enter the spatial regression model: y = ρW y + Xβ + W Xθ + ε, which has been labeled a spatial Durbin model (SDM). In a spatial setting where the dependence structure reflects nearby or neighboring regions, this model indicates that explanatory variable values of neighboring regions may help explain variation in own-region outcomes reflected in the y-vector. More generally, in situations where the matrix W reflects non-spatial dependence relations, the set of explanatory variables defined by W X have been labeled contextual variables, reflecting the notion that a linear combination of explanatory variables from observations in the dependence set describe a context that may impact outcomes in the vector y.

One important issue in the spatial regression literature concerns the appropriateness of the connectivity matrices W to define the underlying simultaneous autoregressive process. An intuitive motivation for reliance on non-spatial weight matrices arises naturally when one considers that spatial dependence relies on the notion of Euclidean distance or spatial proximity to define observations in the dependence set. Euclidean distance can be generalized to define a dependence set of observations based on other metrics. For example, [START_REF] Pace | A method for spatial-temporal forecasting with an application to real estate prices[END_REF] propose a model for prices of homes sold that occur at irregular points in space and time, generalizing distance to include relative locations in time. Related work by Pace et al. (2002) relied on generalized distances that considered the number of bedrooms and bathrooms (of nearby homes) to specify the structure of dependence between the selling prices of homes, with the motivation that appraisers determine selling price estimates based on homes comparable in these two metrics.

Once we open the door to non-spatial metrics as a way to specify dependence between cross-sectional observations, a host of issues arise, which are discussed in LeSage and Pace (2011), and Debarsy and LeSage (2018a). Blankmeyer et al (2011, p.94) point out that "a single weight matrix, based on a multivariate similarity criterion (generalized distance) requires a norm to prevent scale differences from influencing the weight placed on the various measures of similarity. (This is unlike the case of spatial proximity where Euclidian distance provides a natural scaling)".

Multiple connectivity matrices

An approach that avoids relative scaling issues involves extensions of the SAR specification to include more than one spatial lag (see, among others [START_REF] Badinger | Estimation of higher-order spatial autoregressive cross-section models with heteroscedastic disturbances[END_REF][START_REF] Lacombe | Does econometric methodology matter? An analysis of public policy using spatial econometric techniques[END_REF][START_REF] Lee | Efficient GMM estimation of high order spatial autoregressive models with autoregressive disturbances[END_REF][START_REF] Elhorst | On model specification and parameter space definitions in higher order spatial econometric models[END_REF].

These models extend our SAR specification to include multiple spatial or simultaneous lags, each one based on a different connectivity matrix, in an effort to capture different types of dependence between cross-sectional observations using the simultaneous autoregressive specification. The additional simultaneous lags of the dependent variable vector y have been labeled 'higher-order' terms, with this type of extended SAR model shown in Eq.

(1), where L simultaneous lag terms are introduced.

y = ( L ∑ ℓ=1 ρ ℓ W ℓ ) y + Xβ + ε (1)
In Eq. ( 1), y is an N × 1 vector of dependent variable outcomes, X is an exogenous N × K explanatory variables matrix, with β the associated K × 1 parameters, ρ ℓ , ℓ = 1, • • • , L are scalar dependence parameters measuring the strength of each type of dependence modeled by the

N × N dependence matrices W ℓ , ℓ = 1, • • • , L.

Convex combinations of multiple weight matrices

LeSage and Pace (2011) point notably to a number of estimation issues that arise for higher-order models of the type in Eq. (1). In addition, [START_REF] Elhorst | On model specification and parameter space definitions in higher order spatial econometric models[END_REF] point to estimation complications that arise relating to the parameter space for the dependence parameters ρ ℓ , ℓ = 1, • • • , L. The parameter space for ρ 1 depends on values taken by the other parameters ρ j , j = 2, • • • , L, and similarly, the parameter space for ρ 2 depends on values taken by ρ 1 , ρ j , j = 3, . . . ,L, and so on. They provide a solution based on a trigonometric approach to this problem for the case of two W -matrices, that would allow maximum likelihood estimation to be carried out. However, they note that extension of their approach to models containing more than two weight matrices becomes increasingly difficult in situations where the number of weight matrices increases.

An alternative approach that relies on convex combinations of different dependence matrices to form a single weight matrix, was first explored by [START_REF] Pace | Semiparametric maximum likelihood estimates of spatial dependence[END_REF]. They proposed use of a model shown in Eq. ( 2), where a convex combination of multiple N ×N primitive spatial basis matrices B (i) , i = 1 . . . p, each based on a differing sequence of individual doubly stochastic nearest neighbor weight matrices is used to form a single weight matrix, W c = ∑ p i=1 a i B (i) . The model parameters are subject to restrictions set forth in (3) and ( 4). (A doubly stochastic matrix is one whose rows and columns sum to unity.) This leads to a spatial autoregressive model in Eq. ( 2).

y = ρW c y + Xβ + ε (2) W c = p ∑ i=1 a i B (i) 0 ≤ a i ≤ 1, i = 1 . . . p (3) p ∑ i=1 a i = 1 (4)
The model specification in Eq. ( 2) overcomes the interpretation problem noted for the higher-order spatial autoregressive model because the estimates for the parameters ∑ p i=1 a i = 1 can be used to draw inferences regarding the relative importance of the different spatial basis matrices B (i) in the model. It also remedies the problem regarding the parameter space for ρ because W c is a row-normalized weight matrix when individual weight matrices B (i) are row-normalized, which is of course true of the doubly stochastic weights B (i) . We note that a convex combination of row-normalized weight matrices is also row-normalized. This means the parameter space for ρ would obey the standard lower and upper bounds of 1/µ min and 1, where µ min is the smallest eigenvalue of the matrix W c . In typical estimation algorithms for spatial autoregressive models the value of 1/µ min is set to -1, since this ensures that the matrix inverse: (I N -ρW c ) -1 exists, and avoids the need to calculate the minimum eigenvalue of the N × N weight matrix W c . [START_REF] Hazir | The role of R&D collaboration networks on regional innovation performance[END_REF] extend this approach to the more general case of a convex combination of L row-normalized N × N weight matrices based on alternative underlying dependence sets for the N observations reflected by underlying matrices W ℓ , so that

W c (Γ) = ∑ L ℓ=1 γ ℓ W ℓ , with 0 ≤ γ ℓ ≤ 1, ℓ = 1 . . . L, ∑ L ℓ=1 γ ℓ = 1 and Γ = (γ 1 , . . . , γ L ) ′ .
The resulting convex combination of underlying weight structures, W c (Γ), can be used to specify dependence between N observations based on a convex combination of L different types of dependence between observations. The scalar parameters γ ℓ indicate the relative importance assigned to each type of dependence in the cross-sectional dependence scheme. The two sets of constraints imposed reflect the fact that this approach relies on a convex combination. When each W ℓ , ℓ = 1, • • • , L, is row-normalized, then W c (Γ) obeys the conventional row-normalization, avoiding the need to re-normalize the matrix for new values of Γ. We note that the approach based on a convex combination of weight matrices also avoids the issue of scaling for different metrics that arises for generalized measures of distance by casting the problem as one of relative distance/proximity, inherent in conventional spatial regression weight matrices.

To estimate the model parameters ρ, β, Γ, σ 2 , [START_REF] Hazir | The role of R&D collaboration networks on regional innovation performance[END_REF] propose evaluating the likelihood numerous times over a grid of values for the parameters γ ℓ , ℓ = 1 . . . L, and using the vector of parameters (that we denote Γ) associated with the highest likelihood function value as point estimates for these parameters. For example, in the case of two weight matrices where we have values of γ 1 and γ 2 = (1 -γ 1 ), they would define a grid of values for γ 1 = 0, 0.1, 0.2, . . . , 1, and perhaps refine the grid based on an initial estimate from the course grid of (say) 0.2 to be a grid defined over 0.19, 0.20, 0.21, . . . , 0.29. The value of γ1 from this grid of values would be used in the second step of the estimation procedure. The second step would produce estimates for the other parameters ρ, β, σ2 along with measures of dispersion for these parameters by simply fixing the parameters Γ and generating

W c ( Γ) = γ1 W 1 + γ2 W 2 + . . . + γL W L .
Treating W c ( Γ) as fixed, conventional maximum likelihood algorithms are used to estimate ρ, β, σ 2 and associated measures of dispersion. This approach has the virtue of simplicity, since it allows use of conventional maximum likelihood estimation software to estimate the parameters Γ in the first step, and the parameters ρ, β, σ 2 in the second step. An undesirable aspect of this approach is the lack of dispersion estimates for the point estimates of the parameters Γ on which to draw inferences regarding whether individual γ ℓ parameters are significantly different from zero.

Bayesian estimation of convex combinations of weight matrix models

To overcome the problem of inference regarding individual parameters γ ℓ , Debarsy and LeSage (2018a) propose Markov Chain Monte Carlo (MCMC) estimation of the SAR model with convex combination of connectivity matrices. MCMC estimation involves sampling for the complete sequence of conditional distributions of the model parameters. A large sample of draws from the conditional distribution of each parameter is used to construct the joint posterior distribution for the model parameters. In the case of only two weight matrices, Γ = (γ 1 , (1 -γ 1 )) ′ , this involves sampling from the following sequence of conditional distributions: (ρ|β, σ 2 , γ 1 ), (β|ρ, σ 2 , γ 1 ), (σ 2 |ρ, β, γ 1 ), and (γ 1 |ρ, β, σ 2 ). A virtue of MCMC sampling is that the conditional distributions for each parameter given values of all other model parameters typically takes a simple form that is computationally easy to sample from. As an example, the conditional distribution for the K × 1 parameter vector β|ρ, σ 2 , γ 1 takes the form of a K-variate normal distribution with a mean vector and variance-covariance matrix shown in Eqs. ( 5) and ( 6) (when no prior distributions are assigned to the model parameters). When we assign no prior distributions to the model parameters, MCMC sampling serves as simply an estimation method for producing parameter estimates that will be equivalent to those from maximum likelihood estimation. The motivation for use of MCMC estimation in the case of no Bayesian prior information is simply to exploit the computationally simple forms taken by the conditional distributions to make estimation a simpler task than that of optimizing a likelihood function with respect to a large number of parameters. This approach is sometimes labeled simulated maximum likelihood estimation.

β|ρ, σ 2 , γ 1 ∼ N ( β, Σβ ) β = (X ′ X) -1 X ′ [I N -ρW c (γ 1 )]y (5) Σβ = σ 2 (X ′ X) -1 (6)
For the case of only two weight matrices involving the parameter γ 1 (and γ 2 = 1 -γ 1 ), the conditional posterior for γ 1 (given β, ρ, σ 2 ) can be written as in Eq. ( 7), where |A| denote the determinant of the matrix A.

p(γ 1 |ρ, β, σ 2 ) ∝ |A|exp(-1/2σ 2 )(Ay -Xβ) ′ (Ay -Xβ) (7) A = [I N -ρW c (γ 1 )]
A Metropolis-Hastings (M-H) approach can be used to draw samples of the parameter γ 1 (see LeSage and Pace, 2009, chap. 5). This involves comparing the conditional distribution in Eq. ( 7) based on the current value of γ 1 to a proposed value, and accepting or rejecting the proposal according to the M-H rules. The restriction that 0 ≤ γ 1 ≤ 1 is imposed by rejecting values proposed for γ 1 that do not obey the restriction.

A similar conditional distribution arises for the parameter ρ, shown in Eq. ( 8), which is also amenable to M-H sampling. Here again, the restriction -1 < ρ < 1 can be imposed by rejecting proposal values that do not obey the restriction.

p(ρ|γ 1 , β, σ 2 ) ∝ |A|exp(-1/2σ 2 )(Ay -Xβ) ′ (Ay -Xβ) (8) A = [I N -ρW c (γ 1 )]
LeSage and [START_REF] Lesage | Introduction to Spatial Econometrics[END_REF] set forth a random-walk procedure for producing M-H proposal values for the parameter ρ, which can also be used for the parameter γ 1 . Debarsy and LeSage (2018a) actually take a different approach to sampling the parameters ρ, γ 1 , labeled a griddy Gibbs sampler by LeSage and Pace (2009, chap. 5), which should produce equivalent MCMC draws from which the joint posterior estimates for the parameters γ 1 , ρ can be constructed.

Simple extension of the approach set forth above to cases involving more than two weight matrices encounters some difficulties. One is that each parameter γ i must be sampled conditional on all other parameters γ j , j ̸ = i. Another is that the adding up constraint ∑ L ℓ=1 γ ℓ = 1 cannot be imposed when we treat γ j , j ̸ = i as known and fixed. Quite simply a proposal value for γ i requires that we would need to adjust γ j , j ̸ = i to be consistent with the summing up constraint, which rules out treating these as fixed. Of course, this issue of the adding up constraint does not arise in the case of two weight matrices where γ 2 = (1 -γ 1 ), since the constraint holds for all values 0 ≤ γ 1 ≤ 1. Further, the Jacobian of the transformation to be assessed at each pass through the sampler cannot be pre-computed for a grid of values for ρ and Γ. Debarsy and LeSage (2018b) propose an approach that generalizes to specifications involving more than two matrices. Their approach involves rewriting the convex combination SAR model presented in Eq. ( 2) as in Eq. ( 9), which they argue has a number of computational advantages. One advantage is that the parameter vector ω isolates the parameters ρ, γ ℓ , ℓ = 1, . . . ,L from sample data matrices. Separating model parameters from vectors and matrices that contain N sample data observations allows calculation of matrices such as the N × (L + 1) matrix ỹ once, prior to beginning the MCMC sampling loop. More regarding this later.

ỹω = Xβ + ε (9) ỹ = (y, W 1 y, W 2 y, . . . , W L y) ω =        1 -ργ 1 -ργ 2 . . . -ργ L        = ( 1 -ρΓ ) , Γ =      γ 1 γ 2 . . . γ L      , 0 ≤ γ ℓ ≤ 1, ℓ = 1 . . . L
2 Computational issues in estimating convex combination of weights models Debarsy and LeSage (2018b) point to three computational challenges that arise for the convex combination of weight matrices model. One is the need to impose the adding up constraint and the non-negativity requirement for the parameters ∑ L ℓ=1 γ ℓ = 1, which would require use of a constrained optimization algorithm. In the case of MCMC sampling, there is the issue of sampling a single parameter γ i conditional on other parameters γ j , j ̸ = i, while maintaining the adding up constraint and the non-negativity requirement. A second is the need to evaluate a log-determinant term that arises in the likelihood and the conditional distributions for the spatial dependence parameters (ρ, Γ) of the model. The log-determinant depends on both ρ as well as the vector of parameters Γ, so changes in the values of these require calculation of new values for the log-determinant |A| = |I N -ρW c (Γ)| on each parameter draw from the conditional distributions of ρ, Γ. A third problem they note, involves calculating simulated (empirical) measures of dispersion for the partial derivatives ∂y/∂x used to interpret estimates from the model. We discuss the nature of these problems in more detail along with solutions proposed by Debarsy and LeSage (2018b) in the next three sections.

MCMC sampling of the parameters Γ

Debarsy and LeSage (2018b) propose block sampling the entire vector of parameters Γ, using a Metropolis-Hastings (M-H) approach. This means that a vector of all γ ℓ parameters, labeled Γ p is proposed and either accepted or rejected using the M-H rule. If rejected, the current vector of all γ ℓ parameters, labeled Γ c , is retained for the next pass through the MCMC sampler. The virtue of this block sampling approach is that the proposed vector Γ p can be constructed to ensure the adding up constraint and nonnegativity requirements are met.

M-H sampling requires that we evaluate the conditional distribution for the parameters Γ at the current and proposed values for the vectors Γ p , Γ c . Debarsy and LeSage (2018b) show that by integrating out β and σ 2 , we can write the conditional distribution of Γ conditioning only on ρ. As such, the M-H procedure to accept a new block Γ p uses the conditional distribution of γ, assuming ρ is given or known. An acceptance probability is calculated based on evaluating the expression shown in Eq. ( 10) and used to decide if the proposed vector Γ p is accepted as the new set of parameters.

ψ M -H (Γ c ,Γ p ) = min (1,exp[(log p(Γ p |ρ) -log p(Γ c |ρ)]) (10) log p(Γ|ρ) ∝ log|I N -ρW c (Γ)| - N 2 log[ω ′ F ω] (11) F = (ỹ -Xβ d ) ′ (ỹ -Xβ d ) β d = (X ′ X) -1 (X ′ ỹ) (12)
We note that F consists of only sample data, so this expression can be calculated prior to MCMC sampling, leading to a computationally efficient expression reflecting a quadratic form: log(ω ′ F ω), that can be easily evaluated for any vector of dependence parameters ω. Since ω involves both ρ and Γ, we can also sample the conditional distribution for ρ assuming the current vector of values Γ c is known using the same computationally efficient expression involving the quadratic form.

Of course, the conditional distribution also involves the log-determinant log|I N -ρW c (Γ)| which depends on both ρ and the vector Γ, and needs to be calculated for any given values of these parameters in order to evaluate the M-H acceptance probability in Eq. ( 10), which is the subject of the next section.

An approximation to the log-determinant

Pace and LeSage (2002) set forth a Taylor series approximation for the logdeterminant log|I N -ρW c (Γ)|, which appears in the conditional distribution needed to sample both ρ and the parameters in the vector Γ. They show that for a nonnegative weight matrix W c (Γ) with eigenvalues µ min ≥ -1, µ max ≤ 1, and 1/µ min < ρ < 1:

log|I N -ρW c (Γ)| = - ∞ ∑ j=1 ρ j trW c (Γ) j j (13) ≃ - q ∑ j=1 ρ j tr(W c (Γ) j ) j ( 14
)
where tr represents the trace operator. Golub and van Loan (1996, p.566) provide the expression in Eq. ( 13), while [START_REF] Pace | Semiparametric maximum likelihood estimates of spatial dependence[END_REF] note that due to the linearity of the trace operator we have the qth order approximate expression in Eq. ( 14). We note that the 1st-order trace, tr(W c (Γ)), is zero for any convex combination of weight matrices that have zero diagonal elements. Since the underlying weight matrices W ℓ have zero diagonal elements, a convex combination of these will also have this property. The second-order trace can be computed as in Eq. ( 15), and there are similar expressions for higher-order traces shown in Eqs. ( 16) and ( 17). Debarsy and LeSage (2018b) carry out Monte Carlo experiments to show that estimates based on q = 4th order trace approximation have desirable properties, including low bias and good coverage of the true parameter values.

tr(W c (Γ) 2 ) = L ∑ i=1 L ∑ j=1 γ i γ j tr(W i W j ) (15) tr(W c (Γ) 3 ) = L ∑ i=1 L ∑ i=1 L ∑ k=1 γ i γ j γ k tr(W i W j W k ) (16) tr(W c (Γ) 4 ) = L ∑ i=1 L ∑ i=1 L ∑ k=1 L ∑ ℓ=1 γ i γ j γ k γ ℓ tr(W i W j W k W ℓ ) (17)
The computationally desirable aspect of these expressions is that the parameters γ ℓ are separated from the trace calculations involving the N × N dimensional matrices. This allows the trace calculations to be carried out prior to MCMC sampling draws. During MCMC sampling, new values for the parameters Γ and ρ are used in conjunction with the pre-calculated traces to rapidly evaluate the log-determinant term.

Efficient calculation of partial derivatives for the model

In addition to producing estimates for the underlying model parameters ρ, β, Γ, σ 2 , we need to compute the reduced form of the model to compute the partial derivatives of y with the relevant exogenous determinant to interpret its effect. For the standard simultaneous autoregressive model, LeSage and [START_REF] Lesage | Introduction to Spatial Econometrics[END_REF] have derived these partial derivatives and further develop scalar summary measures to interpret them. For the simultaneous autoregressive model with a convex combination of connectivity matrices, the partial derivative of the outcome variable with respect to the r th determinant is shown in Eq. ( 18) while the scalar summary measures, which represent own-and cross-partial derivatives that they label Eq. direct and Eq. indirect effects, are shown in Eqs. ( 19) and ( 21).

∂y/∂x r = = S r (W c (Γ)) (18) S r (W c (Γ)) = (I N -ρW c (Γ)) -1 β r = I N β r + ρW c (Γ)β r + ρ 2 W c (Γ) 2 β r + . . . M (r) direct = N -1 tr(S r (W c (Γ))) (19) M (r) total = n -1 ι ′ N S r (W c (Γ))ι N (20) M (r) indirect = M (r) total -M (r) direct (21) W c (Γ) = L ∑ ℓ=1 γ ℓ W ℓ
While expressions in Eqs. ( 19), ( 20) and ( 21) produce point estimates for the scalar summary measures of effects (own-and cross-partial derivatives) used to interpret the impact of changes in explanatory variables on dependent variable outcomes, we also require measures of dispersion for the purpose of statistical tests regarding the significance of these effects. Use of an empirical distribution constructed by simulating the non-linear expressions in Eq. ( 18) using (say 1,000) draws from the posterior distribution of the underlying parameters ρ, β r , γ ℓ ,ℓ = 1, . . . ,L is suggested by [START_REF] Lesage | Introduction to Spatial Econometrics[END_REF].

A naive approach to such a simulation-based empirical distribution would require calculation of the N × N matrix inverse S r (W c (Γ)) a large number (say 1,000) times, for varying values of the parameters ρ, β r , γ ℓ , ℓ = 1, . . . ,L, which would be very compute intensive. For the SAR model with only one connectivity matrix, W , LeSage and [START_REF] Lesage | Introduction to Spatial Econometrics[END_REF] show that the required quantity for constructing the empirical distribution of the effects is tr(S r (W )), which can be estimated without a great deal of computational effort. Debarsy and LeSage (2018b) extend this approach using estimated traces in conjunction with the four first orders traces evaluated at each pass through the sampler to approximate the log-determinant term. This allows the simulation-based empirical distribution of the direct, indirect and total effects estimates to account for uncertainty in Γ, the parameters of the convex combination. It is important to take into account the stochastic uncertainty that arises as a result of these parameters and the additional variation/uncertainty they contribute to the effects estimates because ignoring it leads to under-coverage of the other parameters (too narrow standard errors).

Bayesian model averaging for multiple models

The convex combination of weights model specification raises questions regarding which weight matrices should be used in these models. For a modeling situation where L weight matrices are considered as candidates for the model specification, there are M = 2 L -L -1 possible ways to employ two or more of the L weight matrices in alternative model specifications. When L = 5, we have M = 26 possible models involving two or more weight matrices, and for L = 10, M = 1,013. This makes a method for comparing model estimates based on alternative combinations of weight matrices an important issue. Debarsy and LeSage (2018b) take a Bayesian approach to deriving the log-marginal likelihood for these models. Given the log-marginal likelihood for a model M i , that we denote LogM i , we can calculate its associated posterior probability, namely prob

(M i ) = exp(LogM i / ∑ Q i=1 LogM i ) (in the case of Q different models
). The highest prob(M i ) answers the questionwhich model is most consistent with the sample data (y,X,W 1 , W 2 , . . . , W L )? -unconditional on any parameter value/estimates. The answer is unconditional with respect to the parameter values because the parameters have been integrated out of the joint likelihood to produce the log-marginal likelihood and associated model probabilities.

Debarsy and LeSage (2018b) combine analytical and numerical strategy to produce estimates of the log-marginal likelihood. This involves analytically integrating out the parameters β, σ 2 from the joint likelihood to produce a joint posterior distribution for the dependence parameters ρ, Γ. They then rely on Metropolis-Hastings guided Monte Carlo integration during MCMC estimation of the models to produce estimates of log-marginal likelihoods for each model.

The Bayesian solution to the problem of uncertainty regarding which model specification is correct involves model averaging, where estimates from all candidate models are weighted by the posterior model probability. This means that model averaged estimates incorporate uncertainty regarding the model parameters as well as uncertainty regarding the model specification.

To illustrate this approach, we use a sample of 856 nursing homes located in Texas taken from [START_REF] Blankmeyer | Peer-group dependence in salary benchmarking: A statistical model[END_REF]. They argue that CEO compensation depends on that of Chief Executive Officers (CEO's) from peer nursing homes, where peers can be defined using alternative weight matrices. Specifically, we consider a spatial weight matrix based on 3-nearest spatial neighbors (W space ), a size similarity matrix based on three nursing homes with the most similar square foot area (W size ), and a payroll similarity matrix constructed from the three nursing homes with the most similar payroll magnitudes (W pay ). We note that peer nursing homes defined using the size and payroll measures could be located anywhere in the state of Texas, which comprises the sample of 856 nursing homes.

The explanatory variables are the share of medicaid and private pay patients, the occupancy rate, and dummy variables indicating for-profit nursing homes and nursing homes that are part of a chain operation. [START_REF] Blankmeyer | Peer-group dependence in salary benchmarking: A statistical model[END_REF] indicate these are criterion on which the executive board evaluates the CEOs managerial performance.

Given L = 3 different candidate weight matrices, we have M = 4 possible models involving two or more weight matrices. The log-marginal likelihoods, model probabilities and estimates for the parameters ρ, γ 1 , γ 2 , γ 3 are shown in Table 1 for the M = 4 possible models, along with Bayesian model averaged estimates (BMA).

The highest model probability of 0.659 points to Model 3 that includes W payroll (three most similar payroll peers) and W size (square foot area sizebased peers), placing no weight on the model that includes W space (the 3-nearest spatial neighboring nursing homes). The second highest logmarginal likelihood value is associated with Model 1 that has a model probability of 0.196 and includes W space and W payroll only, excluding the area-based size peers W size . The BMA estimates will be based on MCMC draws from Model 1, Model 3, and Model 4 which have non-zero posterior probabilities (0.196, 0.659, 0.145), respectively. These are reported in Table 2. From the table, we see evidence of significant dependence on peer institutions in explaining variation in CEO compensation, the dependent variable, as evidenced by the estimate for ρ = 0.5225, which is significant. The most important peer nursing homes for determining CEO compensation are other nursing homes with similar payroll magnitudes indicated by the estimate for γ 2 = 0.8644. This indicates that management of staff and nursing personnel is viewed as the most important basis for determining a peer institution by the executive board that evaluates and determines the CEO's compensation. Larger payrolls point to a larger CEO role in managing personnel at the nursing home.

The estimate of γ 3 = 0.1031 for the weight matrix associated with the matrix W size based on the physical size of the nursing home is the next most important basis for determining a peer institution, with a very small the role played by spatial proximity, since γ 1 = 0.0329, with a lower 0.01 credible interval of 0.0074, very close to zero.

Turning to the effects estimates, we see that a larger number of Medicaid and private pay patients in the nursing home (which are considered less desirable than Medicare patients) has a negative direct impact on CEO compensation. By direct impact, we mean that the own-institution characteristic impacts the own-institution CEO compensation, in a negative way for this variable. Higher occupancy rates have a positive impact on CEO compensation, and [START_REF] Blankmeyer | Peer-group dependence in salary benchmarking: A statistical model[END_REF] point to low occupancy rates as a problem for Texas nursing homes. The dummy variable selecting For profit institutions has a positive and significant impact on CEO compensation, indicating that these institutions have higher compensation levels than non-profit institutions, many of which are affiliated with religious institutions. Finally, the dummy variable selecting nursing homes that are part of a chain of homes administered by a single company, has a positive direct The indirect effects in this model reflect peer effects, such that compensation of CEO's at nursing home i is influenced by compensation levels for CEOs at peer nursing homes j ̸ = i. The median estimates for indirect effects are similar in magnitude to the direct effects, indicating that compensation of other nursing home CEOs plays an equally important role as the own-institution characteristics in determining compensation for the typical CEO in our sample. The total effects are of course the sum of direct plus peer/indirect effects.

Conclusion

The simultaneous autoregressive model specified using a convex combination of dependence matrices represents one solution to the important problem of selecting an appropriate dependence set of observations for use in modeling cross-sectional dependence. Using a specification based on a convex combination of dependence matrices allows formal Bayesian tests for the relative importance of alternative definitions of the dependence set of observations. The specification also allows a combination of several weight matrices that avoids scaling issues that arise in specifications based on multivariate similarity criterion such as those set forth in [START_REF] Blankmeyer | Peer-group dependence in salary benchmarking: A statistical model[END_REF]. Another advantage is avoidance of complex parameter space constraints inherent in higher order simultaneous autoregressive models. Debarsy and LeSage (2018b) propose a solution to the computational challenges inherent in simultaneous autoregressive model specifications based on convex combinations of weight matrices and further develop a Bayesian model averaging approach that allows for model uncertainty regarding alternative definitions of the dependence set. For example, one can draw inferences regarding which matrices should enter the specification and which should be excluded.

This approach should be useful in applied research where there is little theoretical guidance regarding an appropriate dependence set and associated weight matrix. An important caveat is that practitioners should attempt to consider definitions of dependence sets that convey differing types of potential cross-sectional dependence relationships. Use of alternative weight matrices that convey similar cross-sectional dependence relationships can encounter difficulties when simultaneous lags of the dependent variable vector (W 1 y, W 2 y, . . . , W L y), disturbances (W 1 u, W 2 u, . . . W L u), or explanatory variables (W 1 X, W 2 X, . . . W L X) exhibit high levels of correlation.

LeSage and Pace (2014) discuss the issue of highly correlated spatial lags that can arise in spatial autoregressive model specifications, and difficulties that arise in distinguishing between models based on similar definitions of dependence sets of observations. They label the mistaken belief that small changes in the definition of dependence sets can lead to large changes in model estimates and inferences as the the biggest myth in spatial econo-metrics. The same caveat applies to simultaneous autoregressive model specifications described here. Practitioners should strive for inclusion of alternative weight matrices in the convex combination of weights that reflect distinctly different types of dependence. Use of specifications where multiple weight matrices result in similar dependence sets of observations will lead to poor estimates and inferences regarding the relative importance of each set of weights.

Table 1 :

 1 Model probabilities for models involving two or more weight matrices

	Models	LogM Prob	ρ	γ 1 , W space γ 2 , W payroll γ 3 , W size
	Model 1 -455.822 0.196 0.506	0.102	0.898	-
	Model 2 -560.499 0.000 0.144	0.363	-	0.637
	Model 3 -454.611 0.659 0.518	-	0.871	0.129
	Model 4 -456.126 0.145 0.567	0.090	0.792	0.117
	BMA	-455.068 1.000 0.523	0.033	0.865	0.102

Table 2 :

 2 Model averaged estimates

	Variables	Lower 0.01 Lower 0.05 Median Upper 0.95 Upper 0.99
	Constant	3.8810	4.0613	4.6649	5.2588	5.4205
	Medicaid+Private	-0.2191	-0.1978 -0.1178	-0.0417	-0.0180
	Occupancy rate	0.7193	0.7494	0.8234	0.9009	0.9224
	For profit	0.0342	0.0452	0.0807	0.1185	0.1296
	Chain	0.0044	0.0130	0.0460	0.0813	0.0933
	ρ	0.4544	0.4696	0.5225	0.5797	0.5973
	γ 1 , W space	0.0074	0.0110	0.0329	0.0557	0.0629
	γ 2 , W payroll	0.7773	0.7965	0.8644	0.9366	0.9496
	γ 3 , W size	0.0215	0.0355	0.1031	0.1658	0.1847
	Direct	Lower 0.01 Lower 0.05 Median Upper 0.95 Upper 0.99
	Medicaid+Private	-0.2366	-0.2131 -0.1276	-0.0454	-0.0196
	Occupancy rate	0.7839	0.8117	0.8904	0.9711	1.0014
	For profit	0.0373	0.0489	0.0871	0.1283	0.1401
	Chain	0.0047	0.0138	0.0499	0.0878	0.1016
	Indirect	Lower 0.01 Lower 0.05 Median Upper 0.95 Upper 0.99
	Medicaid+Private	-0.2344	-0.2039 -0.1192	-0.0434	-0.0150
	Occupancy rate	0.6346	0.6804	0.8416	1.0602	1.1607
	For profit	0.0362	0.0462	0.0825	0.1297	0.1466
	Chain	0.0057	0.0135	0.0477	0.0880	0.1006
	Total	Lower 0.01 Lower 0.05 Median Upper 0.95 Upper 0.99
	Medicaid+Private	-0.4571	-0.4141 -0.2471	-0.0907	-0.0337
	Occupancy rate	1.4662	1.5341	1.7364	1.9897	2.1080
	For profit	0.0764	0.0960	0.1709	0.2569	0.2808
	Chain	0.0113	0.0272	0.0982	0.1735	0.1973
	impact on CEO compensation.

Cross References

LeGallo, J. , Chapter 1. Cross-Section Spatial Regression Models.

Mills, JA and Parent O. Chapter 3. Bayesian MCMC Estimation.

LeSage, JP and Pace, RK Chapter 8. Interpreting Spatial Econometric Models.


Ali Hodroj

Marc Ibrahim

Yassine Hadjadj-Aoul

Bruno Sericola

Multi-armed bandit algorithms over DASH for multihomed client

Keywords: DASH, Multi-homed, video streaming

There is an expanding request from the mobile clients of video traffic which has come to over the most recent two years the greater part of all the mobile data traffic. The challenge of providing the clients excellent video streaming is expected principally to the constrained transfer speed which makes it important to misuse the rising decent variety of mobile video spilling and get to systems. For enhancing the video quality received by "multi-homed client" (i.e., clients with multiple network interfaces), a network selection algorithm based on the Multi-Armed Bandit heuristic is proposed on top of the most broadly utilized standard for video streaming Dynamic Adaptive Streaming over HTTP (DASH). For choosing the best network at each progression powerfully, DASH gives mobile video quality dependent on the apparent exhibition from the pre-owned system association through the Adaptive Bitrate Rules (ABR) without investigating the system conditions through the other network(s) which could give better quality. Subsequently, a few alterations for DASH ABR is required to enhance the video quality. Two of the MAB algorithms (UCB and Epsilon Greedy) were embraced for improving MPEG-Dash. The investigations are performed through a proving ground execution to show that UCB surpasses Epsilon Greedy, in stable system conditions, regarding goodput received by the Dash customer. Additionally, UCB can discover the harmony between investigating new choices, and abusing the triumphant variation.

I. INTRODUCTION

Mobile traffic will be multiplied in the coming scarcely any years contrasted with 2016 while 82 percent of the information traffic will be for video streaming as indicated by Cisco forecasting [START_REF]Cisco. Cisco visual networking index: Forecast and trends[END_REF]. The boundless utilization of smartphones and the developing mainstream video stages have contributed to this remarkable development. Such expanding utilization of video applications is trying from various outlooks from the framework suppliers to the substance suppliers and not least the spilling stages (e.g., Netflix). To manage the variances of the system execution, the Dynamic Adaptive Streaming over HTTP (DASH) convention [START_REF] Mueller | Dynamic adaptive streaming over http/2.0[END_REF] was proposed to adjust the imagined substance to the system execution, so as to upgrade the nature of experience (QoE) received by the end clients. Fundamentally, the video content is made out of numerous segments having diverse quality levels which are portrayed in an a media presentation description (MPD) record. The DASH customer demands, through HTTP protocol, the video segments with different quality levels that fit with the present status of network conditions and adaptive bitrate streaming rules [START_REF] Lee | Enhanced quality adaptation scheme for improving qoe of mpeg dash[END_REF].The outcomes are that dispersed network clients will watch similar video content from a similar server in various quality relying upon end-to-end throughput on the comparing system ways with the server and the present cushion level [START_REF] Go | Hybrid tcp/udp-based enhanced http adaptive streaming system with multi-homed mobile terminal[END_REF]. Then again, the decent variety of access network interfaces (e.g., Wi-Fi and cell arrange) on smartphones these days, gives such multihomed gadgets the ability to change starting with one network connection then onto the next dependent on the current network conditions level and to improve the video streaming quality. Nevertheless, there are two issues that can be distinguished in the plan. First issue is that the way asymmetry of heterogeneous network can have various attributes and for this situation exchanging between the various networks may debase the general video quality. The subsequent issue is that the network exchanging could negatively affect the mobile technique of DASH which thinks about the present cradle level and throughput to choose the nature of the following segment. Therefore, we suggest in this article an adaptive network selection algorithms based on the Multi-Armed Bandit algorithms which permit a multihomed Dash use to switch between different networks (WiFi, 3G, and 4G). This solution aims to provide a new adaptive bitrate rule that aims to reduce the effect of network selection and achieve better video streaming quality perceived by the mobile users. This 19 is achieved via the following functions of the proposed solution:

• Applying it at the application level which does not require any modification at the device kernel or server sides.

• Implementing the MAB algorithms in the DASH player of multihomed users which affords an adaptive network selection mechanism that takes benefit of the diversity of the network. • Leveraging the availability of a couple of resources as video streaming mirror servers to address the excellent degradation that could be caused by the visitors congestion at server side. [START_REF] Ma | Understanding the power of smartrouter-based peer cdn for video streaming[END_REF] • Adjusting the DASH ABR rules to address the channel diversity. This journal is prepared as follows. Section II provides related studies. Then, our proposed fashions for MAB implementation over DASH is elaborated in Section III. Section IV discusses the measurement scenarios and evaluates our models. Finally, realization and future works are presented in Section V.

II. RELATED WORK

Given that our work falls in the domain of multihoming , multipathing and multi-source over DASH protocol in Content Delivery Networks (CDN), there is no contribution covering every one of these angles together at the application layer. Nonetheless, many arrangements were proposed at the different layers for the presentation improvement of video streaming applications received by the multihomed clients through the simultaneous multipath transmission of media content. We present in this area such commitments in the wake of arranging them into two classes which are application and transport layer arrangements.

A. Application Layer

Sub-frame level (SFL) scheduling approach is introduced in [START_REF] Huang | Notice of violation of ieee publication principles: Group-aware delay-constrained video transmission over multihomed device-to-device networks[END_REF] which splits large-size video frames to optimize the delay performance of HD video streaming over heterogeneous wireless networks. It uses a water-filling algorithm for HD video frame splitting and streaming rate allocation over multiple wireless access networks. Also it uses a FEC coding to mitigate wireless channel losses taking the advantage of the flow rate allocation algorithm. This work was the first work to perform the frame-level data protection and transmission of delay-stringent HD video streaming in heterogeneous wireless networks with multi-homed terminals.

In [START_REF] Chen | Msplayer: Multi-source and multi-path video streaming[END_REF], MSPlayer is proposed as a client-based video streaming ,that doesn't require any alteration at TCP stack,and exploits various video sources just as network paths through various interfaces. MSPlayer intends to decrease start-up latency and provide high-quality video streaming. MSPlayer streams videos from multiple YouTube video servers through two interfaces (e.g., WiFi and LTE) at the same. They thought about the exhibition of standard YouTube players and MSPlayer dependent on pre-buffering and re-buffering stages. MSPlayer outperforms both single-path TCP over WiFi and LTE for different specified amounts of pre-buffered video data. Additionally, MSPlayer decreases the start-up delay by 12%, 21%, 28%, individually, for 20/40/60 second buffers.Nonetheless, they don't report the outcomes on how MSPlayer gives heartiness video conveyance in mobile situations.

The author of [START_REF] Wu | Loss tolerant bandwidth aggregation for multihomed video streaming over heterogeneous wireless networks[END_REF] propose an application-layer transmission framework based on a new Loss Tolerant Bandwidth Aggregation (LTBA) scheme. LTBA consist of is spreading the FEC packets for overcoming the burst channel losses to minimize end-to-end video distortion and delay. The results show that LTBA enhances the average video received quality and guarantee that the video frames are delivered on time even in high loss wireless channel.

GreenBag [START_REF] Bui | Greenbag: Energy-efficient bandwidth aggregation for real-time streaming in heterogeneous mobile wireless networks[END_REF] is an energy-efficient bandwidth aggregation middleware that supports real-time data-streaming services over asymmetric wireless links, requiring no modifications to the existing Internet infrastructure and servers. GreenBag is devised as a middleware operating on the mobile device, requiring neither a proxy server nor any modification to the existing Internet infrastructure and servers. It provides energyefficient bandwidth aggregation service for real-time video streaming over LTE and Wi-Fi interfaces, adapting dynamically to network fluctuations. It is a multi-link data streaming middleware for supporting reliable and energy efficient realtime data-streaming services, adapting dynamically to varying network conditions. GreenBag aims to provide an optimal mode switch between single-link and dual-link modes, in terms of minimizing energy consumption while meeting the required Quality of Experience (QoE).

B. Transport Layer

A well-known transport layer protocol for multipath bandwidth aggregation is Multi-Path TCP (MPTCP) which works on balancing the load over the different network paths between the server and the user's network connections by splitting the data stream across the multiple network interfaces [START_REF] Cao | Towards efficient parallel multipathing: A receiver-centric cross-layer solution to aid multipath tcp[END_REF]. It consists of multiple independent TCP connections (subflows). The concurrent use of these links for a TCP/IP session would improve resource usage within the network and, thus improve the user experience through higher throughput and improved resilience to network failure. However, MTCP has a problem that quality leads to degrading due to different characteristics (bandwidth, latency, packet loss rate, etc.) of available asymmetric paths which cause packet out-of-order problem. To solve this problem many intelligent scheduling algorithms at the sender have been proposed to reduce the number of disordered packets. [9] [5].

Kim [START_REF] Wu | Streaming highquality mobile video with multipath tcp in heterogeneous wireless networks[END_REF] proposed a scheduling algorithm at the sender to reduce the receive buffer. The main idea is to estimate packet arrival time, by calculating the time that every path to transmit data packet before the transmission of a data packet and schedule packets accordingly. The sender maintains a per path timetable including calculated values of receiving time at the receiver from now for each packet. When the sender has the opportunity to send a new packet, it chooses the path owned smallest transmission time to transmit this packet.

The authors in [START_REF] Agarwal | Cognitive multihoming: Maximizing network utility over cr-assisted cellular network[END_REF], [START_REF] Preethi Natarajan | Concurrent multipath transfer during path failure[END_REF], [START_REF] Wu | Quality-aware energy optimization in wireless video communication with multipath tcp[END_REF]- [START_REF] Wu | Content-aware concurrent multipath transfer for high-definition video streaming over heterogeneous wireless networks[END_REF] propose different transport level solutions for concurrent multipath video streaming for multihoming users. In [START_REF] Wu | Streaming highquality mobile video with multipath tcp in heterogeneous wireless networks[END_REF], Jiyan Wu propose ADMIT a quality-driven Multi-path TCP that use Forward Error Correction (FEC) coding and rate allocation to reach best quality for real-time streaming. ADMIT uses a rate allocation algorithm to select the appropriate access network and then uses a FEC coding adaptation scheme to maximize the video quality and balance between delay and loss performance. However, it suffers from delay performance degradations that may occur due to the retransmission mechanism in MPTCP when packet loss occurs. A Distortion aware concurrent multipath transfer scheme, proposed in [START_REF] Agarwal | Cognitive multihoming: Maximizing network utility over cr-assisted cellular network[END_REF], that aims to reduce the end-to-end video distortion rate by relying on path status estimation, flow rate allocation, and re-transmission control. The available bandwidth on each path is estimated by sending ACK packets one each path. A simple queue is used to characterize packet delay over each component link. In [START_REF] Wu | Content-aware concurrent multipath transfer for high-definition video streaming over heterogeneous wireless networks[END_REF], a content aware scheduling technique is proposed to improve the video quality by assessing the video content parameters (e.g., video distortion level, data loss rate, and end to end delay of the frame) in an appropriate route for scheduling the video frame using a decision procedure at server side. A TCP friendly congestion control algorithm is used and packets are distributed to different paths on frame level. In [START_REF] Wu | Quality-aware energy optimization in wireless video communication with multipath tcp[END_REF], the authors propose energy-distortion mechanism to provide quality-guaranteed video streaming with minimization energy consumption of mobile devices over heterogeneous networks. The minimization of energy consumption was reached though video flow rate allocation algorithm using an analytical energy distortion framework under the target video quality constraint. The result shows that for the same video quality EDAM can with lower energy consumption compared to the energy consumed using the traditional MPTCP. A path failure (PF) detector is proposed in [START_REF] Preethi Natarajan | Concurrent multipath transfer during path failure[END_REF] were paths that experience multiple timeout events are classified as "Potential failed" (PF) and not used for data transmission, while the data will be delivered using the active paths. Concurrent multipath transfer (CMT)-PF reduces the detection latency of link failures and improves CMT's throughput. The current solutions enable the flows or paths to be setup between a client and a single server without exploiting service replication, which is important to prevent server overload issue because of the high demand for high quality video. Besides, transport layer solutions require kernel modifications at both server and client sides, in addition to the limitation that many network operators do not allow MPTCP traffic to pass through their middleboxes. Concerning the current solutions at the application level, they require server side modification without taking the upsides of dynamic adaptive video streaming over HTTP when using the multipath communication. Subsequently, there is a need to build up an efficient video streaming solution from multisources to multihomed mobile users implemented over the DASH standard.

III. PROPOSED SOLUTIONS

We develop a selection mechanism between different networks to exploits multihoming feature to maximizing the QoE of the DASH standard clients. Basically, the use of DASH permits the video content to be stored as segments of different qualities known as representations. The DASH player start streaming the content by loading the Media Presentation Description (MPD) which is a manifest file that describe the video content (video and segment length, video qualities and etc.). Then the player start downloads the video segments in sequence over HTTP where an Adaptive Bit rate (ABR) is executed to permit the client to choose different qualities of downloaded segment according to its current network conditions [START_REF] Vlaović | Evaluation of adaptive bitrate selection algorithms for mpeg dash[END_REF]. Thus, the client can watch the video at the highest available quality, and continuously without any freezes. The client has the choice to decrease the bitrate when the network throughput degrades and increase the bitrate when the network throughout improves. In Dash design, there is a buffer at client side permit to store a number of segments then played out when the user requires. Rebuffering occurs when the buffer is empty and not more segments to play out, thus the video playback must freeze until receiving the nest segments. ABR algorithms are classified into three types: throughput, buffer, and hybrid rules. Throughput Rule like PANDA [START_REF] Li | Probe and adapt: Rate adaptation for http video streaming at scale[END_REF] and Festive [START_REF] Jiang | Improving fairness, efficiency, and stability in http-based adaptive video streaming with festive[END_REF] focus on the available network throughput between the player and the server where the requested quality of next downloaded segment is based on the estimation of current network throughput. Buffer Rule like Bola [START_REF] Spiteri | Bola: Near-optimal bitrate adaptation for online videos[END_REF] and BBA [START_REF] Huang | A buffer-based approach to rate adaptation[END_REF], use the buffer level as indictor of network throughout. It chooses higher quality when the buffer level is higher and lower when the buffer level is lower. Hybrid Rule like MPC [START_REF] Xiaoqi Yin | A control-theoretic approach for dynamic adaptive video streaming over http[END_REF] and ELASTIC [START_REF] De Cicco | Elastic: A client-side controller for dynamic adaptive streaming over http (dash)[END_REF] use both throughput and buffer level to take the advantages of both. A key goal of our model is to provide higher QoE in terms of higher bitrate and fewer rebuffers, to a multimomed mobile dash client taking the advantage of CDN where the video can be replicated at multiple servers. Multi-Armed bandit algorithms are implemented on client DASH protocol, permitting the user to watch the video in highest quality by switching between different available networks.

A. Model

• t C : Segment duration.

• D: Video duration.

• C: Number of segments , which is equal to D t C . • k: Time is divided into steps with a duration of K. The duration of the time step is equal to c×t C , where c is the fixed number of requested segments to be downloaded in a time step.

• I = {I i } i=1..N I : Available networks N I • n i (k): Number of times that network i is selected until time step k. • d(k):
The selected network at time step k (i.e. decision).

• r(k): The reward received at time step k, which is the average bitrate measured at network d(k).

• Q = {q j } j=1..N Q is the set of available N Q video bitrates.

B. MAB Algorithm-Epsilon Greedy

Algorithm 1 MAB Algorithm-Epsilon Greedy

1: q init = DEFAULT 2: for i = 1 to N I do bootstrap 3:
r i , q init = getNSegments(c, q init , i) r(k), q init = getNSegments(c, q init , d(k))

12:

r d(k) = αr(k) + (1 -α)r d(k)
13:

k = k + 1
MAB algorithms [START_REF] Zhang | Cooperative relay selection for load balancing with mobility in hierarchical wsns: A multi-armed bandit approach[END_REF] aim mainly to select the network which maximize the QoE. we start by the -greedy algorithm which is one on the main algorithms behind science of decision, and symbolizes the balance of exploration versus exploitation. The algorithm keeps track of the cumulative average reward r i of each network i in I, that is given by:

r i = αr(k) + (1 -α)r i
At each time step k, exploitation or explorations may occur.Exploitation consist of choosing the network that has the best average return based on what we know of the circumstances with a probability 1-. exploration consist of choosing a random network with probability (see lines 6 to 9) using the function the latter case, a discover phase is allowed by selecting a random network, using the getRandNetwork(I).

The parameter represent the percentage of time/trials dedicated for exploration. It is typical be around 10%, but in our experiments we study the impact of this parameter.For the selected network i, the function getNSegments(c, q init , i) is used to download c segments, q init being the quality of the first segment. The qualities for subsequent segments will be adaptively determined by the ABR algorithm. The funtion getNSegments(c, q init , i) returns the reward, which is the average bit rate measured at network i, and the average measured quality. In this section, the quality of requested segments is determined according to the ABR rules used by default in DASH as described before (Throughput and Buffer rules).

C. MAB algorithm -Epsilon Greedy with ABR modification

According to abr rules, playing the video at the highest available quality and continuously without any freeze are the main factors that improve the QoE. In the case of multihoming, the presence of multiple networks may cause some false choices of quality due to different network characteristics. Request a segment at a high bitrate, based on the assumption on first network conditions, through a different network that may have lower throughput than the first network, will cause the client to rebuffed because the player is unable to fetch this required segment at the quality requested. In this real scenario, the characteristics of different available networks should be an additional indicator to ensure that ABR is well working and improve the QoE. A simple modification is proposed in Algorithm 2, to solve this problem, consist of setting the quality of the first requested segment on the switched network equal to the last quality measured on the same network.

r(k), q init (d(k)) = getNSegments(c, q init (d(k)), d(k)) 11: r d(k) = αr(k) + (1 -α)r d(k)
12:

k = k + 1
In the next paragraph, we study the performance of MAB-Epsilon Greedy without abr modification in different networks scenarios. Then in case the described problem appear, we compare the results to MAB-Epsilon Greedy with abr modification as mention in Algorithm 2.

D. MAB algorithm-Upper Confidence Bounds (UCB) with ABR modifictaion

Upper Confidence Bounds (UCB) is a deterministic algorithm that focuses on exploration and exploitation based on a confidence boundary that the algorithm assigns to each network on each time of exploration. UCB is based on choosing the lever that promises the highest reward under optimistic assumptions (i.e. the variance of your estimate of the expected reward is high, therefore you pull levers that you don't know that well) [START_REF] Boldrini | Introducing strategic measure actions in multi-armed bandits[END_REF]. UCB is an allocation strategy achieving logarithmic regret uniformly over n and without any knowing in advance about the reward distributions. It can be briefly said that UCB's main idea is: the principle of optimism in the face of uncertainty in spite of our lack of knowledge in what actions are best. Our optimistic guess will quickly decrease and we'll be urged to switch to a different action (network) in case the guess was wrong. But if we pick well, we'll be able to experience little regret. UCB algorithm is used in multi armed bandit problem where given a set of arms, our goal is to maximize the reward got after using every arm. UCB explores all arms then starts exploiting. The client is multihomed, has multiple interfaces refer to multiple networks or connections. The arms are the networks. These networks can have different throughput's which ensures different rewards. The reward is inferred from the quality of the requested segment. Since the video player does not switch the quality immediately, UCB start by exploring all the networks then try to maximize the quality of the played video. For each available network, an upper-confidence bound is defined that shows the highest guess at the payoff for that arm. With the UCB algorithm, the arm (network) with the highest UCB calculated index is chosen for each transmission. The arm with the highest probability of being "the best arm" is being played. The reward for each network is defined Equation ( 1) in shows the highest guess at the payoff for that network as well as the number of times the network has been chosen.

r i = r i + (2 log 2 ( C c n i (k) )) (1)
As described in Algorithm 3, the bootstrap phase consists of download c chunks at each available network i and calculate the reward of each network. After all the networks have been used, it is the UCB values (reward) turn to determine the network to be selected. After every time step(k), the function getNSegments(c, q init , i) is used to download c segments using the selected network i, q init being the quality of the first segment. Then, the UCB value (reward) is updated for every network. we set the quality of the first requested segment on switched network equal to the last quality detected on the same network. Finally, select the network with the highest UCB and switch to. r i , q init (i) = getNSegments(c, DEFAULT, i)

n d(k) = +1 8: r(k), q init (d(k)) = getNSegments(c, q init (d(k)), d(k)) 9: r d(k) = r d(k) + (2 log 2 ( C c ni(k) ))
10:

k = k + 1
IV. ANALYTICAL PART Moving from theoretical algorithms to applied successful algorithms requires a set of experiments that take into consideration all possibilities. In our proposal, as mention in the Introduction part, we should discover the performance of the proposed algorithms in channel diversity different network conditions between different channels) and intrachannel diversity (network variations inside the same channel). The performance evaluation of our proposed model is achieved through the implementation of an experimental test-bed in order to avoid the dynamicity of the Internet and be able to generate the same network configurations suitable for validating the experimentation scenarios under question. Wondershaper tool is used to limit the throughput of networks, it allows us to set the maximum download rate and/or maximum upload rate of an interface (refer to a network) from the command line.

A. Reference Measurements

we start our experiments by a set of reference measurements.

1) Dash performance in different network speeds: The objective of these measurements is to evaluate how the DASH protocol deal with network speed variations. In order to evaluate our proposed algorithms, we need to compare them to a reference basic scenario presented in Figure 1. In this scenario, a normal DASH client watches a video content of duration 10 minutes stored in a server with 10 different bitrates. In order to evaluate the key of using DASH protocol and standard ABR rule, the same scenario is repeated with different network throughputs controlled by the Wondershaper tool which permit us to set different download/ upload rate on each scenario. As mention before, in order to evaluate the QoE of the player, two main factors are required: the average video rate and re-buffering time. In each scenario, we measure the average video rate which is the video average quality (bitrate) received. Figure 2 presents the results showing that the highest quality could be reached in a fast network throughout and permit us to map the quality to the network throughout. The second indicator of the QoE is the rebuffering issue, resulting in counting the time of empty buffer occurrences and the duration of each occurrence. Figure 3 shows how many time we notice empty buffer and the duration of each appearance in different network throughout's. 2) MAB-EPS greedy varitations: In these experiments, our objective is to discover the best value of eps when using the Epsilon Greedy algorithm. In this second experiment set, presented in Figure 4, the client is Mutlithomed DASH player occupied with two network connections and using the same video content but replicated in two mirror server which permit us to cover the concept of CDN network. We implement algorithm 1 at the DASH client side. The Epsilon Greedy algorithm consists of a tradeoff between exploitation and exploration based on the "epsilon" parameter eps. This parameter represents the percentage of time/trials dedicated for exploration, and it is also typical to do random exploration. 81 scenarios were built and run (Figure 5) in order to determine the best value of eps. each scenario consists of watching the video content with different throughput's for network connection two (varies from S 0 to S i ) while fixing the throughput of the first connection to 5120 kbps. Based on the result of the first scenario, we choose 9 different sets of network throughput's based on the value of throughput that affects the average video bitrate. Then the same scenarios were repeated using 9 different valuess of eps. In Figure 6, the average video bitrate received is compared between different values of eps. The highest quality is reached when the value of eps is 0.1 which means we have 10% exploration while the quality decreased from 5Mbps to 3Mbps when we have 90% exploration (eps equal to 0.9).

B. MAB-Epsilon Greedy performance

In order to evaluate the performance of Epsilon Greedy as network selection protocol for a multi-homed Dash client, a new set of measurements were run using 81 different combination values of network throughput's as shown in Figure 8. For each measurement, we calculate the bitrate enhancement E (Equation 1) by measuring the average bit rate Q i when using network connection i on each scenario.

E = (E 1 + E 2 ) 2 = Q1-QMAB Q1 + Q2-QMAB Q2 2 (2)
Figure 9 illustrates the Cumulative Distributed Function (CDF) of the enhancement metric estimated from the performed scenarios when compared to standard Dash player using each one of these two networks alone and to MAB proposed without abr modification. The results show that using multihoming improve the QoE of user expect in 4 scenarios. Different enhancement ranges were observed, up to 80% enhancement in 70% from the measurements and 20% in the remaining measurements. C. DASH-Abr performance when using MAB-Epsilon Greedy

The same measurement set 2 scenarios are repeated after modifying the ABRrule in DASH in order to solve the freeze effected occurred due to different network throughput's between different networks. The results (Figure 11) show the client plays the video continuously without any freezes and the QoE is enhanced up to 85% in 80% of the scenarios. Thus, providing good QoE in the case of multihoming requires some modification in standard ABR rules.

D. MAB-Epsilon Greedy Vs MAB-UCB performance in intranetworks variation

Based on results in the previous paragraph, we use the same abr modification in Dash in both MAB algorithms (considered before in Algorithm 3). In this section, the same experiment environment presented in Figure 4 is used,we compare the performance of Epsilon Greedy algorithm to UCB algorithm, while the throughput of network 1 is fixed to 5120 kbps and running different scenarios with different speeds of the second connection. The results of Figure 12 show better quality is attained when applying the UCB algorithm especially when there is a big gap between both server's bitrates. For example, in the case where Server 1 bandwidth is 5 Mbps and Server 2 bandwidth is 512 Kbps, the average quality obtained is double with UCB used. We notice that when using the UCB method, the average quality doesn't vary much even when Server 2 bandwidth differs from 512 Kbps to 5 Mbps. This result is anticipated since the UCB chooses the most network speed available that has the highest reward. Furthermore, the video experience is affected by the empty buffer occurrence and the duration of the buffer's duration. Figure 13 show that UCB has lesser number of occurrence than Epsilon Greedy which can be interpreted with a small delay and rebuffering. Adding to that, the buffer's duration in UCB is minor compared to Epsilon Greedy which incurs fewer frequent freezing and annoying interruptions. The QoE of the user is also determined based on the count of freezes occurrence with different periods. Figure 13 shows that a maximum re-buffering duration of 17 seconds in UCB appears to be more desirable than the 66 seconds in Epsilon Greedy. In these experiments, our objective is to ensure the reliability and effectiveness of the performance evaluated between Epsilon Greedy and UCB algorithm during different internetwork variations, which permit us to ensure and validate the results. Each scenario consists of watching the video file with variable network speeds for network connection two (varies on different times as showing in Table 1), while fixing the speed of the first connection to 5 Mbps. Each scenario is repeated using MAB-Eplison Greedy and UCB algorithm. As illustrated in Figure 14, the average qualities perceived in approximately all scenarios were closed to each other in all cases expect in scenario 2. In scenario 2, the average bitrate in UCB was double the average bitrate in Epsilon Greedy. This scenario is when the server 2 speed has increased twice the server 1 speed. In our case, higher bitrate could be achieved when using UCB. Furthermore, as we notice, the average quality obtained in Experiment 1 is much greater than the one in Experiment 2 which assures that network variations affect negatively on the overall output and performance for the Dash client and degrade his QoE. As visualized in Figure 15 , the number of occurrences in Epsilon Greedy is superior to the one in UCB. We can also tackle the issue that especially in scenarios 3 and 4, where server 2 speed decreases. The duration gets higher in both cases, but UCB is much higher than Epsilon Greedy especially for the set of scenarios [START_REF] Cao | Receiverassisted partial-reliable multimedia multipathing over multi-homed wireless networks[END_REF][START_REF] Chen | Msplayer: Multi-source and multi-path video streaming[END_REF][START_REF]Cisco. Cisco visual networking index: Forecast and trends[END_REF][START_REF] De Cicco | Elastic: A client-side controller for dynamic adaptive streaming over http (dash)[END_REF] when more fluctuations in network Speed appear. Therefore, Epsilon Greedy's performance is much annoying to the user when the network flutters. For different speeds between different connections but stable network conditions on each connection,the results show that UCB exceeds Epsilon Greedy in terms of goodput perceived by the Dash client but empty buffer occurrences still exist in both situations.The -greedy algorithm takes the best action most of the time, but does random exploration occasionally while UCB ensures we select optimal actions as t → ∞ marking smart exploration. Thus, the UCB algorithm is able to find the balance between exploring unfamiliar options, and exploiting the winning variant. But for unstable network conditions on each connection, the results show that -greedy is more stable than UCB. Our future work will focus how use the availability of multiple connection to introduce the parallel streaming over DASH, which will be a new case of Multihoming to enhance the video quality perceived.

4: k = N I + 1 5 :

 5 while k < C c do while video not finished 6:I * = argmax i=1..N I r i

Algorithm 2 do while video not finished 5 : 6 :

 256 MAB Algorithm-Epsilon Greedy with Abr modification 1: for i = 1 to N I do bootstrap 2:r i , q init (i) = getNSegments(c, DEFAULT, i)3: k = N I + 1 4: while k < C c I * = argmax i=1..N I r i if rand() < 1 -then

Algorithm 3

 3 MAB Algorithm-UCB Algorithm with Abr modification 1: for i = 1 to N I do bootstrap 2:

3: k = N I + 1 4 :

 4 while k < C c do while video not finished 5: I * = argmax i=1..N I r i 6: d(k) = I * select the best network 7:

Fig. 1 :

 1 Fig. 1: Reference scenario parameters

Fig. 2 :

 2 Fig. 2: Video quality variation in function of network speed variation

Fig. 3 :

 3 Fig. 3: Duration of empty buffer occurrences as a function of network throughput variation

Fig. 4 :

 4 Fig. 4: Experimental environment

Fig. 5 :

 5 Fig. 5: Measurement Set 1

Fig. 6 :Fig. 7 :

 67 Fig. 6: Average Video quality vs Eps value

Fig. 8 :

 8 Fig. 8: Measurement Set 2 of Algorithm 1 implementation

Fig. 9 :

 9 Fig. 9: Average Bitrate Enhancement of Algorithm 1

Fig. 10 :

 10 Fig. 10: Empty Buffer Occurrence and Duration of Algorithm 1

Fig. 11 :

 11 Fig. 11: Average Bitrate Enhancement of Algorithm 2

Fig. 12 :

 12 Fig. 12: Average quality perceived using UCB and Epsilon Greedy

Fig. 13 :

 13 Fig. 13: Empty Buffer Occurrence and Duration of UCB and Epsilon Greedy

Fig. 14 :

 14 Fig. 14: Average quality perceived using UCB and Epsilon Greedy for testing scenarios

Fig. 15 :

 15 Fig. 15: Empty Buffer Occurrence and Duration of UCB and Epsilon Greedy for testing scenarios

TABLE I :

 I Network Variation Scenarios

	Scenario	Speed 2	Time(min)
	1	5 Mbps 10M bps	5
	2	5 Mbps 18M bps	5
	3	5 Mbps 1M bps	5
	4	5 Mbps 512Kbps	5
	5	5 Mbps 1M bps 4M bps 8M bps	2
	6	5 Mbps 8M bps 2M bps 4M bps	2
	7	5 Mbps 8M bps 10M bps 4M bps	2
	8	5 Mbps 2M bps 1M bps 4Kbps	2