A test is parametric if its goal is to test parameters of a known/unknown distribution.

Note: we call population the total data from which a sample is extracted. A statistical test aims at finding information about the distribution the sample is extracted from.

Procedure:

1) find the test to perform

2) find the right estimator to use

3) deduce the reject region

4) compute the test statistic

5) retrieve quantiles of known distributions

See the associated notebook for numerical examples.

Example 1 (Z-test):

$X_1,…,X_n~(iid)\sim \mathbb{P_\theta}$

We want to know the mean of an unknown distribution from which we have a sample.

Requirements/assumptions:

• we need the population to be normally distributed. Note: this may not be needed when sample size is high (> 30) as CLT will anyway give normality

• we need the standard deviation of the population

1) find the test to perform

2) find the right estimator to use

Since we are testing the mean, we choose the empirical mean as estimator $\widehat{\theta}=\frac{1}{n}\sum{X_i}$

3) deduce the reject region

We fix $k$ for a rejection level $\alpha$. The rejection region is \(Z= \{\widehat{\theta} \ge k\}\)

We look for $k$ defined as such:

$\mathbb{P}_{\theta \in \Theta_0}(\widehat{\theta} \ge k)=\alpha$ => under $\mathcal{H}_0$, we reject the hypothesis when our estimator $\widehat{\theta}$ is above $k$

Intuitively, we want to keep our hypothesis if it’s verified in most of the cases => under our hypothesis, there is a low probability that we are in the rejection region.

Thus, if in real life we have a result that makes the hypothesis unverified, we reject the hypothesis. However, we have a risk of $\alpha$ that our hypothesis was correct and that we ended up in the rejection region by mistake.

4) compute the test statistic

We center and reduce the estimator in order to get the Gaussian law and thus end up with known quantiles:

\(\mathbb{P}_{\theta = a}(T \ge \frac{\sqrt{n} (k-a)}{\sqrt{\sigma^2}})=\alpha\) with $T \sim_{n \to \infty} \mathcal{N}(0,1)$

$T$ is the test statistic (a test statistic is a random variable for which we know the law under $\mathcal{H}_0$)

5) retrieve quantiles of known distributions

Finally, $\frac{\sqrt{n} (k-a)}{\sqrt{\sigma^2}}=q_\alpha$ => we can find $k$ telling us when rejecting $\mathcal{H}_0$

Why not looking at the average directly?

=> the average can be influenced by the outliers and thus doesn’t take into consideration extreme events.

How about the median?

=> the median doesn’t take into account the distribution / tendency of the values.

The p-value is the lowest error probability we want to make when rejecting our hypothesis.

The lower the p-value is, the less error we make in rejecting our hypothesis so the more significant the rejection is.

When performing OLS, our hypothesis is $\theta_{x1}=0$ so we don’t reject it if the pvalue column is higher than our threshold. In the below OLS result, pvalues are displayed in column $P>|t|$. All variables are significant.

Example 2 (T-test): when the variance is not known.

Say we want to test whether a coefficient is zero:

1) find the test to perform

2) find the right estimator to use

$\widehat{\theta_j} = (X^TX)^{-1}X^TY$

3) deduce the reject region

$Z={k_1 \leq \widehat{\theta_j} \leq k_2}$

4) compute the test statistic

$T_j = \frac{\widehat{\theta_j}-\theta_j}{\sigma_{\theta_j}} = \frac{\widehat{\theta_j}}{\sigma_{\theta_j}} \sim \mathcal{N}(0,1)$ with $\sigma_{\theta_j} = \sigma \sqrt{(X^TX)^{-1}}$ (recall that $\sigma = \sigma_{\epsilon}$)

Since we don’t know $\sigma$, we can use the Cochrane theorem to remove this value:

$T_j = \frac{ \frac{\widehat{\theta_j}}{\sigma \sqrt{(X^TX)^{-1}}} \sim \mathcal{N}(0,1)}{\sqrt{\frac{\widehat{\sigma}^2 (n-p-1)}{\sigma^2} \sim \mathcal{X}_{n-p-1}}} \sim \mathcal{T}(n-p-1)$ with $\widehat{\sigma}^2 = \frac{1}{n-p-1} \Sigma \epsilon ^2$

