PySDTest: a Python/Stata Package for Stochastic Dominance Tests

PySDTest: a Python/Stata Package for Stochastic Dominance Tests

Introduction

Stochastic dominance (SD) is an ordering rule of distribution functions, originally suggested by hadar1969rules and hanoch1969efficiency. It is based on the expected utility paradigm and gives a uniform ordering of prospects that does not depend on a specific preference structure (or utility function) of a decision maker. Because of its generality and nonparametric feature, academics and practitioners in various areas, such as economics, finance, insurance, medicine and statistics, have regarded SD ordering as a fundamental concept of decision-making. They have empirically examined SD relations in diverse areas. Examples include the portfolio choice problem post2003empirical,bali2009bonds, income inequality and welfare analysis amin2003does, anderson2003poverty, financial market efficiency cho2007there, auction bids de2003empirical, firm productivity and export decision delgado2002firm, obesity inequality pak2016measuring, and agricultural productivity langyintuo2005yield among many others. The number of applications has been growing along with recent advances in inference methods for testing SD hypotheses.

There exists a big literature for statistical tests of SD; see whang2019econometric for an extensive survey. mcfadden1989testing pioneered statistical tests of SD, proposing a Kolmogorov-Smirnov type test for the first and second order SD. Other early works are klecan1991robust, kaur1994testing, anderson1996nonparametric, davidson2000statistical. barrett2003consistent introduce resampling procedures for consistent testing of arbitrary SD order between two distributions under i.i.d. and mutually independent sampling schemes. \citet*[LMW]{linton2005consistent} propose a consistent subsampling method which is valid under general sampling schemes, including serially dependent observations and mutually dependent prospects. The LMW test also allows for prospects that may depend on some unknown finite-dimensional parameters, enabling comparison between residuals after controlling for observed confounders. \citet*[LSW]{linton2010improved} introduce the contact set approach to improve the power performance of SD tests and show large sample uniform validity of the test. Also, they allow the prospects to depend on unknown infinite dimensional parameters. donald2016improving provide an alternatvie method to enhance the performance by suggesting the selective recentering approach. However, to our knowledge, comprehensive statistical software for conducting SD tests is not available.

This paper introduces a Python package PySDTest and Stata command pysdtest for testing SD hypotheses, hopefully to help practitioners implement SD tests in a convenient and flexible way. The package implements various SD tests developed in the recent literature. Supported testing procedures include, but are not limited to, the least favorable case approach of BD, subsampling approach of LMW, contact set approach of LSW, and selective recentering method of DH. In addition, we include the numerical delta method (NDM) based on the results of dumbgen1993nondifferentiable, fang2019inference, and hong2018numerical in a more general context to compute critical values, as the NDM may have computational advantages in some contexts. This is a novel result since there has been no practical software tailored to SD testing. We provide guidance on how to use our package in Python and Stata environments with practical examples.

PySDTest also allows for flexibly testing a complex and advanced null hypothesis based on SD relations. For example, practitioners can compare more than 2 prospects simultaneously. To be specific, given $K (\geq 3)$ prospects, they can perform a joint test whether there exists at least one prospect stochastically dominating some of the others, which is called stochastic non-maximality among multiple prospects klecan1991robust. If they reject the hypothesis, then they may confer that the set of prospects is stochastically maximal, in the sense that no prospect in the set is dominated by any of the others. This practice is possible because the package offers additional features such as diverse resampling methods and functional type of test statistics. We describe them in detail with an illustrative example.

As an empirical illustration, we apply PySDTest to the portfolio choice problem in the cryptocurrency and stock markets. We investigate the SD relation between the daily returns of Bitcoin price and S&P 500 index. Our testing result shows that the S&P 500 returns second order stochastically dominate the Bitcoin returns. This implies that a risk-averse economic agent with a monotone increasing utility function would prefer S&P 500 index to Bitcoin as a financial asset.

This paper is organized as follows. In Section 2, we introduce the concepts of SD and their testing methods. In Section 3, we provide a full description of the Python package PySDTest and Stata command pysdtest. In Section 4, we apply PySDTest to compare the daily return distributions of Bitcoin price and S&P 500 index. In Section 5, we conclude the paper.

Concepts and Tests of Stochastic Dominance

Concepts of SD

Let $X_{1}$ and $X_{2}$ be two (continuous) random variables with cumulative distribution functions (CDF) $F_{1}$ and $F_{2}$, respectively. Let $Q_{k} = \inf \{x: F_{k}(x) \geq \tau \}$ denote the quantile function of $X_{k}$ for $k=1,2$, respectively. In addition, let $\mathcal{U}_{1}$ be the class of all monotone increasing utility (or social welfare) functions. If the functions are differentiable, then $\mathcal{U}_1$ can be written as:

This definition can be understood by considering an investor's portfolio choice problem. Let the random variables be returns of some financial assets. The condition $F_{1}(x) \leq F_{2}(x)$ for all $x \in \mathbb{R}$ can be denoted as:

which means, for all $x \in \mathbb{R}$, the probability of attaining returns larger than $x$ is larger under $F_{1}$ than $F_{2}$. An investor with a monotone increasing utility function would therefore prefer $F_{1}$ to $F_{2}$. If the two CDFs intersect, however, FSD does not hold. This means that there exist some $u$ and $v$ both in $\mathcal{U}_{1}$ such that $E[u(X_{1})] > E[u(X_{2})]$, but $E[v(X_{1})] < E[v(X_{2})]$.

The second-order stochastic dominance can be defined in a similar manner. Let $\mathcal{U}_{2}$ be the class of all monotone increasing and concave (utility or social welfare) functions. If functions are differentiable, then we can write the class $\mathcal{U}_2$ as:

