Eric Gautier
email: gautier@ensae.fr

Bayesian Estimation of Inequalities with Non-Rectangular Censored Survey Data

Keywords: Inequality, Wealth distribution, Survey methodology, Bayesian statistics, Monte-Carlo Markov chains AMS 2000 Subject Classification: 62F15, 62D05, 65C40, 65C05

Synthetic indices are used in Economics to measure various aspects of monetary inequalities. These scalar indices take as input the distribution over a finite population, for example the population of a specific country. In this article we consider the case of the French 2004 Wealth survey. We have at hand a partial measurement on the distribution of interest consisting of bracketed and sometimes missing data, over a subsample of the population of interest. We present in this article the statistical methodology used to obtain point and interval estimates taking into account the various uncertainties. The inequality indices being nonlinear in the input distribution, we rely on a simulation based approach where the model for the wealth per household is multivariate. Using the survey data as well as matched auxiliary tax declarations data, we have at hand a quite intricate non-rectangle multidimensional censoring. For practical issues we use a Bayesian approach. Inference using Monte-Carlo approximations relies on a Monte-Carlo Markov chain algorithm namely the Gibbs sampler. The quantities interesting to the decision maker are taken to be the various inequality indices for the French population. Their distribution conditional on the data of the subsample are assumed to be normal centered on the design-based estimates with variance computed through linearization and taking into account the sample design and total nonresponse. Exogeneous selection of the subsample, in particular the nonresponse mechanism, is assumed and we condition on the adequate covariates.

Introduction

Approximately every six years the French statistical office INSEE collects a cross sectional wealth survey on households. The last dataset was collected in 2004. Several aspects can be studied focusing for example on holdings or particular 1 types of assets like the professional wealth or intergenerational transfers. One natural question concerns the nature of the distribution of wealth and its allocation in the various possible holdings. However, it is known to be a difficult one as well [START_REF] Juster | Improving the Quality of Economic Data: Lessons From the HRS and AHEAD[END_REF] due to the difficulty to have good measurements and to possible selection biases. Questions on income and wealth are particularly sensitive and the nonresponse probability is very likely to be related to the value itself, resulting in possible endogeneous selection and biases. Also, it is for example particularly difficult to give a precise amount for the market value of one's real estate piece of property unless people assessed it recently, say in order to sell it. Thus, amounts are usually collected in bracketed format and imputation methods (see experiment based on the knowledge of the true income distribution in [START_REF] Lollivier | D'une variable discrète à une variable continue : la technique des rsidus simuls[END_REF] are used in practice at the French institute.

For some variables the brackets are defined by each household, they give upper and lower bounds for the amount based on their evaluation. For other variables the households choose among a predefined system of brackets. The method allows to replace missing data, impairing inference due to selection bias and loss of efficiency implied by a reduced sample, by censored data. In the current article we focus on the evaluation of inequality indices on the total wealth for the whole French population. A specific question at the end of the survey: "Suppose you had to sell everything, how much do you assess the value of your total wealth including durable goods, artwork, private collections and jewelry" allows to measure the total wealth. The values of the last items were not collected. It is troublesome to ask such information because the pollster comes to the household's home and they could be suspected of theft in case a robbery occurs after the visit. The system of brackets for the question collecting the total wealth has an unbounded last bracket. The threshold for the higher bracket is 450.000 which is pretty low. Also in order to improve, in principle, the precision of design based estimates, certain categories have been over-sampled: self employed, executives, retired people and people living in rich neighborhoods. These variables, available from the census, are indeed correlated with the wealth and sampling more a priori wealthy people improve the precision of design based estimates of inequality indices sensitive to the top of the distribution. But due to the censoring, a billionaire is equivalent to an household whose total wealth is 451.000 . Over-sampling has increased the number of household for which we measure wealth very imprecisely. Thus, though we are interested in the particular concept of total wealth collected in this final question, we aim to gather more information in order to better estimate the inequality indices. Due to the high amount of censoring for the wealthiest, the less wealthy contribute the most to the likelihood. Moreover, some of the indices are quite sensitive to misspecification of the model. Though it is pretty usual for wages and even for wealth to specify linear models for the logarithm with normals residuals, some of the assumptions like the distributional assumption for the residuals are not testable in the absence of pointwise measurements. Other distributions like Pareto are also quite popular in the literature on wealth inequalities [START_REF] Lollivier | D'une variable discrète à une variable continue : la technique des rsidus simuls[END_REF]. It is also possible that the influence of certain covariates is not additive, for example the contribution of the income might differ for low incomes from that of high incomes, and, due to the censoring of high wealth, it is possible that we do not capture well this mixture (see for example [START_REF] Duclos | Polarization: Concepts, Measurement, Estimation[END_REF] and Esteban and Ray 1994 for theoretical foundations of polarization). Another possible source of misspecification is heteroscedasticity. If for example the variance of the residuals increase with income, censoring and specification of an homoscedastic model will imply lower inequality indices. Therefore, we gather more information in order to recover better knowledge of the top of the distribution of wealth and better estimate the total population indices. We use bracketed information on components of the total wealth as well as bracketed information involving several components: the total wealth (sum of the components) and the information on the imposition on the Solidarity Tax on Wealth (ISF) obtained by matching with data from the tax declarations.