$T_j = \frac{\widehat{\theta_j}}{\Sigma \epsilon ^2 \sqrt{(X^TX)^{-1}}}$

5) retrieve quantiles of known distributions

Finally,

$\mathcal{P}_{\theta_j=0}(\frac{k_1}{\Sigma \epsilon ^2 \sqrt{(X^TX)^{-1}}} \leq T_j \leq \frac{k_2}{\Sigma \epsilon ^2 \sqrt{(X^TX)^{-1}}}) = \alpha$

Thus, $\frac{k_1}{\Sigma \epsilon ^2 \sqrt{(X^TX)^{-1}}} = t_{\frac{\alpha}{2}}$ (same for $k_2$)

Example 3 (T-test with forward selection):

Concept: regress all variables one by one on the most significant variable’s residual, remove the most significant variable after each full round.

(x-axis is the order in which we selected variables)

We can then select only the most significant variables based on p-values on variables from list \(sel\_variables\)

Note: since \(pval = 2*(1-cdf(T)) = 2*\frac{1-(1-\alpha)}{2}\), choosing the biggest T-stat is equivalent to choose the smallest p-value

Example 4 (F-test):

When several variables are correlated (often the case in practice), the student test is not efficient enough since it does not take the correlation into account. F-test allows to test global significativity.

Let’s say we have 4 variables and we want to check the significativity of 2 of them.

$SSR = sum~squared~residuals = \Sigma (\widehat{y_i} - y_i)^2$

$F = \frac{(SSR_C - SSR_{NC})/(p_{NC} - p_C)}{(SSR_{NC})/(n-p_{NC})} \sim \mathcal{F}(p_{NC} - p_C, n-p_{NC})$

NC: not constraint model

C: constraint model

Method:

• OLS on not constraint model => computation of $SSR_{NC}$

• OLS on constraint model => computation of $SSR_{C}$

• Computation of the Fisher stat => computation of p-value (using complementary cdf as above)

# Non constraint model
X0=np.column_stack((educ, exper, tenure, const))
model=sm.OLS(y,X0)
results = model.fit()
u=results.resid
SSR0=u.T@u

# Constraint model
X=np.column_stack((const, educ, tenure))
model=sm.OLS(y,X)
results = model.fit()
u=results.resid
SSR1=u.T@u

# Computation of Fisher stat
n=np.shape(X0)[0]
F=((SSR1-SSR0)/1)/(SSR0/(n-4))
f.sf(F,1,n-4) # p-value

Example 5 (AB-test):

A/B tests are typically used in digital marketing. Two groups of people are randomly selected:

People in group A (control group) see version A and people in group B (test group) see version B. The assessed metric is the conversion rate i.e. the proportion of users who clicked on the “Buy” button.

We can define two random variables $X_A$ and $X_B$ to modelize the clicks. Those variables are binary: 0 (no click) or 1 (click).

If we repeat the experiment many times (with different samples), the two variables follow a Binomial distribution:

Below is an example of results following the experiment:

The two hypothesis are:

Where $p_A$, $p_B$ are the true conversion rates.

Since the variables are binary, we can estimate the conversion rates with the empirical mean:

Our final estimator (according to $\mathcal{H}_0$) is thus $\hat p_A-\hat p_B$. Thanks to the CLT,

Thus,

Note: $p_A$ and $p_B$ in the variance are not known. We may use estimated values (to be verified).

We can now plot the estimator under both hypothesis:

Significance level: probability of false positives i.e. the probability of wrongly rejecting the null hypothesis. Using $\alpha$ as a significance level, $\mathbb{P}_{\mathcal{H}_0}(\hat d >k)= \alpha/2$.

Thus $k = q_{1-\frac{\alpha}{2}}\sigma_{\hat d} + 0$. We can deduce the confidence interval $[-q_{1-\frac{\alpha}{2}}\sigma_{\hat d}, q_{1-\frac{\alpha}{2}}\sigma_{\hat d}]$

Power of the test: probability of true positive i.e. the probability of correctly rejecting the null hypothesis. Intuitively, the more the two shapes are separated, the more true positives we have.