Consider the portfolio choice problem of an investor again. The concavity of utility function reflects risk-aversion of the investor. If $X_{1} \succeq_{2} X_{2}$, which represents $E[u(X_{1})] \geq E[u(X_{2})]$ for all $u \in \mathcal{U}_{2}$, the risk-averse investor would prefer $X_{1}$ to $X_{2}$.

Higher-order SD relations also can be defined. Define integrated CDF and integrated quantile function as davidson2000statistical:

The higher SD order corresponds to the smaller class of utility (or social welfare) functions. That is, for $s \geq 1$, the class of utility functions is defined as:

For example, the third-order stochastic dominance ($s=3$) reflects preference toward positive skewness of a distribution (whitmore1970third, whitmore1978stochastic, and levy2016stochastic). That is, investors would prefer positively skewed distributions to negatively skewed distributions.

Tests of SD

Hypothesis of interest

For a given stochastic order $s$, consider the following hypotheses of interest:

The null hypothesis means that random variable $X_1$ $s$-th order stochastically dominates $X_2$ and the alternative hypothesis is the negation of the null hypothesis.\footnote{This type of hypotheses has mainly been considered for statistical inference of the SD relations in the literature. There are two other types of hypotheses that are also considered. The first one is the null hypothesis of non-dominance and the alternative hypothesis of dominance. The other type is the null hypothesis of equivalent distributions and the alternative hypothesis of (strict) stochastic dominance. In this paper, we do not cover these hypotheses. For further discussions, see Section 2.1 of whang2019econometric.}

BD Test: Least Favorable Case Approach

Let $\{X_{1,i} : i = 1,2,...,N_1\}$ and $\{X_{2,i} : i = 1,2,...,N_2\}$ be random samples from $F_1$ and $F_2$, respectively. We assume $N_1 = N_2 \equiv N$ for simplification.\footnote{PySDTest allows for the case with different sample sizes $N_1 \neq N_2 $.} Let $\mathcal{X}$ denote the common compact support of $X_1$ and $X_2$.

We first introduce and define frequently used notions in the SD testing literature. The integrated CDF, $F_k^{(s)}(\cdot)$, can be written through integration by parts as davidson2000statistical:

where $\left[\cdot \right]_{+} = \max\{\cdot,0\}$. The natural empirical counterpart, empirical integrated CDF, of (ref) is given by

Define the difference between integrated CDFs and empirical integrated CDFs as:

BD propose the supremum (or Kolmogorov-Smirnov) type test statistic

to test the hypothesis (ref). Note that we can write (ref) as:

where $\mathcal{J}_s(x,F_k)$ is a linear and continuous functional with respect to $F_k$. Under i.i.d. random sampling and compact support $\mathcal{X}$, the following weak convergence result holds:

as $N\rightarrow \infty$, where $\mathcal{B}_1$ and $\mathcal{B}_2$ are independent standard Brownian bridges and $\Rightarrow$ denotes the weak convergence. This result comes from applying Donsker theorem to the case $s=1$ and continuous mapping theorem for the case $s \geq 2$. The limit process, $\mathcal{J}_s(\cdot, \mathcal{B}_k \circ F_k)$, is a Gaussian process with the covariance function given by

Define

Equation (ref) implies

Then, under the null hypothesis, we have the following result:

Let $c_{1-\alpha}^{BD}$ is the $(1-\alpha)$ quantile of $\sup_{x \in \mathcal{X}}\nu_s(x)$ for $\alpha \in [0,1]$. BD show the following results:

The equality of (ref) holds when the least favorable case (LFC), $F_1 = F_2$, holds. Note that the limit distribution $\sup_{x \in \mathcal{X}}\nu_s(x)$ depends on the true CDF, $F_k$, hence it cannot be tabulated.

BD propose three procedures to compute the critical values. All these methods impose the LFC restriction under the null hypothesis and try to mimic the distribution $\sup_{x\in \mathcal{X}} \nu_s(x)$. The first one is the multiplier method. It first generates two independent random samples $\{U_{k,i}:i=1,2,...,N\}$ from $N(0,1)$ for $k=1,2$. Then, simulate the process that mimics the $F_k-$Brownian bridge $\mathcal{B}_{k}\circ F_{k}(x)$:

The multiplier method generates the following simulated process:

for $k=1,2$. Then, the critical value for the significance level $\alpha$ is defined as the $(1-\alpha)$ quantile of the following simulated distribution:

Note that $BD_{N,MP}^*$ approximates the distribution of $\sup_{x \in \mathcal{X}} \nu_{s}(x)$ conditional on the original sample.

The second method is the pooled-sample bootstrap method. This method draws bootstrap samples $S_{k}^* = \{X_{k,i}^* : i=1,2,...,N\}$ for $k=1,2$ from the pooled original sample $\{X_{k,i}:k=1,2 \text{ and } i = 1,2,...,N\}$ and calculates the empirical CDFs by using the bootstrap samples as:

Then, we compute the following distribution of the bootstrap test statistic:

The critical value for the significance level $\alpha$ is defined to be the $(1-\alpha)$ quantile of the $BD_{N,pool}^*$ distribution. Note that the LFC restriction is imposed by pooling the original samples.

The third procedure is the recentered bootstrap method. This method draws bootstrap samples $S_{1}^*$ and $S_{2}^*$ from $\{X_{1,i}:i=1,2,...,N\}$ and $\{X_{2,i}:i=1,2,...,N\}$, respectively. Then, the empirical CDFs are calculated using each bootstrap sample as in (ref). Define