Our approach relies on Bayesian multiple imputations (Little and Rubin 2002) for sample survey estimation. However, our strategy to produce point estimates or interval estimates is slightly different. We do not rely on particular rules for combining complete-data inferences or rely on imputations which are Bayesianly proper [START_REF] Schafer | Analysis of Incomplete Multivariate Data[END_REF]). We use a hierarchical modeling in order to take into account in the coverage of the interval estimates the uncertainty due to sampling and total-nonresponse and to incomplete knowledge of the censored wealth. Under the assumption that the proper parametric class for the data generating process (DGP) is known, the uncertainty on the value of the censored wealth boils down to uncertainty on the parameters values and to the remaining model uncertainty due to imperfect observation conditional on the knowledge of the parameters. The first model is the model for the quantities interesting to the decision maker [START_REF] Geweke | Contemporary Bayesian Econometrics and Statistics[END_REF] which are here taken to be the various inequality indices on the finite French population. The remaining models, which are standard, are instrumental in the sense that it is not our goal to produce inference on the posterior predictive distribution of the wealth of the sampled households or on the posterior distribution of the parameters, though they could be obtained simultaneously with the numerical procedure. In order to gather information on the components of total wealth, our DGP is multivariate. It allows to account for example for unobserved heterogeneity. From that point of view, our approach is similar to that of [START_REF] Heeringa | Multivariate Imputation of Coarsened Survey Data on Household Wealth[END_REF], where a multivariate model is used in the context of the American Health and Retirement Survey (HRS). Inference is based on a single path of a Gibbs sampler Markov chain that updates the sampled data, parameters and an error term accounting for sampling error. The only mathematical tool is the Ergodic theorem. The model is discussed in Section 2, the non-rectangular censoring is presented in Section 3. In Section 4 we give details on the Gibbs sampler and discuss our strategy to produce point and interval estimates. Results are presented in Section 5 and it is followed by a discussion in Section 6.

The Modeling Assumptions

We denote by U the finite population of size N composed of all the French households and S ⊂ U the sample, where we number the elements of S from 1 to m. Due to total nonresponse, S corresponds to a subset of the initial sample drawn from particular bases of dwellings. The initial sample is stratified and drawn with unequal probabilities where dwellings composed, at the time of the census, of self employed, executives, retired people and of people living in rich neighborhoods have been over-sampled. It implies that probability of selection is related to the wealth but, in principle, the selection is exogeneous. We assume below that the selection mechanism corresponding to the total nonresponse is also exogeneous and that in the models we have included the adequate covariates to be able to ignore the selection mechanism [START_REF] Little | Statistical Analysis with Missing Data[END_REF].

The target quantities or quantities of interest are taken to be inequality indices on the total wealth of the French: the Gini, Theil or Atkinson's indices, quantiles or inter-quantile ratios. They are functions of the distribution of the total wealth t k for households k from 1 to N . Recall that, for example, the Gini is defined by

G = k∈U (2r(k) -1)t k N k∈U t k -1, where r(k) is the rank of t k . A design based estimate is then Ĝ = k∈U (2r(k) -1)w k t k k∈S w k k∈S w k t k -1,
where w k is the weight of household k, r(k) = j∈S w j I {t j ≤ t k } and I{•} denotes the indicator function. In practice, at INSEE, a normal approximation for the design based estimate is used in order to obtain interval estimates. Also, since the variance of the estimate requires in principle as well the data in U \ S, a variance estimate is used. Ĝ being nonlinear in the weights, the estimate of the variance of Ĝ is approximated by that of its linearized version. The variance estimate should take into account the complex design of the sample and the total nonresponse and raking [START_REF] Deville | Generalized Raking Procedures in Survey Sampling[END_REF]). The procedure is well explained in [START_REF] Dell | Mise en oeuvre de calcul de variance par linéarisation[END_REF]. It is however difficult to justify all the approximations rigorously. We do not aim to enter in these details and start off from the approximation:

Ĝ ≈ G + V Ĝ E
where the error term E is a standard centered Gaussian random variable and the variance estimate is denoted by V Ĝ .

Since the total wealth is censored, we are not able to apply the above tools to estimate the target quantities. We rely a priori on a two stage model. But, since for practical issues we have adopted the Bayesian point of view, we have added an additional ladder to the hierarchy of models:

1. The model (I) for the quantities of interest like the Gini, conditional on the relevant data from the sampled household.

2. The model (DGP) for the components of the wealth for sampled households. It is a multivariate model for owned macro-components of the total wealth among: the financial wealth W 1 , the value of the principal dwelling W 2 , of the other real estate including secondary dwellings rented or for leisure and parking lots W 3 , the professional wealth W 4 and the remainder (durable goods, artwork, private collections and jewelry) W 5 , conditional on the value of covariates x l k , k = 1, . . . m, l = 1, . . . , 5 and on the parameters in a certain parametric class of models.

The prior distribution (P) of the parameters Θ of density π(θ).

For simplicity, we assume that every household has some financial wealth (eg. money on a checking account) and some wealth in form of remainder. Therefore, if we make groups according to the type of portfolio in the 5 above components of wealth we have to distinguish 8 groups. We denote by D k = (D l k ) l=1,...,5 the binary vector such that D l k = I W l k > 0 and define the function P which associates to each D k the number i ∈ 1, . . . , 8 of the pattern. In the remaining, we use capital letters for random variables and lowercase letters for realizations. We also use bold characters for vectors.

The model is defined as follows. In the first stage (I) we set, for example for the Gini,

G = Ĝ (t 1 , . . . , t m ) + V Ĝ (t 1 , . . . , t m ) E, E N (0, 1) (2.1)
with the Assumption (A):