Note: increasing the number of observations has the effect of reducing the standard deviation of both estimators ($\mathcal{H}_0$ and $\mathcal{H}_1$) and thus add more space between the two shapes => increase the power of the test.

Example 6 (Z-test (other explanation)):

Procedure:

1) Define the test to perform

2) Find the right estimator to use

3) Bring the estimator to a known distribution

4) See how likely is our real-life example based on the known distribution

1) Define the test to perform.

2) Find the right estimator to use. Since we are testing the mean, we choose the empirical mean as estimator $\widehat{\theta}=\frac{1}{n}\sum{X_i}$

3) Bring the estimator to a known distribution. Thanks to CLT, we know that the estimator follows a normal distribution.

In particular, under $\mathcal{H}_0$, $\widehat{\theta} \sim \mathcal{N}(a, \frac{\sigma}{\sqrt{n}})$ where $\sigma$ is the standard deviation on the population (known).

3) See how likely is our real-life example based on the known distribution.

We can see that the result we got “in real life” is very unlikely to happen if we repeat the experiment under $\mathcal{H}_0$ many times. Thus, we can reject $\mathcal{H}_0$. However, there’s still a chance (<5%) that $\mathcal{H}_0$ is true and that our real life example is just an extreme event.

Example 7 (Z-test (proportions)):

(i) We believe there are the same number of women than men in the population. Our sample has $n=100$ people including $k=45$ females. Is our sample representative enough?

(ii) We toss a coin $n=50$ times and get $k=30$ heads. Is the coin fair?

1) Define the test to perform. $\mathcal{H}_0: p=p_0$

2) Find the right estimator to use. \(\hat p = \frac{n_{success}}{n}\)

3) Bring the estimator to a known distribution. $X \sim \mathcal{B}(n,p)$ with $X$ the number of successes. Thanks to normal approximation, we have (under $\mathcal{H}_0$):

$X \sim \mathcal{N}(np_0, \sqrt{np_0(1-p_0)})$ with $X$ the number of successes.

Thus: $\hat p \sim \mathcal{N}(p_0, \sqrt{\frac{p_0(1-p_0)}{n}})$

4) See how likely is our real-life example based on the known distribution (illustration of example (i) above).

We can see that the result we got “in real life” is likely to happen if we repeat the experiment under $\mathcal{H}_0$ many times. Thus, we can’t reject $\mathcal{H}_0$.

The rejection region is outside the confidence interval. For a level of $95\%$, we don’t reject $\mathcal{H}_0$ if $\hat p \in [p_0 \pm 1.96 \sqrt{\frac{p_0(1-p_0)}{n}}]$.

Note: if $n$ is too small, the normal approximation is not accurate enough. One can use the binomial test instead (if $p_0 \neq 0.5$ the easiest way is to use “Equal distance from expected” described here). Another option is to use an “abaque” as explained here. See example 8 for an illustration on a toy example.

Notes about confidence intervals:

• a confidence interval is always around an estimator.

• the size of the confidence interval decreases with the number of observations. Intuitively, when dealing with a lower number of observations, we need to have more extreme results to safely reject $\mathcal{H}_0$.

Example 8 (Binomial test):

(i) My trading strategy placed 20 trades, 15 were successful. Can I say that my results are statistically significant?

(ii) I met 4 French people, they were all arrogant. Can I deduce that all French people are arrogant?

1) Define the test to perform: $\mathcal{H}_0: p=50\%$ where $p$ is the true proportion of trades that were successful.

2) Find the right estimator to use: $\hat p = \frac{n_{success}}{n}$.

3) Bring the estimator to a known distribution: $X = n \hat p$, $X \sim \mathcal{B}(n,p)$. Under $\mathcal{H}_0$, $X \sim \mathcal{B}(20, 50\%)$.

4) Compute the p-value, using the distribution of $X$ under $\mathcal{H}_0$.

• One-sided test ($\mathcal{H}_1: p > 50\%$, i.e. “my result is better than random”):

$pval = \mathbb{P}_{\mathcal{H}_0}(X \geq 15) = \sum_{k=15}^{20} C_{20}^k(50\%)^k(1-50\%)^{20-k} = \frac{\sum_{k=15}^{20} C_{20}^k}{2^{20}} = \frac{21700}{1048576} \approx 2.07\%$

