Kernel Consistent Univariate Density Equality Test with Mixed Data Types

npunitest implements the consistent metric entropy test of Maasoumi and Racine (2002) for two arbitrary, stationary univariate nonparametric densities on common support.

Usage

npunitest(data.x = NULL,
          data.y = NULL,
          method = c("integration","summation"),
          bootstrap = TRUE,
          B = 399,
          bw.x = NULL,
          bw.y = NULL,
          random.seed = 42,
          ...)

Arguments

Data, Bandwidth Inputs And Formula Interface

These arguments identify the two samples, statistic variant, and any supplied bandwidths.

bw.x, bw.y: numeric (scalar) bandwidths. Defaults to plug-in (see details below).
data.x, data.y: common support univariate vectors containing the variables.
method: a character string used to specify whether to compute the integral version or the summation version of the statistic. Can be set as integration or summation. Defaults to integration. See ‘Details’ below for important information regarding the use of summation when data.x and data.y lack common support and/or are sparse.

Bootstrap Controls

These arguments control bootstrap execution and reproducibility settings.

B: an integer value specifying the number of bootstrap replications to use. Defaults to 399.
bootstrap: a logical value which specifies whether to conduct the bootstrap test or not. If set to FALSE, only the statistic will be computed. Defaults to TRUE.
random.seed: an integer used to seed R's random number generator. This is to ensure replicability. Defaults to 42.

Additional Arguments

Further arguments are passed to the bandwidth-selection routines used by the test.

...: additional arguments supplied to specify the bandwidth type, kernel types, and so on. This is used since we specify bw as a numeric scalar and not a bandwidth object, and is of interest if you do not desire the default behaviours. To change the defaults, you may specify any of bwscaling, bwtype, ckertype, ckerorder, ukertype, okertype.

Bootstrap P-Values

The bootstrap upper-tail p-value is the proportion of bootstrap statistics greater than or equal to the observed statistic. Ties count as at least as extreme, giving the conservative tie convention for discrete statistics. This is the ordinary bootstrap proportion, without an add-one correction; a reported zero is possible with a finite number of replications.

Value

npunitest returns an object of type unitest with the following components

Srho: the statistic Srho
Srho.bootstrap: contains the bootstrap replications of Srho
P: the P-value of the statistic
boot.num: number of bootstrap replications
bw.x, bw.y: scalar bandwidths for data.x, data.y

summary supports object of type unitest.

Book And Method Pointers

npunitest uses an entropy-distance measure to compare univariate density or probability functions. The statistic is based on the normalized Hellinger-type distance described in the entropy-test literature referenced below.

For book-length background, see Racine (2019), Chapter 2 Continuous Density and Cumulative Distribution Functions, especially the entropy and density-test material. For related density and testing background, see Li and Racine (2007), Chapters 1, 12, and 13.

References

Granger, C.W. and E. Maasoumi and J.S. Racine (2004), “A dependence metric for possibly nonlinear processes”, Journal of Time Series Analysis, 25, 649-669.

Maasoumi, E. and J.S. Racine (2002), “Entropy and predictability of stock market returns,” Journal of Econometrics, 107, 2, pp 291-312.

Author

Tristen Hayfield tristen.hayfield@gmail.com, Jeffrey S. Racine racinej@mcmaster.ca

Details

Documentation guide: see np.kernels for kernels, np.options for global options, and plot for plotting options.

npunitest computes the nonparametric metric entropy (normalized Hellinger of Granger, Maasoumi and Racine (2004)) for testing equality of two univariate density/probability functions, \(D[f(x), f(y)]\). See Maasoumi and Racine (2002) for details. Default bandwidths are of the plug-in variety (bw.SJ for continuous variables and direct plug-in for discrete variables). The bootstrap is conducted via simple resampling with replacement from the pooled data.x and data.y (data.x only for summation).

Missing observations are omitted independently from the two samples before bandwidth selection and overlap checks. Each sample must retain at least one observation; automatic bandwidth selection can require more.