E independent of (t 1 , . . . , t m ) (A).
Concerning the model (DGP), we have the following model for pattern i:

   T k = 5 l=1 W l k , log(W l k ) = x l k β l i + U l i , when d l k = 1 and P (d k ) = i U i N (0, Σ i ) (2.2)
where Σ i is of size

p i = 5 l=1 d l k for any k such that P (d k ) = i.
We make the following restriction on the parameters:

Only the coefficient of the constant is group specific (fixed effect), the remaining coefficients of β l i are equal for all i (RP).

For the last model (P), we choose π(θ) proportional to

8 i=1 det (Σ i ) - p i +1 2 . (2.3)
In usual design based inference in survey sampling, the Gini index G has an unknown but fixed value. Hence, Equation (2.1) is not usual since G is now random. In some approaches to survey sampling though, the finite population values correspond to draws in a super-population and it makes sense to assume that the quantities of interest on the finite population are random. It is also usual to revert the Gaussian approximation to obtain interval estimates. We may also think about estimating G in terms of prediction. We present in Table 1 the kind of covariates we have introduced in the model (DGP). Covariates can also improve a priori the coverage of the interval estimates, up to a certain stage since increasing the size of the vector of parameters deteriorates the knowledge on the parameters. The main justification for introducing the covariates is however to justify Assumption (A). Indeed, the survey sample is drawn exogeneously and we have to condition by the corresponding observed covariates in order to estimate the law of the data unconditional of the selection. Total nonresponse is a second stage of selection for which the selection mechanism is unknown and we assume that this mechanism is ignorable (Little and[START_REF] Little | Statistical Analysis with Missing Data[END_REF][START_REF] Gautier | Eléments sur les mécanismes de sélection dans les enquêtes et sur la non-réponse non-ignorable[END_REF]) and that we condition by the adequate covariates to decondition from selection.

The model (DGP), is such that, though we take into account observed heterogeneity in the form of portfolio allocation and through several covariates, there might remain unobserved heterogeneity (eg. in the form of a missing covariate) like the preference for the risk and time that causes the residuals to be dependent. In order to use product specific variables for the principal dwelling we model the value of the good. In contrast in the other models, for which the variables are sums of components collected in the survey, we model the amount of the share that the household possesses and use household specific variables only.

The vector of parameters θ in R d corresponds to the β l i 's and the matrices Σ i where, denoting by dim l the dimension of β i,l ,

d = 5 l=1 (dim l -1) + 8 * 5 + 1 2 5 k=2 k(k + 1).
The prior is a product of usual priors in the context of Gaussian linear models which are limits of normal/inverse-Wishart's (see for example [START_REF] Little | Statistical Analysis with Missing Data[END_REF][START_REF] Schafer | Analysis of Incomplete Multivariate Data[END_REF]. They are often called non-informative. They indeed correspond to a proper objective choice for the prior for the coefficients β l i . The posterior, if the data were observed, is a bona-fide normal/inverse-Wishart probability distribution.

Censoring

In the absence of scalar measurements, intervals are the main information for identification and estimation. As already discussed, we aim to localize as much as possible the missing W k . For that purpose we use two summarizing questions and the information on which household is eligible to the Solidarity Tax on Wealth. The answers to the summarizing questions take the form of brackets for the sum of the collected components of the financial wealth as well as brackets for the total wealth. Recall that the total wealth includes the remainder which is not collected per se in the detailed questionnaire. The information on which household pays the Solidarity Tax on Wealth has been obtained by matching with a data set from the tax department.The condition to pay the Solidarity Tax on Wealth is to have a taxable wealth exceeding 720.000 . This taxable wealth corresponds to a different concept of wealth. Only part of the professional wealth is taken into account. It is possible to deduct the professional wealth used professionally with the exception that if one owns a share in a firm which is too low then it is not deductible. It is possible to have a rebate of 20% on the value of one's principal dwelling. The artworks are not taxed either. Finally, debts are deducted. It is possible to take into account most of the specificities of this tax. For example, by chance, the few households that possessed a share in a firm gave equal upper and lower bounds for its value. However, it is not possible to distinguish the artworks within the remainder. Hence, we produce lower and upper bounds on the taxable wealth which take the form of an extra bracketing condition.

When an household pays the tax, the upper bound for the taxable wealth

W 1 k + 0.8 * W 2 k + W 3 k + min(W 4 k , N DED max,k ) + W 5 k -DEBT k (3.1)
is greater than 720.000 , where N DED max,k is an upper bound of the nondeductible professional wealth obtained using the detailed information and DEBT k the total debts which is deductible. We assume that households always subtract the deductible amounts. When an household does not pay the tax, the lower bound for the taxable wealth

W 1 k + 0.8 * W 2 k + W 3 k + N DED min,k -DEBT k (3.2)
is lower than 720.000 , where N DED min,k is a lower bound of the nondeductible professional wealth obtained using the detailed information.

The above conditions involving several variables allow, by manipulation of lower and upper bounds, to shorten the initial intervals for each component and to obtain intervals for the remainder. The censoring takes form of a hyperrectangle. The final summarizing condition for the total wealth and the eligibility to the Solidarity Tax on Wealth imply censored domains which are subsets of these hyper-rectangles.