• Two-sided test ($\mathcal{H}_1: p \neq 50\%$), i.e. “my result is different from random”:

$pval = \mathbb{P}_{\mathcal{H}_0}(X \geq 15) + \mathbb{P}_{\mathcal{H}_0}(X \leq 5) = 2.07\% + 2.07\% = 4.14\%$ (5 successes out of 20 is as far from 10 as 15 is)

In both cases $pval < 5\%$, so we reject $\mathcal{H}_0$ at the 5% level: the results are significant.

The green lines are the 95% confidence interval under $\mathcal{H}_0$.

Unlike $\alpha$, the power depends mostly on the sample size: the larger $n$, the smaller the edge the test can detect. The power also depends on the true success rate, which we don’t know. The table and the curves below give the power for any true success rate (one-sided test at 5%), and the number of trades needed to reach a power of 80%:

In other words, with 20 trades, the test only reliably detects very profitable strategies (a success rate of 80% or more). A realistic edge for a trading strategy (55–60%) needs hundreds of trades.

The script below computes the p-value, the rejection region, the minimum $n$ needed for the observed success rate to be significant and the power table above. Change n_success, n and p_table at the top; it also draws the two figures above.

import math
import numpy as np
import matplotlib.pyplot as plt
import scipy.stats as scs

# inputs
n_success = 15  # observed successes
n = 20          # number of trials
p0 = 0.5        # success rate under H0 (random)
p_table = [0.55, 0.6, 0.65, 0.7, 0.75, 0.8, 0.9]  # true success rates in the power table
alpha = 0.05    # significance level
target_power = 0.8  # power we want to reach
two_sided = True  # figure 1: show the two-sided test (acceptance interval, both tails)
# the power and the minimum n are computed for the one-sided test, H1: p > p0


def critical_value(n, p0, alpha):
    # smallest k such that P_H0(X >= k) <= alpha; rejection region is X >= k
    # (k = n + 1 means that no result can be significant)
    k = np.arange(0, n + 2)
    return int(k[scs.binom(n, p0).sf(k - 1) <= alpha][0])


def acceptance_interval(n, p0, alpha):
    # two-sided test: H0 is not rejected if low <= X <= high, each tail outside has
    # probability <= alpha / 2 under H0
    k = np.arange(0, n + 1)
    dist = scs.binom(n, p0)
    lower_tail = k[dist.cdf(k) <= alpha / 2]       # P_H0(X <= k) <= alpha / 2
    upper_tail = k[dist.sf(k - 1) <= alpha / 2]    # P_H0(X >= k) <= alpha / 2
    low = int(lower_tail[-1]) + 1 if len(lower_tail) else 0
    high = int(upper_tail[0]) - 1 if len(upper_tail) else n
    return low, high


def p_value(n_success, n, p0):
    return scs.binom(n, p0).sf(n_success - 1)  # P_H0(X >= n_success)


def power(n, p0, p_true, alpha):
    return scs.binom(n, p_true).sf(critical_value(n, p0, alpha) - 1)  # P_p(X >= k)


def min_n_significant(success_rate, p0, alpha, n_max=10_000):
    # smallest n for which the same success rate is significant
    # (number of successes rounded down, so the rate never exceeds the observed one)
    for m in range(1, n_max + 1):
        if p_value(math.floor(success_rate * m), m, p0) <= alpha:
            return m
    return None


def min_n_power(p_true, p0, alpha, target_power, n_max=2_000):
    # smallest n for which the power reaches target_power, for each true success rate
    # (the power is not monotonic in n because X is discrete: this is the first n reaching it)
    ns = np.arange(1, n_max + 1)
    k_crits = np.array([critical_value(m, p0, alpha) for m in ns])
    n_req = []
    for p in p_true:
        reached = ns[scs.binom(ns, p).sf(k_crits - 1) >= target_power]
        n_req.append(reached[0] if len(reached) else np.nan)
    return np.array(n_req, dtype=float)


k_crit = critical_value(n, p0, alpha)
pval = p_value(n_success, n, p0)
pval_two_sided = scs.binomtest(n_success, n, p0).pvalue
n_min = min_n_significant(n_success / n, p0, alpha)