where the LFC restriction is imposed by recentering the bootstrap statistic $\mathcal{J}_s(x, \Bar{F}^*_k)$ by subtracting its mean $\mathcal{J}_s(x, \Bar{F}_k)$ for $k=1,2$. The critical values are defined to be the $(1-\alpha)$ quantile of the bootstrap distribution $BD_{N,recenter}^*$.

LMW Test: Subsampling Apporach for General Sampling Scheme

For a given set of $K$ prospects, LMW considered the following hypotheses:

where $D_{k,l}^{(s)}(x)$ is defined as in (ref) and $\mathcal{X}$ is the union of the supports of $X_k$ for $k=1,2,...,K$. The null hypothesis $H_0$ is the hypothesis of stochastic non-maximality (klecan1991robust). It says that there exists at least one prospect that $s$-th order stochastically dominates some of the others in the set of the $K$-prospects. The alternative hypothesis is that of stochastic maximality, that is there does not exist any prospect in the set which $s$-th order stochastic dominates any of the others.

LMW propose a supremum type test for testing (ref) with the following features: (i) SD of arbitrary order, (ii) Subsampling critical values allowing for general sampling scheme such as serial dependence of observations and arbitrary dependence amongst prospects at a specific time, (iii) Prospects may depend on finite dimensional parameters, enabling conditional ranking. For instance, the test can be applied to residuals from regressions to control for systematic factors.

Let $X_{k,i}(\theta)$ be the prospect that may depend on a finite dimensional parameter $\theta$ for $k=1,2...,K$. If we consider the residual from the linear regression, we may take $X_{k,i}(\theta) = Y_{k,i} - Z_{k,i}^\intercal \theta$. Define

Likewise, we define $\Bar{F}_k^{(s)}(x,\theta)$ as the sample analogue of $F_{k}^{(s)}(x,\theta)$ for $s \geq 1$.

Let $F_{k}^{(s)}(x) := F_{k}^{(s)}(x,\theta_{k,0})$ and $\Bar{F}_{k}^{(s)}(x) := \Bar{F}_{k}^{(s)}(x,\theta_{k,0})$, where $\theta_{k,0}$ denotes the true value, and define

We can represent the hypotheses (ref) and (ref) as:

LMW use the sample analogue of $d_s$ as the basis of their test statistic:

where $\hat{\theta}_k$ is a consistent estimator of $\theta_{k,0}$ for $k=1,2,...,K$. The pointwise asymptotic mull distribution of $LMW_{N}$ is a functional of a mean zero Gaussian process when $d_s = 0$ (the “boundary" case), and degenerates to $-\infty$ when $d_s <0$ (the “interior" case) under some regularity conditions. For illustration, let $K=2$ and there be no parameter to be estimated.\footnote{The function pysdtest.test_sd() in the PySDTest package with the `subsampling' argument implements the routine of the LMW test for a two-sample and no parameter case. However, users can implement the test with multiple prospects and parameters by using the additional features provided by PySDTest.} Then, the test statistic $LMW_N$ can be denoted as:

To describe the asymptotic null distribution, define the contact set to be:

Then, under the null hypothesis $H_0$,

where $\nu_{k,l}^{(s)}(\cdot)$ is a mean zero Gaussian process for $k,l = 1,2$ and $k \neq l$. \footnote{This result is based the pointwise asymptotics and hence does not fully describe the behavior of the test statistic in finite samples. See Section 2.2.4 for a result based on the uniform asymptotics.} Under the alternative hypothesis $H_1$ ($d_s > 0)$, $LMW_N \rightarrow_{p} \infty$ where $\rightarrow_p$ denotes convergence in probability. Note that the asymptotic null distribution depends on the true distribution hence is not pivotal.

LMW suggest a subsampling procedure to compute the critical values.\footnote{The subsampling method has been proposed by politis1994stationary, see politis1999subsampling for a comprehensive review.} This method not only circumvents the non-pivotal issue but also yields consistent critical values under a general sampling scheme. Let $b$ be the size of subsamples, satisfyng the condition $b \rightarrow \infty$ and $b/N \rightarrow 0$ as $N \rightarrow \infty$. Let $\alpha \in [0,1]$ be the nominal significance level. Then the subsampling procedure can be conducted as the following steps:

Under Assumptions 1-4 of LMW, the subsampling test is asymptotically valid in the sense

In practice, we need to pick the subsample size $b$ which might affect the test result. Given that there is no compelling theory on the optimal choice of $b$ in the literature, LMW propose three data dependent methods for choosing $b$: minimum volatility method politis1999subsampling, mean critical value and median critical value. Let $B_N$ be a set of subsampling size candidates. Then, the minimum volatility method suggests practitioners select one in $B_N$ that minimizes local variations of critical values. The mean and median critical value methods recommend to pick the mean and median of critical values, respectively. However, LMW also suggest computing and drawing a plot of p-values, which checks whether p-values are insensitive to the choice of $b$. This is because, if so, inferences would be robust regardless of the specific value of $b$.

Contact Set Approach

The classical approach imposing the LFC restriction to obtain critical values can be too conservative in practice. This is because it tries to estimate the upper bound of the asymptotic null distributions. To see this, one can establish that under the null hypothesis,

where $\mathcal{C}_0$ is the contact set defined in (ref). Clearly, $\sup_{x \in \mathcal{C}_0}\nu_s(x)$ is stochastically dominated by $\sup_{x \in \mathcal{X}}\nu_s(x)$ because $\mathcal{C}_0 \subseteq \mathcal{X}$. Therefore, the asymptotic size of LFC-based tests achieves the nominal level only when prospects of comparison have equal distributions. This leaves the room for power enhancement.

LSW instead propose the contact set approach which estimates the contact set, $\mathcal{C}_0$, and directly mimics the asymptotic null distribution (ref) for inference. This leads to enhanced power property, compared to the LFC-based procedure.

LSW consider an integral-type test for arbitrary SD orders with the hypotheses of interests given by (ref). In addition to these features, LSW allow unbounded support and prospects depending on finite and infinite dimensional parameters and establish uniform size validity. Their tests can be applied to serially independent observations yet allow cross-sectional dependence between prospects.

Consider the one-sided Cram\'er-von Mises type functional:

where $q(\cdot)$ is a bounded weight function and $D_{1,2}^{(s)}(\cdot)$ is defined as in (ref).\footnote{LSW also allow practitioners to use a weight function for empirical analysis focusing on a specific part of distributions. For example, one may be interested in comparing income levels below the median of distributions. In this case, the functional would be defined to be

where $w(x)$ is a nonnegative integrable weight function.} The weight function here is to accommodate unbounded supports, which may entail the integrability issue of $F_k^{(s)}$ for $s \geq 2$. If $\mathcal{X}$ is bounded or $s=1$, $q(x) = 1$ is allowed. Otherwise, LSW suggest

for $z_1 < z_2$ and $a,\delta > 0$ where $a\vee b$ denotes $\max\{a,b\}$. Then the hypotheses of interests (ref) can be represented as

Non/semi-parametric models are allowed with the following specification of prospects:

where $W = (Y,Z) \in \mathbb{R}^{d_W}$ is a random vector and $\varphi(\cdot;\theta,\tau)$ is a real-valued function with parameters $(\theta,\tau)\in \Theta \times \mathcal{T}$. $\Theta$ is a finite-dimensional Euclidean space and $\mathcal{T}$ is an infinite-dimensional space. For example, $X_k := X_k(\theta_0, \tau_0)$ may be specified as the residual from the partially linear regression, i.e. $X_k(\theta,\tau) = Y_k - Z_{k,1}^\intercal \theta - \tau(Z_{k,2})$. Also, we may consider the single-index model such as $X_{k} = Y_k - \tau(Z_{k,1}^{\intercal}\theta)$.

Let $X_{k,i}(\theta,\tau) = \varphi_{k}(W_i;\theta,\tau)$ and $\{W_i : i=1,2,...,N\}$ be a random sample. Define the integrated empirical CDF as

and their differences as

LSW propose the integral type test statistic

where $(\hat{\theta},\hat{\tau})$ is a consistent estimator of $(\theta,\tau)$.

For brevity, we describe the testing procedure under the case when there are no estimated parameters, $q(x)=1$, and $\mathcal{X}$ is bounded. Thus, we consider the following simplified test statistic

where $\Bar{D}_{1,2}^{(s)}(\cdot)$ and $\Bar{F}_{k}^{(s)}(\cdot)$ are defined to be the same as in (ref) and (ref), respectively. Denote the set of probabilities that satisfy the null hypothesis as $\mathcal{P}_0$. Also, let $\mathcal{P}_{00} := \{P \in \mathcal{P}_0 : \lambda(C_0) > 0 \}$ where $\lambda$ is the Lebesgue measure. Put it differently, $\mathcal{P}_{00}$ is the subset of $\mathcal{P}_0$ with positive Lebesgue measure on the contact set.

We shall describe the asymptotic behavior of $LSW_N$ for a given $P \in \mathcal{P}_0$. With some regularity conditions (Assumptions 1-3 of LSW),

where $\nu_{s}(\cdot)$ is a mean zero Gaussian process. The second equality comes from the definition of the contact set. The third equality holds because $\sqrt{N}D_{1,2}^{(s)}(x) \rightarrow -\infty$ on $\mathcal{X}\setminus \mathcal{C}_0$ under the null hypothesis. Lastly, the weak convergence follows from $\sqrt{N}\left[\Bar{D}_{1,2}^{(s)}(\cdot) - D_{1,2}^{(s)}(\cdot)\right] \Rightarrow \nu_{s}(\cdot)$ and the continuous mapping theorem. On the other hand, for $P \in \mathcal{P}_0 \setminus \mathcal{P}_{00}$, $\lambda(\mathcal{C}_0) = 0$ gives us $LSW_N \Rightarrow 0$. Therefore, we can summarize the asymptotic results under the null hypothesis as

LSW call $\mathcal{P}_{00}$ as the set of boundary points and $\mathcal{P}_{0} \setminus \mathcal{P}_{00}$ as the set of interior points. The notable thing is that there is a discontinuity in the asymptotic null distributions, depending on the true probability. Since there is no discontinuity in finite samples, the above discontinuity indicates that the uniform asymptotics, rather than pointwise asymptotics, would provide better approximation and complete asymptotic behaviors.

Let $\alpha$ be the nominal significance level and $E_{P}$ be the expectation under the probability $P$. Then, the definitions are

The notable thing is that LFC-based testing procedures might be asymptotically non-similar hence biased in large samples against a large set of alternatives. The bootstrap procedure proposed by LSW is asymptotically similar on a larger subset of $\mathcal{P}_0$ in which the LFC is a special case. Therefore, the procedure is asymptotically biased against a smaller class of alternatives than LFC-based tests and hence is more powerful.

Let $B$ denote the number of bootstrap repetitions and $\alpha$ be the nominal significance level. The procedure can be described as follows:

Note that $W_i = (X_{1,i},X_{2,i})^{\intercal}$ and step 3-4 should be omitted if there are no estimated parameters.

We shall introduce some asymptotic properties of the contact set approach. Suppose that Assumptions 1-4 of LSW hold. Under $H_0$,

for some $\Tilde{\mathcal{P}}_0 \subset \mathcal{P}_0$ where $c_{\alpha,N,\infty}^*$ is the critical value from the ideal bootstrap ($B = \infty$). In addition, for some $\Tilde{\mathcal{P}}_{00,N} \subset \Tilde{\mathcal{P}}_0$,

Equation (ref) shows that the test has asymptotically vaild size uniformly over a class of “regular” probabilities $\Tilde{\mathcal{P}}_{0} \subset \mathcal{P}_0$. In addition, equation (ref) indicates that test has asymptotically exact size on $\Tilde{\mathcal{P}}_{00,N} \subset \Tilde{\mathcal{P}}_0$. On the other hand, the test is consistent against a fixed alternative. That is, for fixed alternative $P \notin \mathcal{P}_0$ such that $\int_{\mathcal{X}} \max\{q(x)D_{1,2}(x),0 \}^{2}dx > 0$,

In addition, Theorem 4 of LSW shows that the contact set approach is locally more powerful than the LFC-based tests.

Selective Recentering Approach

DH propose the selective recentering approach to enhance the power performance of BD tests. The main idea is similar to the contact set approach: DH try to directly approximate the limit distribution of $BD_N$, $\sup_{x \in \mathcal{C}_0} \nu_s(x)$ defined in (ref), to attain the critical values. For this end, DH extend the recentering method introduced by hansen2005test. It is also closely related to the generalized moment selection approach of andrews2013inference in the moment inequality literature.

DH also consider the hypotheses of interest given by (ref) and they consider the Barrett-Donald test statisitc $BD_N$. They consider a sampling scheme similar to that of BD, but further discuss extension of their method to deal with weakly dependent data and mutual dependence between $X_1$ and $X_2$.

Let us briefly describe the selective recentering approach. Let

Then, we have that the asymptotic null distribution of

is the same as that of $BD_N$ because $\mu(x) = D_{1,2}^{(s)}(x)$ under the null. DH propose a resampling procedure that tries to mimic (ref) using the sample analogue of (ref). That is, first take a sequence $a_N$ of negative numbers such that $a_N \rightarrow -\infty$ and $a_N/\sqrt{N} \rightarrow 0$ and approximate (ref) via the recentering function:

Under regularity conditions, it is a uniformly consistent estimator in the sense that

Write the simulated processes defined in (ref), (ref), and (ref) as

To compute the critical values, DH suggest using the bootstrap distributions

Letting $\hat{c}_{\alpha,k}$ be the $(1-\alpha)$ quantile of $DH_{N,k}^{*}$, DH suggest the following critical value:

where $\eta$ is a small positive number such as $10^{-6}$. $\eta$ is introduced here to accommodate the situation where both the original test statistic and the bootstrap test statistic degenerate to zero, such as the interior case $D_{1,2}^{(2)} (x) < 0$ for all $x \in \mathcal{X}$.

In practice, two tuning parameters, $a_N$ and $\eta$, need to be selected. Based on the simulation experiments conducted by DH, the specific choice of $\eta$ is not very important as long as $\eta$ is a small positive number. However, $a_N$ should be chosen carefully because the choice of $a_N$ may affect the finite sample performance of the test. DH suggest using $a_N = -0.1\sqrt{\log\log(N)}$ which performed well in their simulations and empirical studies.

As regard to the asymptotic properties, DH show that under the null hypothesis,

For a fixed alternative, the test is shown to be consistent:

Numerical Delta Method

As an alternative method for attaining critical values, we can apply an inference method applicable to directionally differentiable functions. fang2019inference show the asymptotic validity of an inference method that employs a consistent estimate of the first order directional derivative when the target parameters are Hadamard directionally differentiable. Complementing their idea, hong2018numerical propose the method using numerical derivatives for approximating the limit distributions which eases the burden of analytical calculation.\footnote{Both fang2019inference and hong2018numerical are extensions of the pioneering work of dumbgen1993nondifferentiable that studies a validity of the bootstrap under a lack of differentiability.} Even though fang2019inference discuss the testing for SD as an application of their methods, they do not provide a formal result let alone the software for the SD tests.

We briefly discuss the NDM, which also has enhanced power properties compared to the classical LFC-based tests. The power enhancement is possible because the NDM also directly approximates the asymptotic null distribution, as in the contact set approach and selective recentering approach. Theoretical validity of applying the NDM to SD tests is discussed by lee2023testing. To be specific, lee2023testing introduce the test for time stochastic dominance (time SD) with application of the NDM and the static SD concept is regarded as a special case of the time SD concept with a single time period.

To describe testing for SD by the NDM, we first introduce the definition of Hadamard directionally differentiability of a map between Banach spaces shapiro1990concepts.

We say $\phi'_{\theta}(h)$ is (first-order) Hadamard directional derivative at $\theta$ in direction $h$. In addition, if $\phi'_{\theta}$ is linear, then $\phi$ is said to be Hadamard differentiable at $\theta \in \mathbb{D}_{\phi}$ tangentially to a set $\mathbb{D}_{0} \subset \mathbb{D}$.

Theorem 2.1 of fang2019inference shows that if $r_N(\hat{\theta}_N - \theta_0) \Rightarrow \mathbb{G}_0$ in $\mathbb{D}$ for some $r_N \rightarrow \infty$ where $\mathbb{G}_0$ is tight, then $r_N\left[\phi(\hat{\theta}_N) - \phi(\theta_0)\right] \Rightarrow \phi'_{\theta_{0}}(\mathbb{G}_0)$. On top of this result, Theorem 3.2 of fang2019inference shows that $\phi'_{\theta_{0}}(\mathbb{G}_0)$ can be approximated by $\phi'_{N}(\mathbb{Z}_N^*)$ where $\mathbb{Z}_N^*$ and $\phi_{N}'$ are consistent estimators of $\mathbb{G}_0$ and $\phi'_{\theta_{0}}$, respectively. hong2018numerical suggest the NDM which numerically approximates $\phi'_{\theta_{0}}(\mathbb{G}_0)$ hence does not need analytical calculation of $\phi'_{\theta_0}$. hong2018numerical suggest the estimation of $\phi'_{\theta_0}$ by using the following numerical derivatives

where $\epsilon_N$ is a step size satisfying $\epsilon_N \rightarrow 0$ and $r_N\epsilon_N \rightarrow \infty$ as $N \rightarrow \infty$.

The NDM can be applied to the SD testing problem. Consider the testing with the $L_1$-type test statistic and our hypothesis of interest is the same as (ref). First, denote

For any $\theta \in l^{\infty}(\mathcal{X}) \times l^{\infty}(\mathcal{X})$, $\phi : l^{\infty}(\mathcal{X}) \times l^{\infty}(\mathcal{X}) \rightarrow \mathbb{R}$ is defined as

We then may write

where $ \theta_0 = (F_{1}^{(s)}(x), F_{2}^{(s)}(x)) \in l^{\infty}(\mathcal{X}) \times l^{\infty}(\mathcal{X})$. It can be seen that $\phi(\cdot)$ is Hadamard directionally differentiable with its derivative given as

where

Note that $\mathcal{C}_{0}(\theta_0)$ is equivalent to the contact set $\mathcal{C}_0$ in (ref). In addition, we have $\mathcal{C}_{+}(\theta_0) = \emptyset$ under the null hypothesis.

Letting $\hat{\theta}_N = (\Bar{F}_{1}^{(s)}(\cdot),\Bar{F}_{2}^{(s)}(\cdot))$, we may write the $L_1$-type test statistic as $r_N\phi(\hat{\theta}_N)$, where $r_N =\sqrt{N}$. In addition, we know that $r_N(\Bar{F}_{k}^{(s)}(\cdot) - F_{k}^{(s)}(\cdot))$ weakly converges to a tight Gaussian process. If we let $\mathbb{Z}_N^* = r_N(\Bar{F}_{k}^{(s),*}(\cdot) - \Bar{F}_{k}^{(s)}(\cdot))$, where $\Bar{F}_{k}^{(s),*}$ is the bootstrap sample analogue of $\Bar{F}_{k}^{(s)}$, then we can approximate the limit distribution of $r_N\phi(\hat{\theta}_N)$ by employing (ref).

There are some cases in which the first-order derivatives may degenerate, such as the the $L_2$-type test statistic. In this case, we can employ the second-order NDM hong2018numerical. $\phi_{\theta_0}''(h)$ is said to be the second order Hadamard directional derivative at $\theta_0 \in \mathbb{D}_{\phi}$ in the direction $h$ if

for all sequences $\{h_N\} \subset \mathbb{D}$ and $\{t_N\} \subset \mathbb{R}_{+}$ such that $t_N \downarrow 0$, $h_N \rightarrow h \in \mathbb{D}_0$ as $N \rightarrow \infty$ and $\theta + t_N h_N \in \mathbb{D}_{\phi}$ for all $N$. Theorem 4.1 of hong2018numerical shows that if $r_N(\hat{\theta}_n - \theta_0) \Rightarrow \mathbb{G}_0$ for some tight process $\mathbb{G}_0$, then

When the first-order derivatives degenerate, $\phi'_{\theta_0}(h) = 0$, it is possible to approximate $\frac{1}{2}\phi''_{\theta_0}(\mathbb{G}_0)$ by

where $\epsilon_N$ satisfies $\epsilon_N \rightarrow 0$ and $r_N\epsilon_N \rightarrow \infty$ as $N \rightarrow \infty$. In the case of the $L_{2}$-type test statistic, which is the simplified version of $LSW_N$ in (ref), we may take

Then we have $\phi_{\theta_0}'(h) = 0$ under the null hypothesis but

holds. Thus, similarly to the first-order NDM case, we can approximate the limit distribution of $r_N^2 \phi(\hat{\theta}_N)$ using (ref).

The Package: PySDTest

Installation

It is required to install the package, PySDTest, for both Python and Stata environments. The package is listed on the Python Package Index (PyPI) which is a repository of softwares for the Python programming language. If Python is installed on your computer, a user can install PySDTest by entering:

in Windows cmd or Mac (or Linux) terminal.\footnote{Alternatively, users can access source codes of the package through \href{https://github.com/lee-kyungho/pysdtest}{https://github.com/lee-kyungho/pysdtest}} For details about installing a Python package, see \href{https://packaging.python.org/tutorials/installing-packages/}{Python Package User Guide in PyPI}.

The package is written based on Python numpy module with version $\geq 1.26.4$. Therfore, it is also required to install numpy for employing the package.

To install the Stata module, pysdtest, users can use the following command in Stata:

Implementation in Python

This section describes how to use the package PySDTest in the Python environment. Notably, Stata with version $\geq$ 16.0 supports running the package within Stata. For more details on using Python in Stata, please refer to the \href{https://www.stata.com/stata16/python-integration/}{Python integration feature} in Stata.

First, users need to import the package through the following code in the Python environment in order to use the package:

PySDTest implements comprehensive SD tests developed in the literature. In addition, users can freely choose procedures suitable to their objects of study. In this section, we describe the details of how to implement each testing procedure.

We first generate simulated data as:

BD and LMW tests by pysdtest.test_sd()

The function pysdtest.test_sd() enables the user to implement the tests proposed by BD and LMW. The structure of the function pysdtest.test_sd is given as:

The following lines of codes implement the testing procedure by BD:

Users can run the subsampling procedure by giving subsampling as an input for the resampling argument:

Contact Set Approach by pysdtest.test_sd_contact()

The function pysdtest.test_sd_contact() in PySDTest implements the contact set approach proposed by LSW.\footnote{To be specific, the case when there is no estimated parameters and $q(x) = 1$ for all $x$ is implemented by pysdtest.test_sd_contact().} The structure of the function is written as

Arguments sample1, sample2, ngrid, s, resampling, b1, \texttt{b2}, \texttt{nboot}, and \texttt{quiet} are the same as \textsf{pysdtest.test_sd()}.

We need to take a sequence $c_N$ for the contact set estimation in the step 6 of LSW. Following the recommendation by LSW, we set

where $c$ is some positive constant. If two sample sizes are different, $N = (N_1 + N_2)/2$ holds. The argument c designates the constant $c$. We set the default value of $c$ to be $0.75$.

Users can implement the LSW test by:

Selective Recentering Approach by pysdtest.test_sd_SR()

The selective recentering approach of DH with $k=recenter$ (Equation ((ref))) can be implemented by pysdtest.test_sd_SR(). The structure of the function is

Arguments sample1, sample2, ngrid, s, resampling, b1, \texttt{b2}, \texttt{nboot}, and \texttt{quiet} are the same as \textsf{pysdtest.test_sd()}.

The selective recentering approach additionally requires two tuning parameters: $a_N$ and $\eta$. We follow the recommendation by DH and set $a_N$ to be

where $a$ is some positive constant and $N = N_1 + N_2$. Arguments a and eta set the values of $a$ and $\eta$ with the default given by $a = 0.1$ and $\eta = 10^{-6}$.

Numerical Delta Method by pysdtest.test_sd_NDM()

Users can approximate the limit distribution based on the NDM approach using the function pysdtest.test_sd_NDM(). The function has the following structure:

Three different types of functional are available for users: Kolmogorov-Smirnov (supremum), $L_{1}$, and $L_{2}$ type functionals. This can be set by giving input `KS', `L1', and `L2' for the `form' argument, respectively. We apply the second-order NDM for $L_2$ type testing and the limit distribution is approximated by $\frac{\phi(\hat{\theta}_N + \epsilon_N\mathbb{Z}_N^*)- \phi(\hat{\theta}_N) }{\epsilon^2_N}$.

In addition, the NDM requires users to specify the step size ($\epsilon_{N}$) as a tuning parameter. The user can set the step size by the argument `epsilon'. We set $\epsilon = r_{N}^{-1/16}$ as a default value when nothing is specified by the user. The user can implement the NDM for SD tests as

Table (ref) contains common arguments for all testing procedures: BD, LMW, LSW, DH, and NDM. In addition, we summarize additional arguments for each estimation in Table (ref).

Additional Features

Every implementation of testing through PySDTest saves results as an object hence users can save the test results. For example, if the user runs the following code

then the `results' in the above code is the variable in which the Python dictionary of the test result is allocated. The dictionary contains calculated test statistics, resampled test statstics, critical-value and p-value.

In addition, every function for testing provides a feature for plotting integrated empirical CDFs. This functionality may help users to visually compare testing objects. The plotting is based on the Python Matplotlib package. The method plot_CDF() has the following structure

In table (ref), arguments for plot_CDF() are listed.

For instance, plotting can be conducted by

and we show the corresponding result in Figure (ref).

Plotting a higher-order CDF can be conducted by

Figure (ref) shows the corresponding result.

PySDTest provides three resampling procedures: bootstrap, paired bootstrap, and subsampling. The subsampling() function generates and returns subsamples $\mathcal{S}_{i,b}$. The argument b means the length of block, and nsub reflects the number of resampled test statistics. In testing codes, we set nsub as $N-b+1$.

The bootstrap() function generates and returns bootstrap samples. The argument b means the number of observations for one boostrapping, which is used to calculate test statistics. In testing codes, we set b$=N$. The nboot means the number of bootstrapping.

The paired_bootstrap() function conducts bootstrapping in a pair $(X_{1,i},X_{2,i})$, and the description of arguments is the same as bootstrap().

PySDTest provides two functions that help users calculate the test statistics by themselves. The set_grid() function returns the `ngrid' number of grid points which are equally divided from the minimum to the maximum value of samples. For instance,

The CDF() calculates the (integrated) empirical CDF of the sample. For example,

Advanced Hypothesis Testing

Through PySDTest, the user can perform advanced hypothesis regarding SD. For example, consider the case when there are $K$ distributions. Then, the hypothesis of our interest can be $H_{0}^{(s)}:$ There exists at least one distribution that $s$-th order stochastically dominates some of the other distributions. $H_{1}^{(s)}$ is the negation of $H_{0}^{(s)}$. This can be expressed as

We describe testing procedures step-by-step and conduct testing via PySDTest. Before testing, we generate random samples.

First, we need to set options for the testing; the user should specify the SD order, the number of grid points, and grid. In this example, we test first order SD and set the number of grid points as 100.

Then, the user has to compute test statistics. This can be done by using the CDF() function in PySDTest. We consider the supremum-type test statistics, $T_{N}$, which may be written as

For simplicity, we let $K=3$. Then, the sample analogue of the test statistics can be computed as:

Then we conduct resampling. In here, we generate bootstrap samples as

The next step is to calculate bootstrap test statistics and attain bootstrap distribution. We conduct recentered bootstrap as

Then we compute the critical value under the significance level $\alpha$. Let $\alpha=0.05$, and

Lastly, we decide whether to reject the null hypothesis or not. If the computed test statistics is larger than the critical value, reject $H_{0}^{(s)}$.

Stata Command: pysdtest

In this section, we introduce the Stata command pysdtest. The command is based on the python package PySDTest. Therefore, Stata with version $\geq$ 16.0, Python 3 software and installation of the Python package are necessary for using the command. Supported testing procedures include, but are not limited to, tests through the least favorable case approach of BD, subsampling approach of LMW, contact set approach of LSW and selective recentering method of DH. In addition, the NDM can be used to compute the critical value. pysdtest also supports the combination of various resampling methods, test statistics and procedures for approximating the limiting distribution.

Our hypothesis of interest is (ref); for given stochastic order $s$, the null hypothesis $H_0^s$ is that $X_1$ (Sample 1) $s$-th order stochastic dominates $X_2$ (Sample 2) and the alternative hypothesis $H_1^s$ is the negation of $H_0^s$. pysdtest allows arbitrary $s$-th SD order for testing.

Syntax and Options

varlist specifies up to two variables for comparison. If only one variable is specified, a by() option must be used to divide the sample. The command supports several options to customize the testing procedure:

Stored Results

The following scalars, macros and matrices are stored as a result of running pysdtest:

Examples

In this section, we describe example usages of pysdtest in Stata. We first generate the simulated data from the standard normal distribution.

If there is no option given, pysdtest implements Kolmogorov-Smirnov type test statistics as in BD and LMW. The default resampling method is bootstrap and recentering as in BD.

The subsampling method as in LMW can be used by giving “subsampling" as an argument to the option resampling( ). In this case, the user needs to set the number of subsamples by b1( ) and b2( ) options.

In addition, users can employ the paired bootstrap method by giving “paired_bootstrap" as an argument to the option resampling( ). The number of observations for both sample1 and sample2 must be the same.

Also, the contact set approach can be utilized by giving “contact" to approach( ) option. The additional tuning parameter can be chosen by c( ) option. For the contact set approach, the integral type test statistic is employed as in (ref).

pysdtest provides by( ) option to compare distributions divdided by a binary variable. As an illustration, we use automobile data stored in Stata and compare price distributions by the foreign dummy.

If users want to switch the order of samples, they can use the option switch to swap the order of samples.

NDM can be used by giving “NDM" to approach( ) option. We employ Kolmogorov-Smirnov type test statistic in this example. We can also set the number of grid points by ngrid( ) option.

Empirical Illustration

In this section, we compare historical daily returns of Bitcoin price and S$\&$P 500 index and investigate SD relations between them. This paper applies PySDTest to illustrate how portfolio comparison based on SD criteria can be done. SD concept is useful for portfolio comparison because SD provides a general rule for the portfolio selection under certain expected utility paradigm. We can interpret testing results in terms of preference to prospects such as monotonicity and risk aversion of an economic agent. The replication code is shown in Appendix.

In this paper, daily returns, $R_{t}^{k}$, are defined as:

where $P_{t}^{k}$ is price of Bitcoin or S$\&$P 500 index at day $t$.\footnote{We use closing price of Bitcoin and S$\&$P 500. Because the cryptocurrency market does not actually `close', we alternatively use price at 11:59 p.m. provided by the \href{https://coinmarketcap.com/}{CoinMarketCap}} Table (ref) summarizes the time length, number of observations and descriptive statistics of daily return data. Because the cryptocurrency market operates in weekends and holidays, the number of observations of Bitcoin price is more than S&P 500 index.

Standard mean-variance analysis cannot determine superiority of daily returns between Bitcoin and S&P 500 index. Although the mean return of Bitcoin price is larger than that of S$\&$P 500 index, standard deviation of Bitcoin is also larger. In addition, the distribution of Bitcoin price seems to be more dispersed than that of S$\&$P 500 index. Figure (ref) shows cumulative distribution functions of daily return of Bitcoin and S$\&$P 500 index.\footnote{This is plotted by using the function plot_CDF() in PySDTest}

We conduct testing for the following null hypotheses to investigate SD relations:

for $s=1,2$ The alternative hypotheses $H_{1,i}^{s}$ are the negation of $H_{0,i}^{s}$ for $i,s=1,2$, respectively. We use the subsampling method by LMW in order to account for time-series dependence. We set 100 grid points. 1,000 and 900 are set as a size of subsamples for daily returns of Bitcoin and S&P 500 index, respectively. Table (ref) summarizes the testing results.

In Table (ref), all hypotheses of the first SD order ($s=1$) are rejected. This means that an economic agent with monotone increasing utility function strictly prefer neither Bitcoin nor S&P 500 to the other. In contrast, it is shown that daily returns of S&P 500 second order ($s=2$) stochastically dominate those of Bitcoin. This implies that an agent with risk-averse utility function prefers investing in S&P 500 index to Bitcoin.

Conclusion

In this paper, we have introduced a Python package, PySDTest, and a Stata command, pysdtest. The package implements various Stochastic Dominance (SD) tests, including those proposed by BD, LMW, LSW, and DH. PySDTest includes testing using the NDM, which applies modern developments in statistical inference methods to directionally differentiable functions. We provide guidance for using PySDTest and pysdtest, including testing for extended notions of the SD hypothesis. We applied PySDTest to compare the daily return distributions between Bitcoin prices and the S&P 500 index. The results show that the S&P 500 second-order stochastically dominates Bitcoin returns.

Empirical implementation of the SD concept has not been as extensive as it could be despite the fundamental role of the SD concept in diverse research areas and the recent advances in inference procedures. This may be partly due to the lack of publicly available software. In this paper, we provide a comprehensive Python package and Stata command for SD tests to assist empirical practitioners and academics. We look forward to seeing many applications of our software that help practitioners as intended.