Note that further external information has also been used in order to specify upper bounds on the a priori unbounded total wealth. This choice is questionable but these upper bounds are very loose. The motivation is that, since design based inference is the initial goal and each sampled household has a weight that is the inverse of the probability of selection2 and the data is collected once and for all, we want to have reasonable "representativity" of the samples3 . The average weight is around 2.000. Suppose a billionaire is drawn, then it is assumed to represent 2.000 households. In turn, if we were able to draw at random a second time the sample, it is very likely that no billionaire would be drawn. It might result in inequalities which are often too low and sometimes too high. Probably due to the over-sampling of households suspected to be wealthy, an household with a share in a firm of the order of 25.000.000 has been drawn. We have introduced upper bounds on the total wealth based on [START_REF] Cordier | Les inégalités de patrimoine des ménages entre 1992 et 2004[END_REF] and published information on the highest French professional wealth. We have bounded by 50.000.000 the total wealth of the apparently wealthiest household and by 10.000.000 the total wealth of the others. Such restriction might cause under coverage of the interval estimates. Heeringa, Little and Raghunatan (2002) have noted, using a similar modeling and the HRS survey, that introducing such restrictions implies slightly lower means, Ginis or other concentration indices. On the other hand its impact on more robust quantities such as quantiles is minor. [START_REF] Geweke | Contemporary Bayesian Econometrics and Statistics[END_REF] in the non-rectangular context is not standard and also requires proper importance sampling weighting which in the context of our hierarchical modeling does not seem feasible. Accept-reject with instrumental distribution the unconditional distribution is known to be very ineffective, especially as the dimension increases. Robert (1995) for example suggests the use of the Gibbs sampler [START_REF] Arnold | Gibbs sampling[END_REF]). The total procedure is extremely intensive from a computational viewpoint. Moreover, we need to modify slightly the procedure if we want to take into account the uncertainty on the parameters due to the finite distance.

The Gibbs sampler easily adapts to the Bayesian modeling, see for example [START_REF] Mcculloch | An Exact Likelihood Analysis of the Multinomial Probit Model[END_REF] for Bayesian inference for multinomial discrete choice models. It is also popular in missing data problems (Little and[START_REF] Little | Statistical Analysis with Missing Data[END_REF][START_REF] Schafer | Analysis of Incomplete Multivariate Data[END_REF]) and the extension is called data augmentation. The only difference with the usual Gibbs sampler is that the state space of the Markov Chain is augmented in order to include the parameters. Also, in order to infer on the quantities of interest, we augment the state space once more and include the E's. The Gibbs sampler relies on an exhaustive block decomposition of the coordinates of the state space. These blocks are numerated according to a specific order. Starting from an initial value v 0 , the Gibbs sampler simulates a path from a Markov chain (v n ) n≥0 . Given v n , a vector V n+1 decomposed in the above system of blocks is simulated by iteratively updating the blocks and sampling from the distribution of the block conditional on the values at stage n of the future blocks and the value at stage n + 1 of the previously updated blocks. Here V n is taken to be

V n = (Θ, W ′ 1 , . . . , W ′ m , E) ′ .
The sequence is such that we start by updating the covariance matrices, followed by the β l i , then by the components of wealth one by one, household by household and finish with the error term in model (I). The updating for the distribution without truncation is for example explained in [START_REF] Little | Statistical Analysis with Missing Data[END_REF]. Here, we simulate the components of wealth in truncated univariate normals, which is easy and efficient. We update the intervals of truncation for the current variable at each stage of the sequence with the previously simulated components for the same household.

The limit theorems for the Gibbs sampler are given in Tierney ( 2004) and Roberts and Smith (1994). We can also check as in Roberts and Polson (1994), minorizing the transition kernel using that we have introduced upper bounds for the a priori unbounded amounts, that there is uniform exponential L 1 ergodicity. Thus convergence of the laws of the marginals of the Markov chain to the target joint posterior and posterior predictive and distribution 4 of E, which is the invariant probability µ, is very fast. The main result on Markov chains which is useful for the inference here is the ergodic theorem. It states that for g in L 1 (µ), lim

T →∞ 1 T T n=1 g (V n ) = E µ [g (V)] a.s. (4.1)

Posterior Predictions and Posterior Regions

Suppose that the decision maker wants the statistician to give him a single value for each quantity of interest. A natural question is then to ask: "What is the optimal answer one can give?". Once a loss function, say quadratic, is specified, the optimal answer

G = E µ [G (V) |W 1 ∈ D 1 , . . . , W m ∈ D m , X 1 = x 1 , . . . , X m = x m ] , (4.2)
among all answers G * , minimizes the posterior risk

E µ (G * -G (V)) 2 |W 1 ∈ T 1 , . . . , W m ∈ T m , X 1 = x 1 , . . . , X m = x m (4.3)
where the domains T k correspond to the domains of truncation in R p P (D k ) and x k are the matrices of the covariates5 . The expectation giving the posterior prediction in (4.2) could be approximated by a MCMC method [START_REF] Robert | Monte Carlo Statistical Methods[END_REF]): the empirical mean along one path using (4.1). Since T is chosen by the statistician, this approximation of the integral could be as good as one wishes. As usual in MCMC methods we present in Section 5 results with burn-in, i.e. where we have dropped the first B simulations. According to the Cesarò lemma, (4.1) still holds starting the sum from n = B + 1 and replacing T by T -B. Heuristically, it allows to wait for for the chain to stabilize close to the steady state. It only changes the very last decimals here since we have taken T = 20.000, which is very large as attested by Figure 5, and B = 1.000. With the optimality as a goal, simple imputation is no good. The researcher arbitrarily chooses one scenario for G among an infinity of possible scenarii. Also, it is natural that predicting the quantities of interest is different from predicting the unobserved wealth due to the nonlinearity of G in the wealth. Predicting the unobserved wealth is doomed to introduce biases. Suppose that the statistician has convinced the decision maker that it is better to be given an interval estimate. We take in Section 5 symmetric 95% posterior regions but we could have taken highest posterior regions6 or posterior regions of minimal length. The boundaries of the intervals [l, u] can be obtained, without Central Limit Theorem and with a single path, by inverting the functional:

E µ [I {G (V) ∈ [l, u] c }] ≤ 0.05 (4.4)
where the left hand-side could be well approximated for large T by lim Such interval estimates take into account the uncertainty due to sampling and total nonresponse, to imprecise knowledge of the value of the parameters among a parametric class of models due to finite sample size and to the uncertainty due to the imperfect measurement of the components of wealth conditional on the knowledge of the parameters.

T →∞ 1 T T n=1 I {G (v n ) ∈ [l, u] c } (4.

Presentation of the Results

We have applied the methodology of Section 4 to the data set and runned a Gibbs sampler with T = 20.000 and B = 1.000. We have tried to diagnose convergence by plotting the convergence of empirical averages required for the inference. As expected due to exponential ergodicity, the convergence seems to occur very quickly. We have plotted in Figure 1 the convergence of the empirical averages for the Gini. In Table 2 we collect the estimated posterior predictions and confidence regions.

Discussion

Our multivariate model is similar to the model specified in [START_REF] Heeringa | Multivariate Imputation of Coarsened Survey Data on Household Wealth[END_REF]. Al well we consider different models for different types of portfolios. There are still two differences. We introduce much more information on covariates. This information is at least useful for selection issues. We also allow for different covariance matrices in the different groups while they assume the covariance matrices are blocks extracted from a unique matrix which here would be 5 * 5. It seems that it amounts to integrating with respect to components which are not in the portfolio as if they existed but were not observed. We are not able to justify this choice. Also, it is not clear that the posterior of the unique latent covariance matrix is also normal/inverse-Wishart. One relative disadvantage of our approach is that we have many parameters. That is why we have considered only 5 "macro"-components and introduced the restriction on the parameters (RP). It is possible to consider finer decompositions of the total wealth introducing less covariates. Note that it is possible to avoid the introduction of group specific coefficients by specifying one multivariate Tobit model where the residuals of the 5 latent variables are correlated. These latent variables account for both the amount and the decision to invest in the components. We believe that this is a severe restriction. We have also noted that, since the composition of the portfolio is observed, we do not need to model the choice mechanism and can condition on that information. One of the main difficulty not treated in the paper is inconsistency implied by the fact that brackets for the variables are not coherent with the brackets involving several variables. It allowed to detect errors like confusion between old Francs and Euros. Concerning the final question it turns out that it is very little informative on the top of the distribution. This is troublesome for specification issues when we want to use only this last question. It is troublesome as well when we use the detailed components since we have a quite poor information on the remainder. The remainder is a mixture of luxury and durable goods and the bottom of the distribution for which the intervals are more informative is likely to be mainly composed of durable goods. It is thus always important, but difficult due to different selection mechanism especially due to nonresponse and different perception of surveys, to gather information from sources exterior to the survey. It is expected for the future French survey on Wealth to ask for the right to use more auxiliary information from the tax declarations. Also, it is possible that survey sampling use over sampling based on data from these tax declarations. However, it is customary to inform the households that their data will be matched and it sometimes increases the nonresponse rate when it is matched with tax declarations. The threshold for the last bracket for the final question on the total wealth was set to 450.000 in order to be far from the threshold of 720.000 implying eligibility to the Solidarity Tax on Wealth and mitigate nonresponse rates to this question.

Figure 1 :

 1 Figure 1: Convergence of the empirical averages of the Gini, burn-in of the 1.000 first iterations

Figure 2 :

 2 Figure 2: Posterior distribution of the Gini, Theil and Atkinson indices, T=20.000, burn-in of the 1.000 first iterations

Table 1 :

 1 Covariates for the model (DGP) other than the type of portfolio √ √

	Surface and square of the surface	√
	Professional wealth
	wealth used professionally	√
	firm owned	√

  We adopt the Bayesian point of view for practical reasons. From a frequentist point of view, the first stage consists in the estimation of the parameters of the multivariate linear Gaussian model with censored observations. Several methods are at hand among which simulated based methods like the simulated maximum likelihood, simulated scores or simulated method of moments (see[START_REF] Train | Discrete Choice Methods with Simulation[END_REF] for a lively and basic introduction) or a MC-ECM variant of the EM algorithm[START_REF] Little | Statistical Analysis with Missing Data[END_REF]. Also, since the criterion function has in general several local extrema, it is often useful to use a stochastic optimization method or at least to try several initial starting points. It is also useful, as we will see below, to simulate in the multivariate truncated normals in order to finally infer on the finite population inequality indices. The GHK simulator

	4 Numerical Procedure for Point and Interval
	Estimates
	4.1 The Gibbs sampler

  Unlike the original Bayesian multiple imputations (Little and Rubin 2002 and Schafer 2001) we do not require proper Bayesian imputations, i.e. independent sampling, nor rely on formulas to combine multiple imputations. Multiple imputations here is only a tool to infer on the quantities of interest.

	5)
	using (4.1). Again it is possible to use burn-in.
	Remark 4.1

