Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.
81,446 characters · 16 sections · 32 citation commands
Tests for qualitative features in the random coefficients model
}
\paragraph{Keywords:} Gaussian approximation; mode detection; monotonicity; multiscale statistics; shape constraints; Radon transform; ill-posed problems.
In the random coefficients model, $n$ i.i.d. random vectors $(\mathbf{X}_i,Y_i),$ $i=1,\ldots,n$ are observed, with $\mathbf{X}_i=(X_{i,1}, \ldots, X_{i,d})$ a $d$-dimensional vector of design variables and
The unobserved random coefficients ${\boldsymbol{\beta}}_i=(\beta_{i,1}, \ldots, \beta_{i,d}),$ $i=1,\ldots,n,$ are i.i.d. realizations of an unknown $d$-dimensional distribution $F_{\boldsymbol{\beta}}$ with Lebesgue density $f_{\boldsymbol{\beta}}.$ Design variables and random coefficients are assumed to be independent. The statistical task is to recover properties of the joint density $f_{\boldsymbol{\beta}},$ which is assumed to belong to some nonparametric class. In this work, we derive tests for increases and modes of $f_{\boldsymbol{\beta}}.$
For $d=1,$ the random coefficients model simplifies to nonparametric density estimation. For $d>1,$ recovery of $f_{\boldsymbol{\beta}}$ is an inverse problem with ill-posedness depending on the distribution of the design vectors $\mathbf{X}_i.$ If the design is sufficiently regular, the inverse problem is mildly ill-posed. Otherwise, the model can be severely ill-posed or even be non-identifiable. In this work, we study the mildly ill-posed regime and consider in particular the random coefficients model with random intercept
which can be obtained from (ref) setting $X_{i,1} = 1,$ almost surely.
Random coefficients models appear in econometrics and epidemiology and are used to model unobserved heterogeneity in the population. While the standard linear regression model accounts for unobserved heterogeneity only by an intercept that varies across the population, the random coefficients model allows in addition that different individuals have different slopes. Applications in epidemiology are considered by Greenland00, Greenland06. In economics, random coefficients models are frequently used to evaluate panel data, cf. Hsiao14 or Hasio04, Chapter 6, for an overview. Modeling and estimating consumer demand in industrial organization and marketing often makes use of random coefficients BLP:95, Petrin:02, Nevo:01, Berry07, Dube12. In all these works, parametric assumptions on $f_{{\boldsymbol{\beta}}}$ are imposed. Recently, nonparametric approaches for random coefficients became popular in microeconometrics hoderlein2008, Masten15, HHM:15, DHK:17, frequently combined with binary choice Ichimura98, Gautier12, Gautier13, Masten14, DHK:13, Fox16, DHKS:18, among others.
The random coefficients model also includes quantum homodyne tomography. In this case, we observe an angle $\Phi_i$ and
with $(Q_i,P_i)$ i.i.d. random variables which are unobserved and independent of $\Phi_i.$ The angles $\Phi_i$ can be chosen by the experimenter and are typically uniform on $[0,\pi].$ The interest is in reconstruction of the Wigner function which takes the role of the joint density of $(Q_i,P_i).$ Because $P_i$ and $Q_i$ are not jointly observable, the Wigner function can take negative values. For more on quantum homodyne tomography and the Wigner function, see butucea2007.
We propose a nonparametric test for shape information of the joint density $f_{\boldsymbol{\beta}}$ in the random coefficients model. The focus will be on a test for directional derivatives and modes. The nonparametric estimation theory for $f_{{\boldsymbol{\beta}}}$ has been developed in Beran1992,Beran1996,Feuerverger2000,hoderlein2008. Due to the ill-posedness of the problem and the curse of dimensionality induced by $d,$ pointwise estimation rates are slow. The reason is that small perturbations in the signal are indistinguishable given the data. Nevertheless, we can get good detection rates for larger features, such as an accentuated mode or a strong increase in the joint density along some direction. From a practical point of view, the relevant information regarding an unknown density is typically its shape rater than its precise, full reconstruction. It is therefore essential to recover increases/decreases and the modes of a density. If, say, two modes in the joint density of two random quantities are detected, this indicates that two different groups can be identified. Hence, shape information allows to interpret a given dataset.
Larger features of the density will also be discovered by a nonparametric estimator even if it suffers from slow pointwise convergence. There are, however, two important reasons why a testing approach might be more appropriate. Firstly, with a significance test of level $\alpha$ we can conclude that with probability $1-\alpha$ a detected feature is not an artifact. Secondly, for an estimator we need to pick one bandwidth or smoothing parameter while detection of different features might require different bandwidth choices depending on the size of the hidden features themselves. Indeed, a short and steep increase will be best detected on a small scale whereas for finding a longer and less strong increase the choice of a larger bandwidth is beneficial. Using multiple testing methods, it is possible to combine a whole range of smoothness parameters into one test and to adapt to different shapes of features.
We construct a so called multiscale test, aggregating single tests on different scales and directions. Multiscale tests can be viewed as a multiple testing procedure specifically designed for nonparametric models. Given a model, the theoretical challenge is to prove that a multiscale statistic can be approximated by a distribution free statistic which is independent of the observations. This allows us then to compute quantiles and to find approximations for the critical values of the multiscale statistic. So far, qualitative feature detection based on multiscale statistics has been studied for various nonparametric models, including the Gaussian white noise model duembgen2001, density estimation duembgen2008 and deconvolution schmidthieber13. In multivariate settings the classical KMT approximation suffers from the curse of dimensionality which then leads to very restrictive conditions on the usable scales. Instead, very recent results on Gaussian approximations of suprema of multivariate empirical processes developed by Chernozhukov2017 can be used Eckle,proksch16. In this work, we extend these techniques. The main difficulties are twofold. First, we need to derive specific properties of the inverse Radon transform for general dimension $d$. Second, in contrast to the other works on multiscale inference, no distribution free approximation can be obtained and we therefore need to study the approximating process if several unobserved functions are replaced by estimators.
In order to study the power of the multiscale test, a theoretical detection bound and numerical simulations are provided. The theoretical result gives conditions under which a mode can be detected. In a numerical simulation study, we investigate the power of the test for increases/decreases along some direction and mode detection in dependence on the sample size and the design variables. We also analyze real consumer demand data from the British Family Expenditure Survey.
Let us briefly summarize related literature on testing in the random coefficients model. Under a parametric assumption on the density $f_{\boldsymbol{\beta}},$ Beran1993 considers goodness-of-fit testing and swamy1970, andrews2001 test whether some of the random coefficients are deterministic. The only test based on a nonparametric assumption was proposed recently by Breunig16. It allows to assess whether a given set of data follows the random coefficients model.
This paper is organized as follows. In Section (ref), we describe the connection between the random coefficients model and the Radon transform. Rewriting the model as an inverse problem in terms of the Radon transform reveals the ill-posed nature of the model. This allows us to construct and to analyze the multiscale test in Section (ref). In this part we also derive the asymptotic theory of the estimator and obtain theoretical detection bounds. In Section (ref) the test is studied for simulated data. As a real data example, consumer demand is analyzed in Section (ref). Proofs and technicalities are deferred to a supplement. An R package and Python code is available as a supplement as well, \url{https://arxiv.org/abs/1704.01066}.
{\it Notation:} Throughout the paper, vectors are displayed by bold letters, e.g. $\mathbf{X}, {\boldsymbol{\beta}}.$ Inequalities between vectors are understood componentwise. The Euclidean norm on $\mathbb{R}^d$ is denoted by $\|\cdot\|$ and the corresponding standard inner product by $\langle \cdot, \cdot \rangle.$ We further denote by $\mathbf{e}_1,\ldots,\mathbf{e}_d\in\mathbb{R}^d$ the standard ON-basis of the $d$-dimensional Euclidean space, $\mathbb{S}^{d-1}$ denotes the unit sphere in $\mathbb{R}^d$ and we write $\mathcal{Z}$ for the cylinder $\mathcal{Z}=\mathbb{R}\times \mathbb{S}^{d-1}$. Furthermore, we write $\mathbf{v}$ for any direction $\mathbf{v}=\sum_{j=1}^d\mathrm{v}_j\mathbf{e}_j\in \mathbb{S}^{d-1}$. For two positive sequences $(a_n)_n,$ $(b_n)_n,$ $a_n \lesssim b_n$ or $b_n\gtrsim a_n$ mean that for some positive constant $C,$ $a_n \leq Cb_n$ for all $n.$ As usual, we write $a_n \asymp b_n$ if $a_n \lesssim b_n$ and $b_n \lesssim a_n.$
and $\mu_{d-1}$ the surface measure on the $(d-1)$-dimensional hyperplane $\{\mathbf{b}\in\mathbb{R}^d:\langle\mathbf{b},{\boldsymbol{\theta}}\rangle=s\}$. The Radon transform maps therefore a function to all its integrals over hyperplanes parametrized by $(s,{\boldsymbol{\theta}})\in\mathcal{Z}.$ The figure above shows the parametrization in two dimensions.
For the connection between the Radon transform and the random coefficients model (ref) consider the normalized observations
The random vectors ${\boldsymbol{\Theta}}_i$ take values in the $(d-1)$-dimensional sphere $\mathbb{S}^{d-1}$. In the random coefficients model with intercept (ref), ${\boldsymbol{\Theta}}_i$ is always in the upper hemisphere, i.e. the first component of ${\boldsymbol{\Theta}}_i$ is positive. In this case, we extend the distribution of ${\boldsymbol{\Theta}}_i$ to the whole sphere by randomizing the signs of the design variables. For this purpose, we generate independent random variables $\zeta_i$, $i=1,\hdots,n$, with $\mathbb{P}(\zeta_i=1) =\mathbb{P}(\zeta_i=-1) =1/2,$ which are independent of the data $(\mathbf{X}_i, Y_i),$ $i=1,\ldots,n$, and define $S_i := \zeta_i Y_i/\|\mathbf{X}_i\|$ and ${\boldsymbol{\Theta}}_i := \zeta_i \mathbf{X}_i/\|\mathbf{X}_i\|.$ Independent of the symmetrization, we have $S_i = \langle {\boldsymbol{\Theta}}_i , {\boldsymbol{\beta}}_i \rangle$. The conditional distribution of $S_1 | {\boldsymbol{\Theta}}_1$ is therefore \[F_{S|{\boldsymbol{\Theta}}}(x|{\boldsymbol{\theta}})=\mathbb{P}(S_1\leq x|{\boldsymbol{\Theta}}_1={\boldsymbol{\theta}})=\mathbb{P}(\langle{\boldsymbol{\Theta}}_1,{\boldsymbol{\beta}}_1\rangle\leq x| {\boldsymbol{\Theta}}_1={\boldsymbol{\theta}})=\int_{-\infty}^{x}(Rf_{{\boldsymbol{\beta}}})(s, {\boldsymbol{\theta}}) \,ds, \] and the conditional density becomes
Recall that we have access to an i.i.d. sample $(S_i, {\boldsymbol{\Theta}}_i)$ of the joint density $f_{S, {\boldsymbol{\Theta}}}$. This allows for nonparametric estimation of $f_{S|{\boldsymbol{\Theta}}}$ via $f_{S|{\boldsymbol{\Theta}}}(s | {\boldsymbol{\theta}}) = f_{S, {\boldsymbol{\Theta}}}(s , {\boldsymbol{\theta}})/f_{\boldsymbol{\Theta}}({\boldsymbol{\theta}})$. Applying the inverse Radon transform to this estimate gives an estimator for the joint density $f_{\boldsymbol{\beta}}.$ This inversion scheme suffers from two sources of ill-posedness. Firstly, dividing by $f_{\boldsymbol{\Theta}}$ might result in very unstable reconstructions if $f_{\boldsymbol{\Theta}}$ is small. This happens if the normalized design variables ${\boldsymbol{\Theta}}_i =\mathbf{X}_i/\|\mathbf{X}_i\|$ systematically miss observations from some directions. In this case the problem becomes unevenly harder and only logarithmic convergence rates can be obtained Davison1983,Frikel2013,Hohmann2016. When the support of ${\boldsymbol{\Theta}}_i$ does not contain an open ball, $f_{\boldsymbol{\beta}}$ might be non-identifiable. Secondly, even with regularity on the distribution of the design, the Radon inversion is known to be an ill-posed problem with degree of ill-posedness $(d-1)/2$. Hence, regularization of the inversion scheme is necessary.
In this work, we study the mildly ill-posed case where the random directions ${\boldsymbol{\Theta}}_i,$ $i=1,\ldots,n,$ are sufficiently regularly distributed over the sphere and the ill-posedness is only due to the inversion of the Radon transform. The precise assumptions on the design are stated in Section (ref).
Our approach makes use of the following explicit inversion formula of the Radon transform. Define the operator $\Lambda$ via
where $\mathcal{H}_d$ denotes the identity for $d$ odd and the Hilbert transform \[ \mathcal H_d f(u)= \frac{1}{\pi}\lim_{\epsilon\rightarrow 0^+} \int_{(-\infty, u-\epsilon] \cup [u+\epsilon, \infty) } \frac{f(s)}{u-s} ds =\frac{1}{\pi} \text{p.v.} \int_{-\infty}^{\infty} \frac{f(s)}{u-s} ds\quad\big(u\in\mathbb R\big) \] for $d$ even. Let $c_d^{-1}=(-1)^{({d-1})/{2}} 2^{-d}\pi^{{1-d}}$ for $d$ odd and ${c}_d^{-1}= -(-1)^{{d}/{2}}2^{-d}\pi^{{1-d}}$ for $d$ even. If $\varphi$ is a Schwartz function on $\mathbb{R}^d$, then we have the inversion formula
cf. Theorem 3.8 in Helgason2011. The so called back projection operator $R^*$ is the adjoint of the Radon transform with respect to the $L^2$ scalar product. Notice that our constant $c_d$ differs from the constant in Helgason2011 as we use the standard definition of the Hilbert transform and define $R^*$ as the adjoint of the Radon transform (as opposed to the dual transform).
The goal of this work is to derive confidence statements for qualitative features of the joint density of the random coefficients. In particular, we are interested in the detection of modes (local maxima) of the density. Following the approach of schmidthieber13, we express the features in terms of differential operators. To be precise, for a collection of compactly supported, non-negative and sufficiently smooth test functions $\phi_{\mathbf{t},h}$ consider the integral
for some directional vector $\mathbf{v}=(\rm v_1,\hdots,\rm v_d)^\top$ and $d\geq2$. Since there should not be any favored direction, we consider in the following radially symmetric test functions,
with a non-negative and sufficiently smooth kernel $\phi:[0,\infty)\rightarrow [0,\infty)$ with $\int_0^\infty \phi(u) du=1$ and support on $[0,1]$. Moreover, $\operatorname{Vol}(\mathbb{S}^{d-2})$ denotes the volume of the sphere $\mathbb{S}^{d-2}\subset \mathbb{R}^{d-1}$ and $\operatorname{Vol}(\mathbb{S}^0):=2$. Notice that $\phi_{\mathbf{t},h}$ is supported on the ball $B_h(\mathbf{t})$ with center $\mathbf{t}$ and radius $h.$ The normalization for $\phi_{\mathbf{t},h}$ turns out to be convenient but does not entail that $\phi_{\mathbf{t},h}$ integrates to one.
If the integral (ref) is positive, there exists a subset of $B_h(\mathbf{t})$ with positive Lebesgue measure on which $\partial_\mathbf{v} f_{{\boldsymbol{\beta}}}$ is positive. On this subset, $f_{{\boldsymbol{\beta}}}$ is thus strictly increasing in direction $\mathbf{v}.$ Similarly, we can recover a decrease if the integral (ref) is negative. To construct a statistical test for increases and decreases it is therefore natural to use an empirical counterpart of the functional defined in (ref).
Let $\mathcal{T} = \{(\mathbf{t},h,\mathbf{v}) : h\in(0,1], \mathbf a_1+h\leq \mathbf{t} \leq \mathbf{a}_2-h,\mathbf{v}\in \mathbb{S}^{d-1}\}$ for $\mathbf{a}_1\leq\mathbf{a}_2\in\mathbb{R}^d$, where the inequalities for the vectors $\mathbf a_1,\mathbf a_2,\mathbf{t}$ are understood componentwise. For statistical inference regarding the sign of the directional derivatives of $f_{\boldsymbol{\beta}}$, we fix a subset $\mathcal{T}_{n}\subset\mathcal{T}$ and test for all $(\mathbf{t},h,\mathbf{v})\in\mathcal {T}_n$ simultaneously the corresponding hypotheses of the form
and
For constructing global tests, we can now argue as in Eckle. Our main interest are the following three global testing problems (i)-(iii).
{\bf (i) Testing for the presence of a mode at a fixed location.} Tests for the hypotheses (ref) and (ref) can be used for the detection of specific shape constraints such as a mode at a given point $\mathbf{b}_0\in\mathbb{R}^d$. For this purpose, we consider several bandwidths/scales $h$ and for each $h$ consider pairs $(\mathbf{t}_1,\mathbf{v}_1),\hdots,(\mathbf{t}_p,\mathbf{v}_p)$, where $\mathbf{v}_j,$ $j=1,\hdots,p$, are directional vectors and the test locations $\mathbf{t}_j$ are points on the line $\{\mathbf{b}_0+r\mathbf{v}_j:r\geq h\}$ $(j=1,\hdots,p)$ in a neighborhood of $\mathbf{b}_0$. Inference for the presence of a mode at the point $\mathbf{b}_0$ can now be conducted by studying the testing problem
with $h$ ranging over all chosen scales and $j=1,\hdots,p.$ Level and power of the mode test (ref) for different designs are reported in Section (ref). In Section (ref), we also show that it is essential to include several bandwidths/scales $h$ in order to separate modes which are close.
{\bf (ii) A global testing procedure for all modes.} Simultaneous tests for the hypotheses (ref) and (ref) can be used for a global testing procedure to detect all modes of the density on a domain. Compared to the previous case, we search for evidence over a range of different $\mathbf{b}_0$ which then inflates the number of local tests.
{\bf (iii) A graphical representation of the local monotonicity behavior for bivariate densities.} Let $d=2$ and define a subset $\mathcal {T}_n=\{(\widetilde\mathbf{t}_j,h_0,\widetilde\mathbf{v}_j):j=1,\hdots,p\}$ for a fixed scale $h_0$ of the form $\mathcal {T}_n=\mathcal {T}_{\mathbf{t}}\times\{h_0\}\times\mathcal {T}_{\mathbf{v}}$, where $\mathcal {T}_{\mathbf{t}}$ contains the $p/|\mathcal {T}_{\mathbf{v}}|$ vertices of an equidistant grid of width $2h_0$ and $\mathcal {T}_{\mathbf{v}}$ contains the directions. We restrict the testing procedure to one fixed scale for an easy to read graphical representation. We consider four equidistant directions on $\mathbb S^1$ given by $\mathcal {T}_{\mathbf{v}}=\{\mathbf{v}_1,-\mathbf{v}_1,\mathbf{v}_2,-\mathbf{v}_2\}.$ Since $\mathcal {T}_{\mathbf{v}}=-\mathcal {T}_{\mathbf{v}}$, we have symmetry in the hypotheses, i.e. $H_{0,+}^{\widetilde\mathbf{t}_j,h_0,\widetilde\mathbf{v}_j}=H_{0,-}^{\widetilde\mathbf{t}_j,h_0,-\widetilde\mathbf{v}_j}$. Therefore, we test only $H_{0,-}^{\widetilde\mathbf{t}_j,h_0,\widetilde\mathbf{v}_j}$ for all triples $(\widetilde\mathbf{t}_j,h_0,\widetilde\mathbf{v}_j)\in\mathcal {T}_n.$ Figure (ref) displays an example for the test outcome with the hypotheses in (ref) and $\mathcal T_n$ as above. An arrow in a direction $\widetilde\mathbf{v}_j$ at a location $\widetilde\mathbf{t}_j$ represents a rejection of the corresponding hypothesis $H_{0,-}^{\widetilde\mathbf{t}_j,h_0,\widetilde\mathbf{v}_j}$ and provides an indication of a negative directional derivative of $f_{\boldsymbol{\beta}}$ in direction $\widetilde\mathbf{v}_j$ at the location $\widetilde\mathbf{t}_j$. Thus, Figure (ref) provides strong evidence that the density is trimodal with modes close to the locations $(-0.5,-0.5)^\top$, $(1.5,-0.5)^\top$, and $(0.5,1.5)^\top$. A detailed description of the settings used to generate Figure (ref) and an analysis of the results is given in Section (ref).
We now derive an empirical counterpart of the functional (ref) in terms of the Radon transform $Rf_{\boldsymbol{\beta}}$. We make the following assumptions for the inversion of the Radon transform.
\newenvironment{thmbis}[1] { \addtocounter{assump}{-1}
}
Assumption (ref) is too restrictive for quantum homodyne tomography (model (ref)), where the density $f_{\boldsymbol{\beta}}$ is given by the Wigner function. The Wigner function can take negative values and is not compactly supported. In this case, we replace Assumption (ref) by the following conditions.
Under the assumptions above the inversion formula (ref) holds for partial derivatives of the test functions $\partial_{\mathbf{v}}\phi_{\mathbf{t},h}$. This is a direct consequence of Theorem 3.8 in Helgason2011. The following lemma analyzes the structure of a partial derivative of the test function transformed by the operator $\Lambda$ introduced in (ref) and how this transform depends on $h$.
For a given triple $(\mathbf{t},h,\mathbf{v})\in\mathcal T$ we study the statistic
By Lemma (ref), the expectation of this statistic can be written as
with $d{\boldsymbol{\theta}}$ being the surface measure on $\mathbb{S}^{d-1}$, i.e. $|S_{\boldsymbol{\Theta}}|=\int_{ S_{\boldsymbol{\Theta}}}\,d{\boldsymbol{\theta}}$ for any measurable $S_{\boldsymbol{\Theta}}\subseteq \mathbb{S}^{d-1}.$ By an application of the inversion formula introduced in (ref) and Lemma 5.1 in Helgason2011, we obtain
Up to rescaling, $T_{\mathbf{t},h,\mathbf{v}}$ is thus an empirical counterpart of the functional defined in (ref).
The statistic $T_{\mathbf{t},h,\mathbf{v}}$ depends on the density $f_{\boldsymbol{\Theta}}.$ In quantum homodyne tomography this density is known. For many other applications, however, $f_{\boldsymbol{\Theta}}$ needs to be estimated from the data. In this case, we use a standard cut-off kernel density estimator $\widetilde f_{{\boldsymbol{\Theta}}}$ for $f_{\boldsymbol{\Theta}}$ based on an additional sample $(S_i,{\boldsymbol{\Theta}}_i),$ $i=n+1,\hdots,2n$ which is independent of $(S_i,{\boldsymbol{\Theta}}_i),$ $i=1,\hdots,n$. See also Appendix (ref) for the definition of the estimator. We replace $f_{\boldsymbol{\Theta}}$ by its estimator and consider the test statistic
As mentioned in Section (ref), the inverse problem might become severely ill-posed or non-identifiable if the density $f_{\boldsymbol{\Theta}}$ approaches zero for some directions. This section provides conditions on the design which ensure that $f_{\boldsymbol{\Theta}}$ has positive H\"older smoothness and is bounded from below and above. These results are of independent interest.
In the random coefficients model (ref), the density $f_{\boldsymbol{\Theta}}$ can be expressed in terms of the density $f_\mathbf{X}$ via $f_{\boldsymbol{\Theta}}({\boldsymbol{\theta}})=\int_{0}^\infty r^{d-1}f_\mathbf{X}(r{\boldsymbol{\theta}})dr$. To enforce that $f_{\boldsymbol{\Theta}}$ is bounded from below we restrict ourselves to designs where $\int_{0}^\infty r^{d-1}f_\mathbf{X}(r{\boldsymbol{\theta}})dr$ is bounded away from $0$. The formula also allows to relate the smoothness of $f_\mathbf{X}$ to the smoothness of $f_{\boldsymbol{\Theta}}.$
Although, the random coefficients model with intercept (ref) could be viewed as a special case of the more general model (ref) it requires a different set of assumptions. For model (ref), we write $f_\mathbf{X}$ as a function of $\boldsymbol x =(x_2, \ldots, x_d) \in \mathbb{R}^{d-1}$ and obtain
see Appendix (ref) for a proof. A necessary condition to ensure that $\inf_{{\boldsymbol{\theta}}} f_{\boldsymbol{\Theta}}({\boldsymbol{\theta}}) >0,$ is given by $f_\mathbf{X}(\boldsymbol x) \gtrsim \|\boldsymbol x\|^{-d}$ as $\|\boldsymbol x\|\rightarrow \infty$. This corresponds to Cauchy-type tails of the design variables. Thinner tails will increase the ill-posedness of the problem. In order to avoid very technical proofs, we consider in the random coefficients model with intercept only the case where $(X_{i,2}, \ldots, X_{i,d})$ follows a multivariate Cauchy distribution, i.\,e.,
with $ \boldsymbol\mu\in\mathbb{R}^{d-1}$ and $\Sigma \in\mathbb{R}^{(d-1)\times (d-1)}$ a symmetric and positive definite matrix. We can compute $f_{\boldsymbol{\Theta}}$ explicitly using (ref)
where $\mathrm{sgn}(\cdot)$ denotes the signum function. In this case, $f_{\boldsymbol{\Theta}}$ is bounded from above and below and is continuously differentiable on the hemispheres $\mathbb{S}^{d-1}_+:= \{{\boldsymbol{\theta}} \in \mathbb{S}^{d-1}\mid\theta_1 > 0\}$ and $\mathbb{S}^{d-1}_- := \{{\boldsymbol{\theta}} \in \mathbb{S}^{d-1}\mid\theta_1 < 0\}$. In particular, if $(X_{i,2}, \ldots, X_{i,d})$ is standard Cauchy, then the density $f_{\boldsymbol{\Theta}}$ is constant. This leads to the following assumptions on the design.
In quantum homodyne tomography we set a global $\gamma$ equal to the minimum of $\gamma'$ from Assumption 2' and $\gamma$ from Assumption (ref).
It is important to notice that statistical testing in the random coefficients model relies on two unrelated sets of assumptions. Firstly, there are assumptions on the density of the random coefficients as introduced in Section (ref). Restrictions of this kind are common in statistical inference for an unknown density. On the other hand, there are assumptions (see Assumption (ref) above) on the design. These assumptions control the ill-posedness of the problem.
Note that Assumption (ref) (ii) can be weakened to $\int_{-\infty}^\infty r^{d-1}f_\mathbf{X}(r{\boldsymbol{\theta}})dr\geq c>0 \quad \text{for all }{\boldsymbol{\theta}}\in\mathbb{S}^{d-1}$. For example, if the support of $\Theta$ is a hemisphere this condition can hold while Assumption (ref) (ii) is violated. This relaxation can be achieved by multiplying independently generated $\zeta_i$ to $\Theta_i$ and $S_i$ as proposed for model (ref) in Section (ref).
This section presents the main theoretical result of the paper stating that the standardized and properly calibrated test statistic (ref) can be uniformly approximated by a maximum of a Gaussian process. For that we need the definition of a Gaussian process on the cylinder $\mathcal{Z}$. To this end, let $\mathcal{B}(\mathcal{Z})$ be the Borel $\sigma$-algebra on $\mathcal{Z}$. Define the $\sigma$-finite measure
Let $\bigl(\mathcal{B}(\mathcal{Z})\bigr)_{\nu}$ denote the collection of all sets of finite $\nu$-measure and let $W$ denote Gaussian $\nu$-noise on $\bigl(\mathcal{B}(\mathcal{Z}),\nu\bigr)$. For disjoint sets $E_1,E_2 \in\bigl(\mathcal{B}(\mathcal{Z})\bigr)_{\nu}$ this implies
Adler2007. $W$ is a random, finitely additive, signed measure. Integration w.r.t. $W$ can be defined similarly to Lebesgue-integration, starting with a definition for simple functions and an extension to general $f\in L^2(\nu)$ via approximation by simple functions in the $L^2$-limit. Integration with respect to $W$ yields
and $\mathrm{Cov}(W(f),W(g))=\langle W(f),W(g)\rangle_{L^2(\mathbb{P})}=\langle f,g\rangle_{L^2(\nu)}$ for $f,g\in L^2(\nu),$ where $L^k(\mathbb{P})$ denotes the collection of all random variables whose first $k$ absolute moments exist. For more details, cf. Adler2007, Chapter 5.2.
Let us provide some heuristic for the Gaussian approximation of $T_{\mathbf{t},h,\mathbf{v}}$. The process $(\mathbf{t},h,\mathbf{v}) \mapsto \sqrt{n}(T_{\mathbf{t},h,\mathbf{v}}-\mathbb{E}[T_{\mathbf{t},h,\mathbf{v}}])$ has in the important case $\mathbb{E}[T_{\mathbf{t},h,\mathbf{v}}]=0$ the same mean and covariance structure as the Gaussian process
In the proof of Theorem (ref) below we show that the expectation $\mathbb{E}[T_{\mathbf{t},h,\mathbf{v}}]$ is asymptotically negligible in the limit process. The test statistic and the Gaussian process depend, however, on the unknown densities $f_{S,{\boldsymbol{\Theta}}}$ and $f_{{\boldsymbol{\Theta}}}$ which have to be estimated from the second part of the sample. To this end, we use the standard cut-off kernel density estimates $ \widetilde f_{{\boldsymbol{\Theta}}}$ and $\widetilde f_{S,{\boldsymbol{\Theta}}}$ defined in Appendix (ref). For Gaussian $\nu$-noise $W$ that is independent of the data let
and
The Gaussian approximation result for the family of test statistics $\widehat{T}_{\mathbf{t},h,\mathbf{v}}$ holds for a finite subset $\mathcal{T}_{n} \subset \mathcal{T}$. Its cardinality may, however, grow polynomially of arbitrary degree with the sample size. Moreover, the range of bandwidths must be bounded from above and below by $h_{\max}$ and $h_{\min}$, both converging to zero as $n$ goes to infinity. The precise conditions are summarized in the following assumption.
Let $\mathcal {A}_p$ be the set of half-open hyperrectangles in $\mathbb{R}^p$, i.e. every $A\in\mathcal {A}_p$ has the representation $A=\{\boldsymbol x\in\mathbb{R}^p:-\infty< \boldsymbol x\leq\boldsymbol a\}$ for some $\boldsymbol a\in\mathbb{R}^p$. For finite sets $S_n$ and two stochastic processes $(X_{s,n})_{s\in S_n}$ and $(\widetilde X_{s,n})_{s\in S_n},$ which are defined on the same probability space, we write
if $\lim_n \sup_{A \in \mathcal {A}_{|S_n|}} | \mathbb{P}( (X_{s,n})_{s\in S_n} \in A ) - \mathbb{P}( (\widetilde X_{s,n})_{s\in S_n} \in A )| = 0.$
With the previous theorem, we can now construct simultaneous statistical tests for the hypotheses (ref) and (ref). If the constant $c_d$ is positive then the method consists of rejecting the hypotheses $H_{0,+}^{\mathbf{t},h,\mathbf{v}}$ in (ref) for small values of $\widehat{T}_{\mathbf{t},h,\mathbf{v}}$ and rejecting $H_{0,-}^{\mathbf{t},h,\mathbf{v}}$ in (ref) for large values of $\widehat{T}_{\mathbf{t},h,\mathbf{v}}$, and vice versa if $c_d$ is negative. Theorem (ref) is used to control the multiple level of the tests. Let $\alpha\in(0,1)$ and denote by $\kappa_n(\alpha)$ the smallest number such that \[ \mathbb{P}\left(\sup_{(\mathbf{t},h,\mathbf{v})\in\mathcal T_n}\beta_h\bigg( \frac{|\widehat X_{\mathbf{t},h,\mathbf{v}}|}{\widehat{\sigma}_{\mathbf{t},h,\mathbf{v}}} - \alpha_h\bigg) \leq\kappa_n(\alpha)\right)\geq 1-\alpha. \] By Theorem (ref), $\kappa_n(\alpha)$ is bounded uniformly with respect to $n$. Define for $(\mathbf{t},h,\mathbf{v})\in\mathcal T_n$ the quantiles
and reject the hypothesis (ref), if
Similarly, the hypothesis (ref) is rejected, whenever
Based on the previous result, we now propose a method for the detection and localization of modes on a subdomain. For convenience, we only consider the case of a hyperrectangle $[\mathbf a_1,\mathbf a_2] \subset \mathbb{R}^d$. We study the case that there is a mode $\mathbf{b}_0$ in the interior $(\mathbf a_1,\mathbf a_2).$ For the multiscale test to have power we need that the set of local tests is rich enough. This can be expressed in terms of conditions on $\mathcal {T}_n.$ Let $H(\mathcal {T}_n)$ be the set of all scales/bandwidths such that for every scale $h$ in this set all triples $(\mathbf{t},h,\mathbf{v})$ are in $\mathcal {T}_n,$ where $\mathbf{t}$ ranges over all grid points of an equidistant grid with component wise mesh size $h$ in the hyperrectangle $[\mathbf a_1+h,\mathbf a_2-h]$, and $\mathbf{v}$ ranges over all grid points of a grid of $S^{d-1}$ and grid width converging to zero with increasing sample size.
The testing procedure is as follows. For any $\mathbf{b}_0$ in $(\mathbf a_1,\mathbf a_2),$ let $\mathcal {T}_n^{\mathbf{b}_0}$ be the set of all sequences of triples $(\mathbf{t}_n,h_n,\mathbf{v}_n)_n \in\mathcal {T}_n$ such that $ch_n\geq\|\mathbf{b}_0-\mathbf{t}_n\|\geq2h_n$ for some sufficiently large $c>2$ and $\angle (\mathbf{t}_n-\mathbf{b}_0,\mathbf{v}_n) \rightarrow 0$ as $n\rightarrow\infty$, where $\angle$ denotes the angle between two vectors. The previous conditions ensure that several such sequences can be found. If for all triples in $\mathcal {T}_n^{\mathbf{b}_0}$ all local tests (ref) reject the hypotheses (ref), we have evidence for the existence of a mode at the point $\mathbf{b}_0$. By choosing the test locations as the vertices of an equidistant grid no prior knowledge about the location of $\mathbf{b}_0$ has to be assumed. Theorem (ref) below states that the procedure detects all modes of the density with probability converging to one as $n\rightarrow\infty$.
The rate for the localization of the modes is $n^{-1/(2d+3)}$ (up to some logarithmic factor). Since $2d+3=d+2\tfrac{d-1}2+4,$ this rate matches the optimal rate for mode detection in an inverse problem with ill-posedness of degree $(d-1)/2$ over a 2-H\"older class. Assumption (ref) requires, however, $h_{\, \textnormal{max}}\lesssim \log(n)^{-14\gamma/(d-1)-5}n^{2\gamma/(d-1)-1}.$ To be able to include scales of the order $n^{-1/(2d+3)}$ we need $\gamma>(d^2-1)/(2d+3).$ The right hand side is smaller than one for $d=2,3.$
\setcounter{equation}{0} In this section we illustrate the finite sample properties of the proposed test in a bivariate and a trivariate setting. In the bivariate setting we illustrate how simultaneous tests for the hypotheses (ref) and (ref) can be used to obtain a graphical representation of the local monotonicity properties of the density. In the trivariate setting we investigate the performance of the test for modality at a given point $\mathbf{b}_0$ (see the hypotheses in (ref)) and the dependence of its power on the distribution of $\mathbf{X}$.
As test function we consider the simplest polynomial which satisfies the conditions of Assumption (ref) for $d=2,3$, that is, \[ \phi(x)=c(56x^3+21x^2+6x+1)(1-x)^6\mathbbm{1}\{x\leq 1\},\quad x\in[0,\infty), \] with $c$ such that $\int \phi =1.$ Figure (ref) displays the function $\mathcal H_d\widetilde\phi^{(d-1)}$ for $d=2,3.$ Throughout this section the nominal level is fixed as $\alpha=0.05$, and all level and power statements are in percent. Except for Table (ref), none of the simulations in this section assumes knowledge of the design density $f_{\boldsymbol{\Theta}}$ or uses a parametric specification of it.
We follow the multiscale approach in Section (ref) to obtain a graphical representation of the monotonicity behavior for a bivariate density of random coefficients. To test the hypotheses (ref) we use (ref) with $\mathcal T_n=\mathcal {T}_\mathbf{t}\times\{h_0\}\times\mathcal T_\mathbf{v}.$ Here, $h=h_0=0.5$ is fixed and the set of test locations $\mathcal {T}_\mathbf{t}$ is defined as the set of vertices on an equidistant grid in the square $[-1,2]^2$ with width one. Finally, the set of test directions is $$\mathcal {T}_\mathbf{v}=\big\{\mathbf{v}_1=-\mathbf{v}_3=\sqrt{2}^{-1}(1,1)^\top,\mathbf{v}_2=-\mathbf{v}_4=\sqrt{2}^{-1}(-1,1)^\top\big\}.$$ The data are simulated with $f_{\boldsymbol{\beta}}$ the density of the normal mixture $\tfrac 13 {\cal N}((-0.4,-0.57)^\top,0.2I)+ \tfrac 13 {\cal N}((1.5,-0.52)^\top,0.2I) + \tfrac 13 {\cal N}((0.45,1.6)^\top,0.15I).$ The design is chosen such that ${\boldsymbol{\Theta}}_i$ is uniformly distributed on the sphere $\mathbb S^1$. Figure (ref) in Section (ref) displays the monotonicity behavior of the density $f_{\boldsymbol{\beta}}$ based on sample size $n=20000$. Each arrow at a location $\mathbf{t}$ in direction $\mathbf{v}$ displays a rejection of a hypothesis (ref). The map indicates the existence of modes around the points $(-0.5,-0.5)^\top$, $(1.5,-0.5)^\top$, and $(0.5,1.5)^\top$ and thus detects the true modes fairly well.
Given the random coefficients model with $d=3$, we study the power of the test for the existence of a mode at a given location $\mathbf{b}_0$ considering only few local tests. The postulated mode is given by the point $\mathbf{b}_0 =(0,0,0)^\top$ and we take $\mathcal T_n =\{ (\mathbf{t}, h, \mathbf{v}): h=1,\mathbf{t}=\mathbf{v} \in \{\pm \mathbf{e}_i, i=1,2,3\}\}$ with $\mathbf{e}_i$ the $i$-th standard unit vector in $\mathbb{R}^3.$ We conclude that $f_{\boldsymbol{\beta}}$ has a local maximum at the point $\mathbf{b}_0$, whenever all hypotheses $H_{0,-}^{\mathbf{t}_j,h_0,\mathbf{v}_j},$ $j=1,\hdots,6,$ are rejected, i.e.
where $\kappa^{\mathbf{t}_j,h_0,\mathbf{v}_j}_n(\alpha) $ is defined by (ref). Recall that the quantiles $\kappa^{\mathbf{t}_j,h_0,\mathbf{v}_j}_n(\alpha) $ are constructed in such a way that the probability of at least one false rejection within the six tests (ref) is bounded by $\alpha$. However, the mode test detects the presence of a mode whenever all six tests (ref) are rejected at the same time. The multiscale method is therefore rather conservative for the specific task of mode detection. In this simulation, we also study a calibrated version of the test where the quantiles are chosen such that the test keeps its nominal level $\alpha=0.05$ and detects the presence of a non existing mode in about 5 percent of the simulation runs. For the calibration of the test we work under the null hypothesis assuming that $f_{{\boldsymbol{\beta}}}$ is uniform. Therefore, knowledge about the true unknown density $f_{\boldsymbol{\beta}}$ is not required.
Numerical simulations for random coefficients model without intercept: At first, we consider model (ref) with uniform design $\mathbf{X}_i \sim \operatorname{Unif}[-5,5]^3.$ To study the level of the test, we used $f_{\boldsymbol{\beta}}({\boldsymbol{\beta}}) \propto \mathbf{1}({\boldsymbol{\beta}} \in [-5,5]^3).$ For the power we took $f_{\boldsymbol{\beta}}$ as the density of a trivariate standard normal distribution. All results are based on the local tests (ref) and 1000 repetitions. Level simulations with the theoretical quantiles confirm that the multiscale test keeps its nominal level as the percentage of false rejections of at least one of the six hypotheses in (ref) for sample size $n\in\{250,500,1000\}$ is nearly $5$ percent. The results of mode test are reported in Table (ref).
Next, we investigate an asymmetric distribution of the directions ${\boldsymbol{\Theta}}_i$ by sampling $\mathbf{X}_i \sim \mathcal{N}((3,0,0)^\top,2I)$, with $I$ the $3\times 3$ identity matrix. We only consider the calibrated mode test. Results are reported in columns five and six in Table (ref). Compared to the case of uniform design, we observe a decrease in the power of the mode test (ref). The explanation is that the uniform design on the $\mathbf{X}_i$ induces a more uniform distribution of ${\boldsymbol{\Theta}}_i$ on the sphere which makes it simpler to recover information about the joint density as discussed in Section (ref).
Numerical simulations for random coefficients model with intercept: We study model (ref) with $d=3.$ In a first simulation, we sample the random vectors $(X_{i,2},X_{i,3})^\top$ from a standard bivariate Cauchy distribution, such that the density $f_{\boldsymbol{\Theta}}$ is constant. Except for the different design, we consider otherwise the same test settings as above. The simulated level and power of the calibrated version of the test (ref) are reported in Table (ref). To investigate the influence of the estimation of the design density $f_{\boldsymbol{\Theta}}$ on the power of the test we also perform simulations in which we assume that the density $f_{\boldsymbol{\Theta}}$ is known to be constant. These are shown in Table (ref), fourth and fifth column.
Compared to the power approximations for unknown $f_{\boldsymbol{\Theta}}$, we observe only a slight increase in power of the test for known $f_{\boldsymbol{\Theta}}$.
Finally, we consider two designs which do not satisfy Assumption (ref). Table (ref) reports the level and power for the same setting as above except that now $(X_{i,2}, X_{i,3})^\top$ is drawn from a standard normal distribution or $(X_{i,2}, X_{i,3})^\top \sim \operatorname{Unif}[-5,5]^2.$ We observe only a slight decrease in the power of the test for normally distributed design compared to the setting where Assumption (ref) holds. Even under uniform design, the test performs fairly well.
For multimodal densities which have a second mode close to the test location $\mathbf{b}_0$ testing different bandwidths simultaneously can be advantageous to separate the modes. This is illustrated by the following example, where we consider the random coefficients model without intercept with $d=2.$ The data are simulated with $f_{\boldsymbol{\beta}}$ being the density of the normal mixture $$\frac 12 \mathcal {N}\Big((0,0)^\top,\Big( \;
\;\Big)\Big) + \frac 12 \mathcal {N}((2,0)^\top,0.1\cdot I).$$ We consider a design such that ${\boldsymbol{\Theta}}$ is uniformly distributed on the circle $\mathbb{S}^1$. The level $\alpha$ is fixed to five percent. We conducted simultaneously twelve tests with three different scales $h\in \{0.5, 1, 2.5\}$ for the hypotheses \eqref{t5neu2} with $\mathbf{b}_0=(0,0)^\top$. The tests are given by $\{(\mathbf{t}_i, h_j, \mathbf{v}_i): i=1,\ldots,4; h_j \in \{0.5, 1, 2.5\}, \mathbf{v}_i\in \{\pm \mathbf{e}_1,\pm \mathbf{e}_2\}, \mathbf{t}_i=h_j\mathbf{v}_i\}.$ We analyze the outcome of the twelve tests in two ways. Firstly, we investigate the performance of each of the three tests for modality for the bandwidths $h=2.5$, $h=1$ and $h=0.5$ separately. Secondly, we consider the performance of a combined test for modality using two scales $h=1$ and $h=0.5$. The power approximations for sample sizes $n\in\{2000,5000,15000\}$ based on 1000 repetitions are displayed in Table (ref).
Figure (ref) illustrates the results of the twelve tests for the hypotheses (ref) conducted simultaneously. Each arrow at a location $\mathbf{t}$ in direction $\mathbf{v}$ displays a rejection of a hypothesis (ref) and the length of the arrows corresponds to the respective bandwidths.
The results in Table (ref) show that the mode test for $h=2.5$ detects in all cases a mode in a neighborhood of $(0,0)^\top$. However, the bandwidth $h=2.5$ is too large to distinguish between the underlying modes of $f_{\boldsymbol{\beta}}$ at $(0,0)^\top$ and $(2,0)^\top$. The test for the bandwidth $h_0=1$ fails to detect the mode as $h=1$ is still too large to separate the two modes of $f_{\boldsymbol{\beta}}.$ On the contrary, $h_0=0.5$ is too small to detect the decrease with small slope corresponding to the eigenvalue $0.4$ of the covariance matrix in the first mixture component of $f_{\boldsymbol{\beta}}.$ Note that this effect vanishes for increasing sample size. By conducting the tests for the bandwidths $h=1$ and $h=0.5$ simultaneously, we are able to detect the mode at $(0,0)^\top$ in most of the simulations.
It is common in the literature to use a parametric specification of the random coefficients. Usually, ${\boldsymbol{\beta}}$ is assumed to have a multivariate normal distribution. In this section, we compare our nonparametric approach with mode estimation under a parametric specification of the random coefficient distribution. Our findings are similar to what is usually observed when parametric and nonparametric methods are compared. If the parametric specification of the model is consistent with the data generating process, the parametric approach outperforms the nonparametric one in terms of estimation errors, power of tests, and computation time. However, when the parametric specification and the data generating process differ, a model misspecification bias is present in the parametric estimation. This will not vanish, even if the sample size is large. The nonparametric model does not suffer from this problem and can in these cases perform better than the parametric model. We illustrate this effect in the following simulations.
We consider the parametric random coefficients model $Y=\tilde{\boldsymbol{\beta}}^\top \mathbf{X}$ where the density of $\tilde{\boldsymbol{\beta}}$ belongs to some parametric family $f_{\tilde{\boldsymbol{\beta}}} (\mathbf{b};\eta)$ with parameter $\eta$. Note that we can rewrite the equation as a mixed model $Y = {\boldsymbol{\gamma}}^\top \mathbf{X} + \mathbf{u}^\top \mathbf{X}$. Here ${\boldsymbol{\gamma}} = \mathbb{E}(\tilde{\boldsymbol{\beta}})$ are fixed effects and $\mathbf{u} = \tilde{\boldsymbol{\beta}} - \mathbb{E}(\tilde{\boldsymbol{\beta}})$ are random effects with expectation $\boldsymbol 0$. If $f_{\tilde{\boldsymbol{\beta}}} (\mathbf{b};\eta)$ is unimodal with a mode at $\mathbb{E}(\tilde{\boldsymbol{\beta}})$, we just need to derive confidence statements for ${\boldsymbol{\gamma}}$. This is in particular true if $f_{\tilde{\boldsymbol{\beta}}} (\mathbf{b};\eta)$ is a normal density, which is the most common choice. A simple and efficient way to estimate ${\boldsymbol{\gamma}}$ is heteroscedasticity robust linear regression. Obviously, the computational complexity of this method is much smaller than the complexity of our algorithm.
Our test example is the random coefficients model with intercept (ref) with $d=3$ and $(X_{i,2},X_{i,3})^\top$ sampled from a standard bivariate Cauchy distribution. The distribution of the random coefficients ${\boldsymbol{\beta}}_i$ in the data generating process is given by (i) a standard Gaussian, (ii) the normal mixture $0.5\mathcal N(\boldsymbol 0,0.1I)+0.5\mathcal {N}((2,0,0)^\top,0.1I)$ and (iii) an exponential-2 distribution for $\beta_{i,1}$ and $(\beta_{i,2},\beta_{i,3})^\top\sim\mathcal N(\boldsymbol 0,0.1I)$ independent of the first component. Table (ref) reports estimates for ${\boldsymbol{\gamma}}$ obtained by transforming $Y$ and $\mathbf{X}$ to $S$ and ${\boldsymbol{\Theta}}$ and running a heteroscedasticity robust regression.
In column (i) of Table (ref) the parametric assumption holds and the procedure detects the true mode $\boldsymbol 0$ of the density with high precision compared to the bandwidth choice $h=1$ in the simulations presented in Table (ref) in Section (ref). In contrast, for the bimodal density in column (ii) the coefficient vector $(1,0,0)^\top$ does not describe a representative member of the population because the misspecification bias is too large. For the skewed distribution of column (iii) the estimator also fails to detect the mode of the density for the same reason.
By applying our testing procedure in the setting of columns (ii) and (iii) we can show that the OLS results do not represent modes of the density. To this end, we set in (ii) \[\mathbf{t}=(0.5,0,0)^\top, \quad h=0.5, \quad \mathbf{v} =(1,0,0)^\top\] and in (iii) \[\mathbf{t}=(1,0,0)^\top, \quad h=1, \quad \mathbf{v} =(1,0,0)^\top.\] Both tests reject $H_{0,-}^{\mathbf{t},h,\mathbf{v}}$ (in (ii) in $100\%$ and in (iii) in $95.5\%$ percent of 1000 repetitions). Therefore, our procedure shows that neither $(1,0,0)^\top$ in \textit{(ii)} nor $(2,0,0)^\top$ in \textit{(iii)} are modes of the underlying density. Of course, the parametric model would be able to detect the modes in cases \textit{(ii)} and \textit{(iii)} with a different parametric specification. However, this would require considerable a priori knowledge about the data. If we did not interpret the results as estimators for the mode but as estimators for $\mathbb{E}({\boldsymbol{\beta}})$, the parametric method would perform well.
Heterogeneity of consumers is a major challenge in modeling and estimating consumer demand. In several different demand models random coefficients were proposed to account for the heterogeneity in the population of consumers.
In this section we are interested in the almost ideal demand system (AIDS) which was initially proposed by Deaton:80 with fixed coefficients. This model does not explain demand for a product itself but explains the budget share spent on a product by a linear equation. The explanatory variables are log prices and the log of total expenditure divided by a price index. A detailed discussion of the model is contained in Lewbel:97.
Fixed coefficients in this model mean that all consumers are assumed to react in the same way when the price of a product changes. It is well known that some consumers are very price sensitive and change their behavior significantly with small variations in prices while other consumers are less price sensitive. This type of heterogeneity can be modeled by a random coefficient on log prices which is assumed to vary across the population of consumers. A similar argument suggests a random coefficient on log total expenditure. Recently, applications of the AIDS using a nonparametric random coefficient specification instead of fixed coefficients were presented in HHM:15 and Breunig16.
We apply our multiscale test to detect modes in the random coefficients for budget shares for food at home (BSF)
Food expenditure is a large fraction of total expenditure and is roughly about 20%.
We analyze the data of the British Family Expenditure Survey which ran from 1961 to 2001. It reported yearly cross sections for household income, expenditure and other characteristics of roughly 7000 households. We use data of the years 1997--2001 only which gives a sample size of about 33000. Budget shares of food are generated by dividing the expenditure for all food by total expenditure. Food prices are reported as relative prices in comparison to a general prize index. The variable $TotExp$ is normalized to January 2000 real prizes.
Assumption (ref) and the numerical simulations in Section (ref) suggest that our test has more power when the normalized regressors are approximately uniform on the sphere. We can achieve this by symmetrizing the design in model (ref) as follows:
The relation of the modified model to the random coefficients in (ref) is $\beta_{i,1} = \tilde \beta_{i,1} - 5 \tilde \beta_{i,2} - 0.3 \tilde \beta_{i,3}$, $\beta_{i,2} = \tilde \beta_{i,2}$, $\beta_{i,3} = 25 \tilde \beta_{i,3}$. Observations of the new variable $\ln(TotExp_i)-5$ lie between $-5$ and $3.7$. The observations of $25\ln(FoodPrice_i) - 0.3$ range from $-1$ to $1.3$.
For a first evaluation of the data we assumed fixed coefficients in model (ref) and estimated the model with ordinary least squares (OLS).
In order to find modes of the density we conducted simultaneously tests on the $5\%$ level of the form (ref) on the two scales $h_1=0.75$ and $h_2=0.5$. Recall from Section (ref) that our testing procedure also performs well when the $\mathbf{X}_i$ are not Cauchy distributed. Therefore, to obtain a testing procedure which is more flexible with respect to the design, we use the nonparametric density estimator $\widehat{f}_{\boldsymbol{\Theta}}$ instead of a parametric estimation procedure. We were testing for modes on the equidistant grid covering $[-1,1]^3$ with grid width $1$. Hence, the grid had 27 nodes. For every grid point $\mathbf{b}\in\mathbb{R}^3$ tests of the hypotheses (ref) were conducted for the directions and locations
where $\mathbf{e}_1,\mathbf{e}_2,\mathbf{e}_3 \in \mathbb{R}^3$ denote the standard unit vectors of $\mathbb{R}^3$. We detected a single mode in the neighborhood of the grid point $(0,0,0)^\top$ for the tests with bandwidth $h_1=0.75$. The test for the bandwidth $h_2=0.5$ did not detect a mode.
In the following we use nonparametric density estimation to motivate hypotheses for the testing procedure. It is important that this estimate and the test are independent, otherwise the testing procedure would be biased and could not guarantee a bound on the error rate. We meet the requirement by splitting the sample in two independent equally sized sub-samples. The first sub-sample is used for nonparametric estimation of the random coefficient density in model (ref) with the estimator in hoderlein2008. Figure (ref) gives contour plots for the joint densities of $f_{\tilde\beta_1,\tilde\beta_2}$, $f_{\tilde\beta_1,\tilde\beta_3}$, and $f_{\tilde\beta_2,\tilde\beta_3}$ based on about 16500 observations. We chose the smoothing parameters $h$ and $g$ in the estimator in hoderlein2008 equal to 0.05 and 0.1, respectively. Note that these bandwidth choices do not affect the level of the test performed below as the test does not depend on this estimator. The nonparametric estimate suggests that the random coefficient density of $f_{\tilde\beta_1,\tilde\beta_2,\tilde\beta_3}$ has one (well-pronounced) mode close to
This is consistent with the results of the test above which found a mode close to $(0,0,0)$. Since the marginal densities of $\beta_1,\beta_2,\beta_3$ are nearly symmetric it is also consistent that the mode is close to the OLS estimates given in Table (ref). With a significantly skewed or with a multimodal random coefficient density location of modes would differ from OLS.
With the second sub-sample we studied whether a mode can be found for the location given by (ref) on a smaller scale than in the test above. This would then indicate that the location of the mode is not $(0,0,0)^\top$. Our primary interest is in the coefficients on $\tilde\beta_2$ and $\tilde\beta_3$ on total expenditure and food prices. In order to see if the mode is indeed in a location where $\tilde\beta_2<0$ and $\tilde\beta_3>0$, we conduct two mode tests (ref) simultaneously for the bandwidths $h_1=0.07$ and $h_2=0.02$, that is, we test twelve hypotheses (ref) for the following locations and directions $(\mathbf{t}_i, \mathbf{v}_i).$ On scale $h_1=0.07$, we tested
and on scale $h_2=0.02,$
The test rejected all local hypotheses on scale $h_1=0.07$ but not on scale $h_2=0.02.$ This gives evidence that the mode is in a location where $\tilde\beta_2$ is negative but we cannot decide whether $\tilde\beta_3$ is positive at the mode.
Let us return to the initial model (ref). The results of our test give evidence that a mode exists close to \[ (\beta_1,\beta_2,\beta_3) = (0.5,-0.07,0.5) \] with strong evidence that $\beta_2$ is indeed negative. This vector of coefficients describes a representative member of the majority of consumers. It suggests that in the majority group food budget shares decrease with increasing log total expenditure. The nonparametric estimate in Figure (ref) shows that there is considerable variance among consumers around this representative member.
J. Schmidt-Hieber was partially funded by a TOP II grant from the Dutch science foundation. F. Dunker acknowledges support by the Ministry of Education and Cultural Affairs of Lower Saxony in the project Reducing Poverty Risk. K. Proksch acknowledges financial support by the German Research Foundation DFG through subproject A07 of CRC 755. K. Eckle has been supported by the Collaborative Research Center “Statistical modeling of nonlinear dynamic processes” (SFB 823, Project C4) of the German Research Foundation (DFG). We are very grateful to two referees and an associate editor for their constructive comments, which led to substantial improvement of an earlier version of this manuscript.
\setcounter{page}{1}