where \(Y_{it}\) and \(W_{it}\) are binary. The bitmap implementation stores panel histories in fixed-width integer words and computes \(X'X\) and \(X'Y\) from unit and time margins.
The computational reduction is:
\[
\text{full FE regression}
\rightarrow
\text{Mundlak sufficient statistics}
\rightarrow
\text{direct moment computation}
\rightarrow
\text{direct bitmap moment computation}.
\]
The Mundlak normal equations reduce the regression to a small set of counts and intersections. For binary panels those counts can be computed with bitwise AND and popcount over packed machine words.
The public estimator API is:
from duckreg import BitmapMundlak, BitmapMundlakEventStudystatic_fit = BitmapMundlak.from_long(data).fit()event_fit = BitmapMundlakEventStudy.from_long(data).fit()
The from_long(...) constructors are convenience wrappers: each call packs the long panel before fitting. When the same panel feeds multiple estimators or many bootstrap draws, pack once with LongPanelBitPacker().fit_transform(data) and pass the resulting BitmapPanel to BitmapMundlak and BitmapMundlakEventStudy. Bootstrap replications then resample unit multiplicities against the cached layout; they do not pay the long-panel packing cost again.
The lower-level packing classes and display helpers live in duckreg.bitmap_utils. They are used directly below to show the physical layout and to reuse a packed panel across repeated moment queries.
A bitmap stores a Boolean vector as bits inside an integer word. For example, a treatment path such as
W_i = [1, 0, 1, 1, 0, 0, 1, 0]
can be stored in one byte, and for \(T \le 64\) it fits in one uint64. The operation
\[
\sum_t W_{it}Y_{it}
\]
is then computed as
\[
\operatorname{popcount}(W_i \& Y_i),
\]
where & keeps only positions where both variables are one and popcount counts the one bits. One AND plus one popcount replaces up to 64 scalar row-level comparisons and additions.
The binary case stores each unit’s panel history by concatenating the time series into fixed-width words. For unit \(i\), bit position \(t \bmod 64\) in chunk \(\lfloor t / 64 \rfloor\) stores the value of \(W_{it}\) or a bitplane for \(Y_{it}\). If \(T \le 64\), each binary variable is one word per unit. If \(T = 365\), each binary variable is six words per unit.
The unit-major table has \(N \lceil T / 64 \rceil\) rows of words. The logical panel still contains \(NT\) scalar observations, and the physical representation stores those observations as packed words:
\[
\text{bytes per binary variable in unit-major layout}
\approx
8N \lceil T / 64 \rceil.
\]
For models with time margins, the implementation also keeps a time-major transpose. That table stores bit position \(i \bmod 64\) in unit chunk \(\lfloor i / 64 \rfloor\) for each time period. With binary \(Y\) and binary \(W\), keeping both unit-major and time-major layouts costs approximately
\[
8 \cdot 2
\left(
N \lceil T / 64 \rceil
+ T \lceil N / 64 \rceil
\right)
\approx
NT / 2
\text{ bytes},
\]
up to identifier columns and table overhead. At \(N = 500{,}000\) and \(T = 30\), this is only a few megabytes for the raw words. Layout construction is a one-time cost when the packed representation is reused across point estimates, event studies, and bootstrap draws. The benchmark below reports layout construction separately from estimator runtime.
Horizontal And Vertical Scans
A row scan loops over every logical observation:
\[
O(NT)
\]
scalar row operations. A unit-major bitmap scan loops over units and 64-period chunks:
machine-word operations. For \(T \le 64\), many unit-level statistics become \(O(N)\) word scans. With \(N = 500{,}000\) and \(T = 30\), a long-table scan touches \(15{,}000{,}000\) rows, while the unit-major bitmap scan touches \(500{,}000\) words. The logical scan count falls by a factor of 30; each word has unused bits because the panel has fewer than 64 periods. With \(T = 365\), the reduction is \(365 / \lceil 365 / 64 \rceil \approx 60.8\).
Time margins use the transposed layout. Unit-major packing is good for quantities such as
The estimator works from a small list of marginal counts and intersections. Those counts are the native operations of the packed layout.
Low-cardinality outcomes extend the same way. If \(Y_{it} \in \{0, \dots, K\}\), store a bit mask for each outcome level, or for each nonzero level when zero can be implicit. For any subset \(A\) represented by a treatment, cohort, event-time, or design mask,
The bitplanes also give \(\sum_{(i,t) \in A} Y_{it}^2\) by replacing \(k\) with \(k^2\). Binary outcomes are the one-bit special case. Small integer supports such as \(0,\dots,9\) require a modest number of masks. Exact real-valued outcomes are outside this bitmap-only moment strategy unless they are quantized first.
The asymptotic class is still linear in the logical panel size. The gain comes from changing the unit of work from scalar rows to machine words, while also reducing memory traffic, branches, and interpreter overhead. The bitmap kernel is most attractive when the layout is already available or amortized over repeated queries; packing a long table first is itself an \(O(NT)\) operation.
Five-Period Example
Start with a tiny long panel with five time periods. The bitmap convention here stores time period 0 in bit position 0. In the display table below, the leftmost character is period 0, so 00111 means untreated in periods 0 and 1, then treated in periods 2, 3, and 4.
The packed table has one row per unit because five periods fit in one 64-bit word. With 365 periods, the same layout has six rows per unit.
Direct Bitmap Query
The static Mundlak bitmap query computes the seven scalar margins needed for \(X'X\) and \(X'Y\):
```sql
WITH unit_treatment AS (
SELECT
unit_index,
SUM(bit_count(w_bits & valid_mask))::DOUBLE AS c
FROM unit_bits
GROUP BY unit_index
),
unit_y AS (
SELECT
u.unit_index,
SUM(y.outcome_value::DOUBLE * bit_count(y.y_bits & u.valid_mask)) AS y_sum,
SUM(y.outcome_value::DOUBLE * bit_count(y.y_bits & u.w_bits & u.valid_mask)) AS wy_sum
FROM unit_bits AS u
JOIN unit_outcome_bits AS y USING (unit_index, chunk)
GROUP BY u.unit_index
),
unit_stats AS (
SELECT
t.unit_index,
t.c,
COALESCE(y.y_sum, 0.0) AS y_sum,
COALESCE(y.wy_sum, 0.0) AS wy_sum
FROM unit_treatment AS t
LEFT JOIN unit_y AS y USING (unit_index)
),
time_treatment AS (
SELECT
time_index,
SUM(bit_count(w_bits & valid_mask))::DOUBLE AS d
FROM time_bits
GROUP BY time_index
),
time_y AS (
SELECT
t.time_index,
SUM(y.outcome_value::DOUBLE * bit_count(y.y_bits & t.valid_mask)) AS y_sum
FROM time_bits AS t
JOIN time_outcome_bits AS y USING (time_index, chunk)
GROUP BY t.time_index
),
time_stats AS (
SELECT
t.time_index,
t.d,
COALESCE(y.y_sum, 0.0) AS y_sum
FROM time_treatment AS t
LEFT JOIN time_y AS y USING (time_index)
)
SELECT
(SELECT SUM(c) FROM unit_stats) AS sum_w,
(SELECT SUM(c * c) FROM unit_stats) AS sum_c2,
(SELECT SUM(y_sum) FROM unit_stats) AS sum_y,
(SELECT SUM(wy_sum) FROM unit_stats) AS sum_wy,
(SELECT SUM(c * y_sum) FROM unit_stats) AS sum_cy,
(SELECT SUM(d * d) FROM time_stats) AS sum_d2,
(SELECT SUM(d * y_sum) FROM time_stats) AS sum_dy
```
The query runs on four bitmap tables:
unit_bits: treatment words by unit and 64-period chunk;
time_bits: treatment words by time period and 64-unit chunk;
unit_outcome_bits: outcome bitplanes by unit, 64-period chunk, and nonzero outcome value;
time_outcome_bits: outcome bitplanes by time period, 64-unit chunk, and nonzero outcome value.
The query returns the scalar moments needed to solve the Mundlak normal equations.
Point-Estimate Runtime Benchmark
The point-estimate benchmark uses the one-shot adoption dimensions from notebooks/event_study.ipynb:
This produces a 15 million row long panel. The comparison reports three estimators for the same static model:
pyfixest on the full 15 million row long table, estimating Y_it ~ W_it | unit_id + time_id;
DBMundlak, which computes sufficient statistics by grouping the long table’s Mundlak design rows in DuckDB;
direct bitmap moments, which assumes the panel is available in packed unit-major and time-major bitmap layouts.
The main runtime bars exclude data-generation time and one-time layout creation. Each method has a different natural storage assumption, so the second panel reports those layout costs separately.
fig, axes = plt.subplots(1, 2, figsize=(13, 4.8), constrained_layout=True)runtime_colors = ["#4C78A8", "#59A14F", "#E15759"]axes[0].bar(runtime_df["method"], runtime_df["seconds"], color=runtime_colors)axes[0].set_yscale("log")axes[0].set_ylabel("seconds, log scale")axes[0].set_title("Estimator runtime on 15M rows")axes[0].tick_params(axis="x", rotation=20)for i, value inenumerate(runtime_df["seconds"]): axes[0].text(i, value *1.10, f"{value:.3f}s", ha="center", va="bottom")layout_colors = ["#9C755F", "#BAB0AC", "#F28E2B"]axes[1].bar(layout_df["stage"], layout_df["seconds"], color=layout_colors)axes[1].set_yscale("log")axes[1].set_ylabel("seconds, log scale")axes[1].set_title("One-time setup from generated long data")axes[1].tick_params(axis="x", rotation=20)for i, value inenumerate(layout_df["seconds"]): axes[1].text(i, value *1.10, f"{value:.3f}s", ha="center", va="bottom")for ax in axes: ax.grid(axis="y", which="both", alpha=0.25) ax.set_axisbelow(True)plt.show()
Show code
runtime_df
method
seconds
0
pyfixest full data
8.156105
1
DBMundlak groupby
0.150643
2
bitmap direct moments
0.174097
Show code
layout_df
stage
seconds
0
data generation
0.779396
1
DuckDB table creation
0.494419
2
long panel to bitmap layout
12.793199
Show code
size_df
quantity
value
0
long panel rows
15,000,000
1
long panel memory MB
120.0
2
unit-major treatment rows
500,000
3
time-major treatment rows
234,390
4
unit-major outcome bitplane rows
499,629
5
time-major outcome bitplane rows
234,390
Inference With Unit-Cluster Bootstrap
The point-estimate benchmark uses the packed panel once. Unit-cluster bootstrap inference reuses it many times. One bootstrap replicate samples exactly \(N\) units with replacement from the \(N\) observed units. The replicate is stored as a vector of unit multiplicities \(m_i\), where \(m_i\) is the number of times unit \(i\) appears in that full-sample resample.
For the static estimator, unit-additive quantities become weighted sums over units. Time margins are recomputed as
The event-study estimator applies the same multiplicities to cohort-time counts and outcome sums. The long table and the long-format groupby compression stay out of the bootstrap loop.
The example below uses a smaller staggered-adoption panel so the render can compare analytic cluster-robust inference from pyfixest, unit-cluster bootstrap inference from DBMundlak, and unit-cluster bootstrap inference from the bitmap estimator in one run.
pyfixest_w_estimate =float(pyfixest_clustered.coef().loc["W_it"])pyfixest_w_se =float(pyfixest_clustered.se().loc["W_it"])mundlak_w_estimate =float(mundlak_boot.point_estimate.reshape(-1)[1])mundlak_w_se =float(np.sqrt(np.diag(mundlak_boot.vcov))[1])bitmap_w_estimate =float(static_summary["point_estimate"].loc["W_it"])bitmap_w_se =float(static_summary["standard_error"].loc["W_it"])static_se_comparison = pd.DataFrame( {"method": ["pyfixest full FE","DBMundlak unit bootstrap","bitmap unit bootstrap", ],"inference": ["analytic CRV1 clustered by unit",f"B={num_bootstraps} full-sample unit bootstrap replicates",f"B={num_bootstraps} full-sample unit bootstrap replicates", ],"W_it estimate": [ pyfixest_w_estimate, mundlak_w_estimate, bitmap_w_estimate, ],"W_it se": [ pyfixest_w_se, mundlak_w_se, bitmap_w_se, ], }).round({"W_it estimate": 4, "W_it se": 4})static_se_comparison
method
inference
W_it estimate
W_it se
0
pyfixest full FE
analytic CRV1 clustered by unit
0.0693
0.0031
1
DBMundlak unit bootstrap
B=200 full-sample unit bootstrap replicates
0.0693
0.0030
2
bitmap unit bootstrap
B=200 full-sample unit bootstrap replicates
0.0693
0.0030
Show code
static_inference_runtime = pd.DataFrame( {"method": ["pyfixest full FE","DBMundlak unit bootstrap","bitmap unit bootstrap", ],"setup_seconds": [0.0, mundlak_setup_seconds, bitmap_setup_seconds, ],"fit_inference_seconds": [ pyfixest_inference_seconds, mundlak_inference_seconds, bitmap_inference_seconds, ], })static_inference_runtime["total_seconds"] = ( static_inference_runtime["setup_seconds"]+ static_inference_runtime["fit_inference_seconds"])static_inference_runtime.round(4)
method
setup_seconds
fit_inference_seconds
total_seconds
0
pyfixest full FE
0.0000
9.4659
9.4659
1
DBMundlak unit bootstrap
0.0234
1.1359
1.1594
2
bitmap unit bootstrap
0.2459
0.0889
0.3349
Show code
fig, ax = plt.subplots(figsize=(8.8, 4.8), constrained_layout=True)runtime_plot = static_inference_runtime.copy()colors = ["#4C78A8", "#59A14F", "#E15759"]ax.bar(runtime_plot["method"], runtime_plot["total_seconds"], color=colors)ax.set_yscale("log")ax.set_ylabel("seconds, log scale")ax.set_title("Complete static estimation plus inference")ax.tick_params(axis="x", rotation=18)for i, value inenumerate(runtime_plot["total_seconds"]): ax.text(i, value *1.10, f"{value:.3f}s", ha="center", va="bottom")ax.grid(axis="y", which="both", alpha=0.25)ax.set_axisbelow(True)plt.show()
Coverage Over 1000 Replications
The coverage experiment uses the one-shot adoption DGP with the sharkfin post-treatment effect path: \(N=10{,}000\), \(T=30\), \(5{,}000\) treated units, and first treatment in period 15. Each Monte Carlo replication uses analytic unit-cluster CRV1 inference for pyfixest and \(B=200\) full-sample unit bootstrap replicates for each compressed estimator.
The DGP has a dynamic post-treatment effect, so the static coefficient is a projection estimand. For each replication, the target below is the coefficient on \(W_{it}\) from the static Mundlak projection with \(Y_{it}\) replaced by the DGP probability \(\mathbb E[Y_{it} \mid W, \alpha_i, \lambda_t]\). The simulation is run directly in this notebook.
The timing columns are per Monte Carlo replication. For the compressed estimators, mean_total_seconds includes table creation or bitmap packing; mean_fit_inference_seconds excludes that layout step and includes all \(B=200\) bootstrap replicates. The final column divides the fit/inference time by \(B\).
The event-study example reports bitmap bootstrap intervals by cohort and event time. It uses the same packed panel as the static bootstrap example.
The point-estimate benchmark separates kernel time from layout construction. pyfixest solves the full-data fixed-effect problem on the long table. DBMundlak follows the current duckreg sufficient-statistic route: create Mundlak-generated regressors, group identical design rows, and solve weighted least squares. The bitmap route starts from packed unit-major and time-major words and solves the moments directly.
The setup panel is part of the result. When data arrive as a normal long panel, this prototype pays a Python-side packing cost before reaching the fast kernel. The intended deployment is a native packed representation for dense binary panels, so point-estimate, event-study, and bootstrap queries reuse the same layout.