Table 2 :

 2

	Quantity of interest	Prediction	Lower bound Upper bound
	Mean ( )	205.003,98	192.879,50	217.647,24
	Median ( )	111.459,26	105.672,32	117.563,45
	P99 ( )	1.584.602,96 1.359.261,98	1.825.362,43
	P95 ( )	690.793,96	636.924,50	746.759,31
	P90 ( )	434.458,13	416.410,51	452.006,79
	Q3 ( )	232.307,50	224.849,03	240.204,86
	Q1 ( )	16.998,67	15.117,08	19.149,60
	P10 ( )	3.959,07	2.870,83	5.070,55
	P95/D5	6,1972	5,7232	6,6808
	P99/D5	14,2175	12,1667	16,4388
	Q3/Q1	13,6847	12,2838	15,1034
	D9/D1	109,9332	80,7615	140,5827
	D9/D5	3,8981	3,7081	4,0858
	Gini	0,6519	0,6328	0,6717
	Theil	0,9044	0,8138	1,0001
	Atkinson (ǫ = 1.5)	0,9063	0,8838	0,9253
	Atkinson (ǫ = 2)	0,9742	0,9549	0,9920
	We finally represent in Figure 2 histograms for the posterior distribution of some
	of the target quantities of interest.

Posterior predictions and posterior regions (lower and upper bound of a 95% symmetric region), T=20.000, burn-in of the 1.000 first iterations

It is in fact estimated due to total nonresponse.

The samples are simulated, as explained later

Recall that it is always independent of the rest of the components.

For example block diagonal with p P (D k ) rows.

HPD regions.

Acknowledgement. Eric Gautier is grateful to Christian Robert for his guidance on Monte-Carlo Markov chains and his colleagues at the French Statistical Office INSEE among which Marie Cordier, Cédric Houdré, Daniel Verger and Alain Trognon for stimulating discussions on Economics of Inequalities and Econometrics. The author also thanks the participants of the Social Statistics Seminar at INSEE and the Econometrics Research Seminar at Yale.


Hassane Aboua Ïssa
email: hassane.abouaissa@univ-artois.fr

Michel Fliess
email: michel.fliess@polytechnique.edu

Cédric Join
email: cedric.join@cran.uhp-nancy.fr

Fast Parametric Estimation for Macroscopic Traffic Flow Model

Keywords: Traffic control, linear systems, parameter identification

The main traffic parameters such as the critical density and the free flow speed of a motorway segments are subject to changes over time due to traffic conditions (traffic composition, incidents, . . . ) and environmental factors (dense fog, strong wind, snow, . . . ). As such parameters have an impact on the performance of the traffic control strategies, they must be estimated on-line. Our approach, which is of algebraic flavor and avoids asymptotic and statistical techniques, yields fast implementable formulae in closed form. Some convincing computer simulations are provided.

INTRODUCTION

Intelligent transportation systems were introduced to ensure an optimal infrastructure utilization and to improve the driving safety, by alleviating congestion phenomena. All proposed control strategies rest on the use of models particularly the macroscopic one (Papageorgiou, Haj-Salem & Blosseville (1991)). However, some of the important parameters of such a model are usually unknown and are subject to frequent changes over time [START_REF] Ozbay | Improved online estimation methods for a feedback-based freeway ramp metering strategy[END_REF]). The success of traffic control strategies depends therefore on the capacity to calibrate the control parameters as soon as one is able to track changes in the main traffic model characteristics [START_REF] Wang | Real-time parallel parameter estimators for a second-order macroscopic traffic flow model[END_REF]). The model parameters calibration is generally done off-line using simulation software as in [START_REF] Zhang | Design and evaluation of roadway controller for freeway traffic[END_REF] or estimated off-line as in [START_REF] Karaaslan | Two proposals to improve freeway traffic flow[END_REF]. Nevertheless, as the parameters are subject to changes over time owing traffic conditions (traffic composition, incidents, ...) and environmental factors (dense fog, strong wind, snow, ...), their estimation must be achieved on-line. In this context, several parameter estimators are available. [START_REF] Bellemans | Model predictive control with repeated model fitting for ramp metering[END_REF] have proposed a nonlinear optimization algorithm that refits the most sensitive traffic flow model parameters. The approach is based on a regular set of data. Whereas the method is very powerful, it requires more computation, and the fitting with one set of data may not guarantee a good results with another set of data. [START_REF] Ozbay | Improved online estimation methods for a feedback-based freeway ramp metering strategy[END_REF] have used an extended Kalman filter (EKF) to estimate the critical density considered as state variable. [START_REF] Wang | Real-time parallel parameter estimators for a second-order macroscopic traffic flow model[END_REF] have developed a parallel adaptive least-squares parametric estimator for a second-order macroscopic model. Several other works (see, e.g., [START_REF] Wang | Real-time freeway traffic state estimation based on extended kalman filter: a general approach[END_REF], [START_REF] Wang | A real-time freeway network traffic surveillance tool[END_REF]) are utilizing the same principle of EKF to estimate together traffic parameters and states. This paper1 is proposing a new fast parametric estimation method resulting from algebraic technics first introduced in [START_REF] Fliess | An algebraic framework for linear identification[END_REF], 2008) 2 . This new setting, which avoids statistical and asymptotic techniques, yields fast implementable algebraic formulae, in closed form. Excellent estimates of the most sensitive traffic flow model parameters, like the critical density and the free flow speed, are obtained in this way.

