   216    The current implementation of :func:`lwdid.lwdid` stores summary
   217    randomization inference information (``ri_pvalue``, ``ri_method``,
   218    ``rireps``, ``ri_valid``, ``ri_failed``, ``ri_seed``) but does not
   219    retain the full array of permuted test statistics. If you need to
   220    visualize the permutation distribution, you must construct and store
   221    the permuted statistics yourself in a custom resampling loop and then
   222    pass that array to standard plotting tools such as :mod:`matplotlib`.
   214 .. note::
   215 
Randomization Module (randomization)
=====================================

The randomization module implements randomization inference (RI) for
non-parametric hypothesis testing in difference-in-differences estimation.

.. automodule:: lwdid.randomization
   :members:
   :undoc-members:
   :show-inheritance:

Overview
--------

Randomization inference provides an alternative to t-based inference that does
not rely on distributional assumptions (normality, homoskedasticity). It is
particularly useful when:

- Sample size is very small (N < 10)
- Normality is questionable
- You want a robustness check for parametric inference

The module implements two RI methods:

- **Permutation** (recommended): Classical Fisher randomization inference
  without replacement
- **Bootstrap**: With-replacement sampling (for backward compatibility)

How Randomization Inference Works
----------------------------------

Basic Procedure
~~~~~~~~~~~~~~~

1. **Compute observed test statistic** (e.g., t-statistic for ATT)
2. **Permute treatment assignment** across units, keeping the number of
   treated units fixed
3. **Re-estimate the model** with permuted treatment

4. **Record the permuted test statistic**
5. **Repeat steps 2-4** many times (e.g., 1000 permutations)
6. **Compute p-value** as the fraction of permuted statistics that are more
   extreme than the observed statistic

Null Hypothesis
~~~~~~~~~~~~~~~

RI tests the **sharp null hypothesis**: treatment has no effect on any unit.

Under this null, the observed outcomes are the same as the potential outcomes
under control, so permuting treatment assignment should not systematically
change the test statistic.

P-Value Calculation
~~~~~~~~~~~~~~~~~~~

For a two-sided test:

p-value = (# of permuted t-statistics ≥ observed t-statistic) / (# of permutations)

The p-value represents the probability of observing a test statistic as
extreme as the observed one under random treatment assignment.

RI Methods
----------

Permutation (Recommended)
~~~~~~~~~~~~~~~~~~~~~~~~~

**Description:** Sample treatment assignments without replacement. Each

permutation is a unique reassignment of the treatment indicator across units.

**Advantages:**

- Theoretically correct for Fisher's exact test
- No duplicate permutations (until all possible permutations are exhausted)
- Recommended in the literature

**Usage:**

.. code-block:: python

   from lwdid import lwdid

   results = lwdid(
       data, 'y', 'd', 'unit', 'year', 'post', 'demean',
       ri=True,
       ri_method='permutation',  # Recommended
       rireps=1000,
       seed=42
   )

Bootstrap
~~~~~~~~~

**Description:** Sample treatment assignments with replacement. The same

permutation can appear multiple times.

**Advantages:**

- Simpler implementation
- Can generate more than C(N, N_treated) permutations

**Disadvantages:**

- Not the classical Fisher RI
- May have duplicate permutations

**Usage:**

.. code-block:: python

   results = lwdid(
       data, 'y', 'd', 'unit', 'year', 'post', 'demean',
       ri=True,
       ri_method='bootstrap',
       rireps=1000,
       seed=42
   )

**Note:** Included primarily for backward compatibility. Use ``'permutation'``

for new analyses.

Practical Considerations
------------------------

Number of Permutations
~~~~~~~~~~~~~~~~~~~~~~

**Recommended values:**

- **1000**: Good balance of precision and speed (default)
- **2000-5000**: More precise p-values for publication
- **10000+**: High precision, but computationally expensive

**Precision of p-value:**

With R permutations, the standard error of the p-value is approximately:

SE(p) ≈ √(p(1-p) / R)

For p = 0.05 and R = 1000: SE ≈ 0.007

Maximum Permutations
~~~~~~~~~~~~~~~~~~~~

The maximum number of unique permutations is:

C(N, N_treated) = N! / (N_treated! × N_control!)

For example:

- N = 10, N_treated = 2: C(10,2) = 45 permutations
- N = 20, N_treated = 5: C(20,5) = 15,504 permutations
- N = 39, N_treated = 1: C(39,1) = 39 permutations

If ``rireps`` exceeds the maximum, permutation method will sample with
replacement after exhausting unique permutations.

Computational Cost
~~~~~~~~~~~~~~~~~~

RI requires re-estimating the model for each permutation, so computational
time scales linearly with ``rireps``:

- 1000 permutations: ~1-10 seconds (typical)
- 10000 permutations: ~10-100 seconds

For large datasets or many permutations, consider:

- Using a smaller ``rireps`` for exploratory analysis
- Running RI only for final results
- Using parallel processing (not currently implemented)

Reproducibility
~~~~~~~~~~~~~~~

Always set a ``seed`` for reproducible results:

.. code-block:: python

   results = lwdid(..., ri=True, seed=42)

Without a seed, results will vary slightly across runs due to random
permutation sampling.

Interpreting RI Results
-----------------------

Comparing with t-based Inference
~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. code-block:: python

   results = lwdid(
       data, 'y', 'd', 'unit', 'year', 'post', 'demean',
       vce='hc3', ri=True, rireps=2000, seed=42
   )

   print(f"t-based p-value: {results.pvalue:.4f}")
   print(f"RI p-value: {results.ri_pvalue:.4f}")

**Interpretation:**

- **Similar p-values:** Parametric assumptions (normality, homoskedasticity)
  are likely reasonable
- **RI p-value > t-based p-value:** Parametric inference may be too liberal
  (rejecting too often)
- **RI p-value < t-based p-value:** Parametric inference may be too
  conservative

Visualizing RI Distribution
~~~~~~~~~~~~~~~~~~~~~~~~~~~~

.. note::
   The current implementation of :func:`lwdid.lwdid` only stores summary
   randomization inference information (``ri_pvalue``, ``ri_method``,
   ``rireps``, ``ri_valid``, ``ri_failed``, ``ri_seed``) but does not retain
   the full array of permuted test statistics. If you need to visualize the
   permutation distribution, you must construct and store the permuted
   statistics yourself in a custom resampling loop and then pass that array
   to standard plotting tools such as :mod:`matplotlib`.
  well-behaved

Limitations
-----------

Sharp Null Hypothesis
~~~~~~~~~~~~~~~~~~~~~~

RI tests the sharp null (no effect on any unit), which is stronger than the
average null (no average effect). If treatment effects are heterogeneous,
RI may reject even when the average effect is zero.

Computational Cost
~~~~~~~~~~~~~~~~~~

RI is computationally intensive for large ``rireps``. For very large datasets
or complex models, RI may be impractical.

No Confidence Intervals
~~~~~~~~~~~~~~~~~~~~~~~~

Standard RI provides p-values but not confidence intervals. Constructing RI
confidence intervals requires inversion of the test, which is not currently
implemented.

Examples
--------

Basic RI
~~~~~~~~

.. code-block:: python

   from lwdid import lwdid
   import pandas as pd

   data = pd.read_csv('data.csv')

   results = lwdid(
       data, 'outcome', 'treated', 'unit', 'year', 'post', 'demean',
       ri=True, rireps=1000, seed=42
   )

   print(f"ATT: {results.att:.4f}")
   print(f"t-based p-value: {results.pvalue:.4f}")
   print(f"RI p-value: {results.ri_pvalue:.4f}")

RI with Robust SEs
~~~~~~~~~~~~~~~~~~

Combine RI with heteroskedasticity-robust standard errors:

.. code-block:: python

   results = lwdid(
       data, 'outcome', 'treated', 'unit', 'year', 'post', 'detrend',
       vce='hc3',
       ri=True, rireps=2000, ri_method='permutation', seed=42
   )

   print(f"HC3 p-value: {results.pvalue:.4f}")
   print(f"RI p-value: {results.ri_pvalue:.4f}")

High-Precision RI
~~~~~~~~~~~~~~~~~

For publication, use more permutations:

.. code-block:: python

   results = lwdid(
       data, 'outcome', 'treated', 'unit', 'year', 'post', 'demean',
       ri=True, rireps=10000, ri_method='permutation', seed=12345
   )

   print(f"RI p-value (10,000 permutations): {results.ri_pvalue:.5f}")

See Also
--------

- :func:`lwdid.lwdid` - Main estimation function with RI support
- :class:`lwdid.LWDIDResults` - Results object containing RI results
- :doc:`../user_guide` - Comprehensive usage guide
- :doc:`../methodological_notes` - Theoretical background