For categorical data, the declared factor levels define the finite domain of the probability comparison, including levels not observed in a sample. The two samples share the union of their declared levels, retained in all bootstrap samples. Plug-in bandwidths use that category count, and the discrete statistic sums over the entire common domain, not only observed categories. Existing kernel normalization is unchanged. Numeric ordered levels retain their numeric scale. Qualitative ordered samples must specify the same level order or an order-consistent subset; otherwise declare a common ordered level set explicitly before calling the test.

The summation version of this statistic can be numerically unstable when data.x and data.y lack common support or when the overlap is sparse (the summation version involves division of densities while the integration version involves differences, and the statistic in such cases can be reported as exactly 0.5 or 0). Warning messages are produced when this occurs (‘integration recommended’) and should be heeded.

For numeric data using the fixed second-order Gaussian density estimator, the integration statistic uses deterministic corrected composite trapezoidal quadrature. The integration domain extends eight bandwidths beyond the pooled observed support and uses four evaluation points per smallest bandwidth. During a default numeric iid bootstrap, resampling counts are evaluated in bounded chunks on the same quadrature grid for integration or on the original support for summation. This avoids repeated adaptive integration and repeated reconstruction of resampled data. The summation density ratios are evaluated by a registered native routine without forming a sample-by-sample kernel matrix.

The accelerated route is also used when ... contains only options explicitly equivalent to the defaults: bwtype="fixed", ckertype="gaussian", and ckerorder=2. Other density options retain the corresponding general estimator and numerical integration route. In particular, generalized and adaptive nearest-neighbour bandwidths use duplicate-preserving resamples and recompute their realized neighbour radii; bootstrap counts are not treated as fixed-bandwidth observation weights.

Usage Issues

See the example below for proper usage.

Examples

if (FALSE) { # \dontrun{
set.seed(1234)
n <- 1000

## Compute the statistic only for data drawn from same distribution

x <- rnorm(n)
y <- rnorm(n)

npunitest(x,y,bootstrap=FALSE)

if (interactive()) Sys.sleep(5)

## Conduct the test for this data

npunitest(x,y,B=99)

if (interactive()) Sys.sleep(5)

## Conduct the test for data drawn from different distributions having
## the same mean and variance

x <- rchisq(n,df=5)
y <- rnorm(n,mean=5,sd=sqrt(10))
mean(x)
mean(y)
sd(x)
sd(y)

npunitest(x,y,B=99)

if (interactive()) Sys.sleep(5)

## Two sample t-test for equality of means
t.test(x,y)
## F test for equality of variances and asymptotic
## critical values
Fstat <- var(x)/var(y)
qf(c(0.025,0.975),df1=n-1,df2=n-1)

## Plot the nonparametric density estimates on the same axes

fx <- density(x)
fy <- density(y)
xlim <- c(min(fx$x,fy$x),max(fx$x,fy$x))
ylim <- c(min(fx$y,fy$y),max(fx$y,fy$y))
if (interactive()) plot(fx,xlim=xlim,ylim=ylim,xlab="Data",main="f(x), f(y)")
lines(fy$x,fy$y,col="red")

if (interactive()) Sys.sleep(5)

## Test for equality of log(wage) distributions

data(wage1)
with(wage1, {
lwage.male <- lwage[female=="Male"]
lwage.female <- lwage[female=="Female"]

npunitest(lwage.male,lwage.female,B=99)

if (interactive()) Sys.sleep(5)

## Plot the nonparametric density estimates on the same axes

f.m <- density(lwage.male)
f.f <- density(lwage.female)
xlim <- c(min(f.m$x,f.f$x),max(f.m$x,f.f$x))
ylim <- c(min(f.m$y,f.f$y),max(f.m$y,f.f$y))
plot(f.m,xlim=xlim,ylim=ylim,
     xlab="log(wage)",
     main="Male/Female log(wage) Distributions")
lines(f.f$x,f.f$y,col="red",lty=2)
rug(lwage.male)
legend(-1,1.2,c("Male","Female"),lty=c(1,2),col=c("black","red"))

})

if (interactive()) Sys.sleep(5)

## Conduct the test for data drawn from different discrete probability
## distributions

x <- factor(rbinom(n,2,.5))
y <- factor(rbinom(n,2,.1))

npunitest(x,y,B=99)
} # }