After recalling in Section 2 a basic macroscopic model, Section 3 demonstrates that the basic parameters may be estimated thanks to this new setting, which is briefly reviewed. Convincing numerical simulations3 are provided in Section 4. We assume for the dynamics of the traffic flow a first order macroscopic model (LWR) [START_REF] Lighthill | On kinematic waves II: A theory of traffic flow in long crowded roads[END_REF]; [START_REF] Richards | Shock waves on the highway[END_REF]), which is briefly presented below. Consider a freeway stretch subdivided into n segments as shown in Figure 1. The macroscopic first order model is based on the following mass conservation law:

∂ρ ∂t + ∂q ∂x = 0
For a space discrete representation, define the traffic density ρ i as the number of vehicles in segment i at time t, divided by the segment length L i . The traffic volume q i-1 in veh/h, is defined as the number of vehicles entering the segment i; q i is the number of vehicles leaving the segment i. The equation of the nonlinear model for each freeway segment reads:

ρi (t) = 1 L i (q i-1 (t) -q i (t) + a i r i (t) -b i p i (t)) (1)
where, a i and b i are binary variables which indicate respectively the presence or the absence of an on-ramp r i (t) and off-ramp p i (t). Note that another space discrete model [START_REF] Mammar | Coordinated ramp metering via second order sliding mode control[END_REF]) that take into account the shock waves phenomena may also be utilized. Eq. ( 1) is supplemented by the following two equations:

q i (t) = ρ i (t)v i (t) (2) v i (t) = v f,i 1 - ρ i (t) ρ m,i (3)
Here, one considers the Greenshields fundamental diagram [START_REF] Greenshields | A study of traffic capacity[END_REF]) based on Eq. ( 3), where, v f,i is the free flow speed or the speed of the vehicle in the segment i when there are no interactions between vehicles. Write ρ m,i the jam density of the segment i. Note that each segment is characterized by a capacity defined by the maximum flow q m,i corresponding to the critical density ρ cr,i . There is a linear relation between this density and the maximal density ρ m,i : [START_REF] Ozbay | Improved online estimation methods for a feedback-based freeway ramp metering strategy[END_REF], the segment capacity is not constant. It may change over time owing to circumstances such as traffic condition and driver behavior, traffic composition and external factors (environmental conditions and segment geometry, . . . ). Strong wind and dense fog, for example, do also affect the flow-density and the speed-flow relationships. This means that both the critical density and the free flow speed are subject to changes, which may be quite abrupt. Due to their impact on the control strategies performances, these two main parameters ought to be estimated on-line.

ρ cr,i = ρ m,i 2 As stated in

ON-LINE PARAMETER ESTIMATION

Let us start with an elementary academic example for explaining the basic principles of our techniques (see [START_REF] Fliess | An algebraic framework for linear identification[END_REF], 2008) for details).

An elementary academic example

Consider the first-order monovariable input-output system ẏ(t) = Ky(t) + u(t) (4) for t ≥ 0, where:

• u(t) and y(t) are respectively the control and output variables, • the constant gain K is unknown.

Rewrite Eq. (4) in the usual operational form (see, e.g., [START_REF] Van Der Pol | Operational calculus based on the two-sided Laplace integral[END_REF]; [START_REF] Yosida | Operational Calculus[END_REF]), which is most often introduced via the classic Laplace transform: sŷ = K ŷ + û + y(0) Apply to both sides the derivative d ds , which corresponds in the time domain to the multiplication by -t (van der [START_REF] Van Der Pol | Operational calculus based on the two-sided Laplace integral[END_REF]; [START_REF] Yosida | Operational Calculus[END_REF]), in order to eliminate the initial condition y(0), the measurement of which might be corrupted by noises:

ŷ + s dŷ ds = K dŷ ds + dû ds (5)
It yields the linear identifiability [START_REF] Fliess | An algebraic framework for linear identification[END_REF], 2008)) of K:

K = ŷ + s dŷ ds -dû ds dŷ ds
In order to get rid of time derivatives, multiply both sides of Eq. ( 5) by s -N , where N ≥ 1 is large enough. Remember that negative powers of s corresponds in the time-domain to iterated time-integrals (van der [START_REF] Van Der Pol | Operational calculus based on the two-sided Laplace integral[END_REF]; [START_REF] Yosida | Operational Calculus[END_REF]). Those iterated integrals are classic instances of low pass filters, which are attenuating corrupting noises, when considered as highly fluctuating phenomena4 .

Remark 3.1. The above calculations were for simplicity's sake assuming noises of zero mean. This is not the case in the numerical experiments of Section 4, where there is a constant mean equal to 150. Let us indicate very briefly how to modify our previous calculations5 . Replace Eq. ( 4) by ẏ(t) = Ky(t) + u(t) + C where C is an unknown constant. Eq. ( 5) becomes sŷ = K ŷ + û + y(0) + C s . We are eliminating y(0) and C by multiplying both sides of the last equation by the differential operator d 2 ds 2 s. The remaining estimation of K is quite similar to what we have already done. Remark 3.2. An excellent estimated value of K, which is derived via iterated time integrals, may be obtained by utilizing quite short time windows [START_REF] Fliess | An algebraic framework for linear identification[END_REF], 2008)). This is why we may speak of a fast or on-line parametric estimation technique.

Traffic flow parametric estimation

The speed v i is easily obtained via Eq. ( 2) by utilizing the measured signals q i and ρ i = 0. In what follows v i is thus supposed to be known. Set θ 1,i = v f,i and θ 2,i = v f,i 2ρcr,i for estimating the parameters v f,i and ρ cr,i , which are linearly identifiable. Eq. ( 3) reads here:

v i (t) = θ 1,i -θ 2,i ρ i (t)
It becomes in the operational domain:

vi (s) = θ 1,i s -ρi (s)θ 2,i
if the unknown parameters θ 1,i and θ 2,i are assumed to be constant. Multiplying both sides of the preceding equation by Π 1 = 1 s 2 d ds s and Π 2 = 1 s yields: vi (s)

s 2 + 1 s d ds vi (s) = -θ 2,i ( ρi (s) s 2 + 1 s d ds ρi (s)) vi (s) s = θ 1,i s 2 -θ 2,i
ρi (s) s where only iterated integrals appear in order to get a good noise attenuation. The well known Cauchy rule yields in the time-domain the following estimators:

[θ 2,i ] e = - T 0 (T -2τ )v i (τ )dτ T 0 (T -2τ )ρ i (τ )dτ [θ 1,i ] e = 1 T [θ 2,i ] e T 0 ρ(τ )dτ + T 0 v i (τ )dτ
where [0, T ] is the time window, which is necessary for achieving reliable calculations. This time window, which is quite short, permits on-line estimations of

[v f,i ] e = [θ 1,i ] e and [ρ cr,i ] e = [θ1,i]e 2[θ2,i]e .
Remark 3.3. See, for instance, [START_REF] Fliess | An algebraic framework for linear identification[END_REF] when the denominators in the above formulae are close to zeros and yield numerical instabilities.

NUMERICAL EXPERIMENTS

First results

Consider a simple section of 1 km length. The nominal traffic flow changes between 1800 and 1400 veh/h. Set v f = 60 km/h for t < 0.4h v f = 72 km/h for t ≥ 0.4h and ρ cr = 60 veh/km for t < 0.7h ρ cr = 48/km for t ≥ 0.7h Figure 2 demonstrates that, without corrupting noises, good estimates are obtained almost immediately after the parameter variations. A very small integration time window, here 10 samples6 , is sufficient.

The output measurements q in , ρ are now corrupted by an additive non-Gaussian noise, of mean 150, which is fluctuating between 0 and 300. Without pre-filtering7 , this noise is degrading the reconstruction of v. Compare Figures 2-(c) and 3-(c). We have to increase the integration time window from 10T e to 600T e , in order to ensure good and robust estimates.

Application to a ramp metering

The equation for the same section as above, with an onramp, reads:

ρ(t) = 1 L (q in (t) -q out (t) + r(t))
where L and r(t) are the length and the control variable.

Combining it with Eq. ( 2), (3) yields:

ρ(t) = 1 L (q in (t) -ρ(t)[v f ] e 1 - ρ(t) 2[ρ cr ] e ) + r(t) (6)
where [•] e denotes the estimated value. The objective of the control, is to maintain the traffic density below a critical one ρ < ρ cr . Let us define the maximal density ρ m = 120veh/km. For a free-flow speed of v f = 60km/h and a critical density ρ cr = 60, the section capacity expressed in term of traffic flow q m = 1800veh/h. The initial density ρ(0) = 10veh/km, which means that the traffic is in an un-congested regime.

As already noticed by [START_REF] Abouaïssa | Integrated control of highway systems using flatness-based concept[END_REF], Eq. ( 6) defines a controllable linear system which is therefore flat (see [START_REF] Fliess | Flatness and defect of non-linear systems: Introductory theory and examples[END_REF]; [START_REF] Fliess | Continuous-time linear predictive control and flatness: A module-theoretic setting with examples[END_REF]; [START_REF] Rotella | Commande des systèmes par platitude[END_REF]; Sira-Ramírez & Agrawal ( 2004)); y = ρ is a flat output. It yields the closed-loop control

r(t) = L ẏ * -K 1 (y -y * ) + y[v f ] e 1 - y 2[ρ cr ] e -q in
where the gain K 1 is to be tuned. We do not consider the queue length that can be formed at the on-ramp and may provide negative ramp flow values. The results reported in Figure 4 Future works will present more realistic scenarios for traffic control. The new on-line nonlinear estimation and identification techniques [START_REF] Fliess | Non-linear estimation is easy[END_REF]) should play an important role.

Fig. 1 .

 1 Fig. 1. A freeway stretch subdivided into segments 2. TRAFFIC FLOW MODELLING

Fig

  Fig. 2. Simulation results with noise-free measurements

See also[START_REF] Abouaïssa | Modélisation macroscopique du flux de trafic basé sur l'estimation paramétrique des modèles bond graphs[END_REF]. See Villagra, d'Andréa-Novel,[START_REF] Sira-Ramírez | Estimation of longitudinal and lateral vehicle velocities: An algebraic approach[END_REF] for applications to another traffic control problem via analogous methods.

See Fliess, Join & Sira-Ramírez (2004) for applications to linear diagnosis. See[START_REF] Fliess | An introduction to algebraic discretetime linear parametric identification with a concrete application[END_REF] for the corresponding parameter identification in the discrete-time case

.3 Lack of space is preventing us to detail the computer implementation. It may be obtained from H. Abouaïssa (hassane.abouaissa@univ-artois.fr) and C. Join (cedric.join@cran.uhp-nancy.fr).

See Fliess (2006) for more details. Let us emphasize that no statistical tools are needed.

See Fliess (2006) and[START_REF] Fliess | An algebraic framework for linear identification[END_REF], 2008) for more details.

The sampling period Te is equal to 1s.

Filtering those two measurements would increase the estimation performances. We did not do it here in order to test our approach.