print(f'observed success rate: {n_success}/{n} = {n_success / n:.1%}')
decision = 'reject' if pval <= alpha else 'do not reject'
print(f'p-value (one-sided):   {pval:.2%}  -> {decision} H0 at {alpha:.0%}')
print(f'p-value (two-sided):   {pval_two_sided:.2%}')
low, high = acceptance_interval(n, p0, alpha)
print(f'acceptance interval (two-sided): {low} <= X <= {high},'
      f' i.e. {low / n:.1%} <= success rate <= {high / n:.1%}')
if k_crit <= n:
    size = scs.binom(n, p0).sf(k_crit - 1)  # actual type I error, <= alpha
    print(f'rejection region:      X >= {k_crit}, i.e. success rate >= {k_crit / n:.1%}'
          f' (actual size {size:.2%})')
else:
    print(f'rejection region:      empty, no result out of {n} is significant at {alpha:.0%}')
print(f'min n to be significant at a {n_success / n:.1%} success rate: {n_min}')

# power table: power with n trials, and n needed to reach target_power
n_needed = min_n_power(p_table, p0, alpha, target_power)
print(f'\ntrue success rate | power with n={n} | n needed for a power of {target_power:.0%}')
for p, m in zip(p_table, n_needed):
    print(f'{p:17.0%} | {power(n, p0, p, alpha):15.1%} | {m:.0f}')

x = np.arange(0, n + 1)
pmf_null = scs.binom(n, p0).pmf(x)

# figure 1: distribution under H0, observed value and p-value
fig, ax = plt.subplots(figsize=(12, 6))
ax.bar(x, pmf_null, alpha=0.5, label='$X$ under $H_0$')
if two_sided:
    # results at least as extreme = as unlikely under H0 as the observed one (as scs.binomtest)
    extreme = pmf_null <= pmf_null[n_success] * (1 + 1e-7)
    ax.bar(x[extreme], pmf_null[extreme], color='red', alpha=0.5,
           label=f'two-sided p-value = {pval_two_sided:.2%}')
    low, high = acceptance_interval(n, p0, alpha)
    ax.axvline(x=low - 0.5, c='green', linestyle='-')
    ax.axvline(x=high + 0.5, c='green', linestyle='-',
               label=f'{1 - alpha:.0%} acceptance interval: {low} <= X <= {high}')
else:
    ax.bar(x[x >= n_success], pmf_null[x >= n_success], color='red', alpha=0.5,
           label=f'$P_(X \\geq {n_success})$ = p-value = {pval:.2%}')
ax.axvline(x=n_success, c='red', linestyle='-', label=f'observed: {n_success} successes')
if n <= 20:
    ax.set_xticks(x)
ax.set_xlabel('number of successes $X$')
ax.legend()
fig.savefig('binomial-test.png', bbox_inches='tight')

# figure 2: power curves and number of trials needed, for any true success rate
p_true = np.arange(p0, 0.951, 0.01)
fig, (ax1, ax2) = plt.subplots(1, 2, figsize=(14, 6))
for m in sorted({n, 50, 100, 200}):
    ax1.plot(p_true, [power(m, p0, p, alpha) for p in p_true],
             label=f'n={m}', lw=3 if m == n else 1.5)
ax1.axhline(y=target_power, c='grey', linestyle='--', label=f'power={target_power:.0%}')
ax1.axhline(y=alpha, c='grey', linestyle=':', label=f'alpha={alpha:.0%}')
ax1.set_xlabel('true success rate $p$')
ax1.set_ylabel('power')
ax1.set_title('power of the one-sided test')
ax1.legend()

p_req = p_true[p_true >= p0 + 0.05]
ax2.plot(p_req, min_n_power(p_req, p0, alpha, target_power), marker='o', ms=3)
ax2.axhline(y=n, c='grey', linestyle='--', label=f'n={n}')
ax2.set_yscale('log')
ax2.set_xlabel('true success rate $p$')
ax2.set_ylabel('number of trials $n$ (log scale)')
ax2.set_title(f'number of trials needed for a power of {target_power:.0%}')
ax2.legend()
fig.savefig('binomial-test-power-curve.png', bbox_inches='tight')

plt.show()