Metadata-Version: 2.4
Name: mllpanel
Version: 0.1.0
Summary: Machine learning based panel data models with fixed effects and cross-sectional dependence (Yang, Long and Cai, 2024)
Author-email: Merwan Roudane <merwanroudane920@gmail.com>
Maintainer-email: Merwan Roudane <merwanroudane920@gmail.com>
License: MIT
Project-URL: Homepage, https://github.com/merwanroudane/mllpanel
Project-URL: Source, https://github.com/merwanroudane/mllpanel
Project-URL: Guide, https://github.com/merwanroudane/mllpanel/blob/main/docs/GUIDE.md
Project-URL: Issues, https://github.com/merwanroudane/mllpanel/issues
Keywords: panel data,machine learning,double machine learning,debiased machine learning,cross-sectional dependence,common correlated effects,fixed effects,total factor productivity,econometrics
Classifier: Development Status :: 4 - Beta
Classifier: Intended Audience :: Science/Research
Classifier: License :: OSI Approved :: MIT License
Classifier: Operating System :: OS Independent
Classifier: Programming Language :: Python :: 3
Classifier: Programming Language :: Python :: 3.9
Classifier: Programming Language :: Python :: 3.10
Classifier: Programming Language :: Python :: 3.11
Classifier: Programming Language :: Python :: 3.12
Classifier: Topic :: Scientific/Engineering :: Mathematics
Requires-Python: >=3.9
Description-Content-Type: text/markdown
License-File: LICENSE
Requires-Dist: numpy>=1.21
Requires-Dist: pandas>=1.3
Requires-Dist: scikit-learn>=1.1
Requires-Dist: scipy>=1.7
Provides-Extra: plots
Requires-Dist: matplotlib>=3.5; extra == "plots"
Provides-Extra: parallel
Requires-Dist: joblib>=1.1; extra == "parallel"
Provides-Extra: dev
Requires-Dist: matplotlib>=3.5; extra == "dev"
Requires-Dist: joblib>=1.1; extra == "dev"
Requires-Dist: pytest>=7.0; extra == "dev"
Requires-Dist: tabulate>=0.9; extra == "dev"
Requires-Dist: openpyxl>=3.0; extra == "dev"
Dynamic: license-file

# mllpanel

[![PyPI version](https://img.shields.io/pypi/v/mllpanel.svg)](https://pypi.org/project/mllpanel/)
[![Python versions](https://img.shields.io/pypi/pyversions/mllpanel.svg)](https://pypi.org/project/mllpanel/)
[![License: MIT](https://img.shields.io/badge/license-MIT-blue.svg)](https://github.com/merwanroudane/mllpanel/blob/main/LICENSE)
[![Tests](https://img.shields.io/badge/tests-148%20passing-brightgreen.svg)](https://github.com/merwanroudane/mllpanel/tree/main/tests)

**Machine learning based panel data models with fixed effects and
cross-sectional dependence.**

A complete, tested Python implementation of

> Yang, B., Long, W. and Cai, Z. (2024). *Machine Learning Based Panel Data
> Models*. Working paper.

`mllpanel` estimates a panel regression in which **the functional form is
unknown**, the **individual effects may be correlated with the regressors**, and
the errors carry an **unobserved factor structure** — using any machine learning
method for the unknown function, and **no proxy variable and no instrument** for
the latent terms.

---

## 📖 Start here: the researcher's guide

### **[→ docs/GUIDE.md — How to write a complete application](https://github.com/merwanroudane/mllpanel/blob/main/docs/GUIDE.md)**

An 18-section walkthrough from a raw panel to manuscript-ready tables and
figures: is the method right for your data, how to diagnose whether you
actually have these problems, how to choose a learner and a specification, how
to check condition (C3) before believing anything, which standard errors to
report, how to read a nonparametric production function, how to benchmark
against Olley–Pakes and Levinsohn–Petrin, how to forecast — plus **the complete
script**, a **12-row pitfalls table** and a **reporting checklist** for
referees.

Supporting documents:

| Document | What it covers |
|---|---|
| **[GUIDE.md](https://github.com/merwanroudane/mllpanel/blob/main/docs/GUIDE.md)** | **The researcher's guide — read this first** |
| [SYNTAX.md](https://github.com/merwanroudane/mllpanel/blob/main/docs/SYNTAX.md) | Complete argument-by-argument reference for every class and function |
| [THEORY.md](https://github.com/merwanroudane/mllpanel/blob/main/docs/THEORY.md) | The derivation equation by equation, a fidelity map from paper equations to code, every deviation with its justification, and an errata list for the published paper |
| [CHANGELOG.md](https://github.com/merwanroudane/mllpanel/blob/main/CHANGELOG.md) | Release notes |
| [examples/](https://github.com/merwanroudane/mllpanel/tree/main/examples) | Six runnable scripts |

---

## Contents

- [The problem this solves](#the-problem-this-solves)
- [The method](#the-method)
- [Install](#install)
- [Sixty seconds](#sixty-seconds)
- [API tour](#api-tour)
- [Real economic data, bundled](#real-economic-data-bundled)
- [Two results the paper does not report](#two-results-the-paper-does-not-report)
- [Reproducing the paper](#reproducing-the-paper)
- [Specifications](#specifications)
- [Learners](#learners)
- [Standard errors](#standard-errors)
- [Faithfulness to the paper](#faithfulness-to-the-paper)
- [Tables and figures](#tables-and-figures)
- [Cost and performance](#cost-and-performance)
- [Troubleshooting](#troubleshooting)
- [Testing](#testing)
- [Citation](#citation)
- [License](#license)

---

## The problem this solves

Suppose you want to estimate a production function, a demand curve, or an asset
pricing relation on panel data, and three things are true at once:

1. **You do not want to impose a functional form.** Cobb–Douglas is a
   restriction, not a finding.
2. **Unit-level heterogeneity is correlated with your regressors.** Firms with
   high unobserved productivity hire more labour, so OLS is biased.
3. **There are shocks that hit everyone at once.** A business cycle, an oil
   price shock, a policy change. These make the errors cross-sectionally
   dependent and — if they also move your regressors — endogenous.

The standard toolkit forces a choice. Nonparametric panel estimators (sieves,
kernels) handle (1) but collapse under the curse of dimensionality beyond one or
two regressors. Machine learning panel methods handle (1) and scale, but
existing ones either pool the data outright or put fixed effects in as
predictors, leaving (3) untouched. Control-function estimators
(Olley–Pakes, Levinsohn–Petrin) handle (2) but require a proxy variable, impose
a parametric production function, and say nothing about (3).

**This method handles all three simultaneously, and needs no proxy.** The
insight — from Mundlak (1978), Pesaran (2006) and Huang (2013) — is that both
latent terms are *spanned by observable averages*, so they can be replaced by
regressors you can compute. What is left is a partially linear model, which the
debiased machine learning of Chernozhukov et al. (2018) estimates with any
learner you like.

**When not to use it:** if you are happy with a linear form, a two-way fixed
effects regression is faster and has cleaner inference. If your application
turns on *heterogeneity* in the factor loadings, this method implicitly assumes
that away — see
[THEORY.md §2](https://github.com/merwanroudane/mllpanel/blob/main/docs/THEORY.md#2-the-transformation-appendix-a).
And it needs both `N` and `T` reasonably large; see the sample-size table in
[GUIDE.md §1](https://github.com/merwanroudane/mllpanel/blob/main/docs/GUIDE.md#1-is-this-the-right-method-for-your-data).

---

## The method

The model — the paper's equation (1):

$$
y_{it} = f(x_{it}) + \alpha_i + \gamma_i'\lambda_t + \varepsilon_{it}
$$

with $f$ unknown, the fixed effects $\alpha_i$ freely correlated with $x_{it}$,
and the common factors $\lambda_t$ also allowed to correlate with $x_{it}$.

**The trick** — the paper's equation (2), derived in its Appendix A by a Taylor
expansion inside a shrinking neighbourhood of $x$. Both latent terms are spanned
by observable averages: the fixed effects by **time averages**, the common
factors by **cross-sectional averages**. So the model becomes partially linear:

$$
y_{it} = f(x_{it}) + \beta' z_{it} + e_{it},
\qquad
z_{it} = \bigl(\underbrace{\bar y_t - \bar y,\ (\bar x_t - \bar x)'}_{\text{spans } \lambda_t},\
               \underbrace{\bar y_i - \bar y,\ (\bar x_i - \bar x)'}_{\text{spans } \alpha_i}\bigr)'
$$

**The estimation** — equations (3)–(6):

1. Partial out $g = E(y\mid x)$ and $m = E(z\mid x)$ by machine learning.
2. **Cross-fit over the time index.** Folds are over *periods, not
   observations*: every unit in period $t$ shares a fold, which is what stops
   $\lambda_t$ leaking across the split.
3. Solve the Neyman-orthogonal moment condition in closed form.
4. Recover $f$ from a second-stage regression of $y - \hat\beta'z$ on $x$.

**Why the orthogonality matters:** it makes $\hat\beta$ robust to first-order
error in the nuisance estimates. That is what lets you use a random forest for
$f$ without regularisation bias contaminating everything else.

---

## Install

From PyPI:

```bash
pip install mllpanel
```

With figures and parallel Monte Carlo:

```bash
pip install "mllpanel[plots,parallel]"
```

From source, for development:

```bash
git clone https://github.com/merwanroudane/mllpanel.git
cd mllpanel
pip install -e ".[dev]"
pytest -q
```

| Extra | Adds | For |
|---|---|---|
| *(none)* | `numpy`, `pandas`, `scikit-learn`, `scipy` | estimation, tables |
| `plots` | `matplotlib` | all figures |
| `parallel` | `joblib` | `n_jobs=` in Monte Carlo and cross-fitting |
| `dev` | the above + `pytest`, `tabulate`, `openpyxl` | tests, Markdown tables, xlsx export |

Requires Python ≥ 3.9. The four real datasets ship inside the wheel, so
**every example in this documentation runs offline**.

---

## Sixty seconds

```python
from mllpanel import MLLPanel
from mllpanel.datasets import load_snmesp

df = load_snmesp()          # 738 Spanish manufacturing firms x 8 years, real data

model = MLLPanel(
    learner="rf",           # random forest for the unknown production function
    spec="crfe",            # cross-sectional AND time averages
    vcov="twoway",          # two-way clustered standard errors
    random_state=0,
).fit(df, y="y", x=["n", "k"], unit="firm", time="year")

print(model.summary())
```

```
====================================================================================
                      Machine Learning Based Panel Data Model
               Yang, Long and Cai (2024) - debiased machine learning
====================================================================================
Dep. variable     y                       No. observations                    5904
Specification     CR + FE                 No. units (N)                        738
Learner (g, m, f) Random Forests          No. periods (T)                        8
Cross-fitting     K=5, block              Panel                           balanced
beta aggregation  pooled                  Covariates (d)                         2
                                          Augmented regressors                   6
Std. errors       Two-way clustered (unit and period)
------------------------------------------------------------------------------------
                    Parametric component of  y = f(x) + b'z + e
------------------------------------------------------------------------------------
Regressor                    Coef.   Std. err.           z       P>|z|  [2.5%  97.5%]
  Cross-sectional averages (common factors)
ybar_t                      1.0020      0.0812     12.3390      0.0000  0.8429  1.1612
xbar_t[n]                  -0.4200      0.4323     -0.9716      0.3313 -1.2673  0.4273
xbar_t[k]                  -0.4202      0.1017     -4.1320      0.0000 -0.6195 -0.2209
  Time averages (fixed effects)
ybar_i                      0.9729      0.0046    210.1484      0.0000  0.9639  0.9820
xbar_i[n]                  -0.6753      0.0336    -20.1247      0.0000 -0.7411 -0.6096
xbar_i[k]                  -0.2874      0.0291     -9.8839      0.0000 -0.3444 -0.2304
------------------------------------------------------------------------------------
Specification tests
  H0: cross-sectional averages jointly zero (no unobserved common factors)
    chi2(3) = 566.9690    p = 0.0000
  H0: time averages jointly zero (individual effects uncorrelated with x; Mundlak)
    chi2(3) = 52015.0492    p = 0.0000
------------------------------------------------------------------------------------
Fit:  R2 = 0.9976   R2 (uncentred, GKX) = 0.9999   RMSE = 0.0677
Nuisance out-of-fold R2:  g = 0.8588   m: min -0.4677, median 0.2330, max 0.9900
====================================================================================
```

**A free correctness check is built into that output.** Appendix A of the paper
implies a coefficient of exactly **one** on both `ybar_t` and `ybar_i` when
loadings are homogeneous, and minus the input elasticity on `xbar_i[·]`. The
estimates are 1.002 and 0.973, with −0.675 and −0.287 against elasticities of
+0.63 and +0.33. If your own fit does not show this pattern, stop and
investigate before interpreting anything.

Then read the economics off the nonparametric part:

```python
ape = model.partial_effects()          # output elasticities, since y and x are logs
print(ape[["APE", "Median"]].round(3))
#             APE  Median
# n         0.627   0.132
# k         0.334   0.226
print("Returns to scale:", round(ape["APE"].sum(), 3))   # 0.961

prod = model.productivity()            # firm-year TFP decomposition
```

---

## API tour

Fifteen modules. Full signatures in
[SYNTAX.md](https://github.com/merwanroudane/mllpanel/blob/main/docs/SYNTAX.md).

### The estimator

```python
from mllpanel import MLLPanel, fit_mllpanel

MLLPanel(
    learner="lasso", spec="crfe", *,
    learner_f=None, learner_m=None,          # different learners per stage
    n_folds=5, fold_scheme="block",          # cross-fitting
    drop_remainder=False, beta_agg="pooled",
    vcov="twoway", vcov_lags=None,
    f_crossfit=False, grand_mean="obs",
    ybar_forecast_learner="ols",
    alpha=0.05, random_state=None,
    store_models=False, n_jobs=None, verbose=False,
)
```

Three interchangeable call styles:

```python
model.fit(df, y="gsp", x=["emp", "pc"], unit="state", time="year")   # DataFrame
model.fit(panel_data_object)                                        # PanelData
model.fit(y=y_arr, x=X_arr, unit=id_arr, time=t_arr)                # arrays
```

Key fitted attributes:

| Attribute | Meaning |
|---|---|
| `beta_`, `vcov_`, `inference_` | $\hat\beta$, its covariance, the full coefficient table |
| `wald_` | block tests: `"common_factors"`, `"fixed_effects"` |
| `f_`, `f_hat_` | the fitted second-stage model, and $\hat f(x_{it})$ |
| `fittedvalues_`, `resid_` | $\hat y$ and $y-\hat y$ |
| `nuisance_fit_` | **out-of-fold $R^2$ per nuisance regression — the practical check on condition (C3)** |
| `z_`, `z_tilde_`, `y_tilde_` | augmented regressors and their cross-fitted residualised versions |
| `r2_`, `r2_uncentered_`, `rmse_` | centred $R^2$ (report this), the paper's uncentred one, RMSE |
| `folds_`, `used_`, `n_dropped_` | the time-index partition and what entered the score |

Methods:

```python
model.predict(df_new, unit="firm", time="year",
              ybar_t_source="auto",     # "model" for honest forecasting
              unseen_units="error",
              return_components=False)
model.predict_f(X)                      # f-hat alone; levels not identified
model.partial_effects(relative_step=0.05)     # elasticities, if y and x are logs
model.partial_dependence("k", n_grid=60, at="average")
model.productivity()                    # f_hat, factor_part, fe_part, tfp_hat
model.summary(); model.summary_frame(); model.info(); model.to_latex()
```

### Diagnostics

```python
from mllpanel.diagnostics import cd_test, dependence_report, residual_diagnostics

cd_test(model.resid_, unit, time)       # Pesaran CD + mean |rho|
dependence_report(df, y, x, unit, time, specs=("pooled","fe","cr","crfe"))
```

### Simulation and Monte Carlo

```python
from mllpanel import simulate_panel, time_split, run_design, run_sensitivity
from mllpanel.montecarlo import (paper_designs, inference_study,
                                 theoretical_beta, estimate_runtime, aggregate)

df = simulate_panel(n_units=10, n_times=100, d=5, c1=1.0, c2=1.0,
                    ftype="type1", gamma=0.5, factor_ar=0.0, seed=0)
parts = time_split(df, fractions=(0.3, 0.2, 0.5))       # consecutive periods

results = run_design(paper_designs(), learners=("lasso","rf","nn"),
                     n_reps=1000, n_jobs=8)
agg = aggregate(results)                 # + Monte Carlo standard errors
```

### Production functions

```python
from mllpanel.tfp import (estimate_pooled, estimate_fe, estimate_op,
                          estimate_lp, estimate_acf,
                          specification_table, comparison_table)

fits = [estimate_pooled(df, "y", ["n","k"], "firm", "year", cluster="firm"),
        estimate_fe(df, "y", ["n","k"], "firm", "year", effects="twoway"),
        estimate_op(df, "y", "n", "k", "log_inv", "firm", "year"),
        estimate_lp(df, "y", "n", "k", "i", "firm", "year"),
        estimate_acf(df, "y", "n", "k", "i", "firm", "year")]
specification_table(fits)                                       # Table 7
comparison_table(fits, df, {"ln M": "i", "ln I": "log_inv"})     # Table 8
```

### Output

```python
from mllpanel.tables import save_table, montecarlo_table, table_to_latex
from mllpanel.plots import (set_style, save_figure, plot_sensitivity,
                            plot_coefficients, plot_partial_dependence,
                            plot_surface, plot_factor_path, plot_tfp,
                            plot_nuisance_quality, plot_dependence_report)

set_style()                                        # 9 pt serif, journal defaults
save_table(tab, "out/table_1", formats=("csv","tex","md"),
           caption="...", label="tab:1", notes="...", highlight_min=True)
save_figure(plot_surface(model, "n", "k"), "out/fig_3", formats=("pdf","png"))
```

### Meta

```python
import mllpanel
mllpanel.__version__
print(mllpanel.deviations())          # every implementation choice, at runtime
print(mllpanel.citation("bibtex"))
```

### Module map

| Module | Contents |
|---|---|
| `mllpanel.core` | `MLLPanel` — the estimator |
| `mllpanel.panel` | `PanelData`, the six averages, `build_z` |
| `mllpanel.crossfit` | sample splitting over the time index, three schemes |
| `mllpanel.learners` | LASSO, random forests, 32-16-8 network, plus five more |
| `mllpanel.tune` | validation-block hyper-parameter selection (the paper's protocol) |
| `mllpanel.inference` | five sandwich estimators, Wald tests on each block |
| `mllpanel.metrics` | $MSE(f)$, in- and out-of-sample $R^2$, centred and uncentred |
| `mllpanel.simulate` | the DGP of equation (7), Type 1 and Type 2, the 30/20/50 split |
| `mllpanel.montecarlo` | Tables 1–6, Figures 1–6, inference study, runtime estimation |
| `mllpanel.tfp` | pooled, two-way FE, Olley–Pakes, Levinsohn–Petrin, ACF, Tables 7–8 |
| `mllpanel.diagnostics` | Pesaran CD test, serial correlation, `dependence_report` |
| `mllpanel.tables` | `booktabs` LaTeX, Markdown, CSV, xlsx |
| `mllpanel.plots` | ten journal-quality figure types |
| `mllpanel.colors` | MATLAB `parula`, Okabe–Ito colour-blind-safe palette |
| `mllpanel.datasets` | four real economic panels |

### Examples

| Script | What it does |
|---|---|
| [`01_quickstart.py`](https://github.com/merwanroudane/mllpanel/blob/main/examples/01_quickstart.py) | the minimum viable application |
| [`02_real_data_tfp.py`](https://github.com/merwanroudane/mllpanel/blob/main/examples/02_real_data_tfp.py) | the full Section 4 exercise: 3 learners, 6 benchmarks, Tables 7–8, 8 figures. **Copy this one.** |
| [`03_replicate_tables.py`](https://github.com/merwanroudane/mllpanel/blob/main/examples/03_replicate_tables.py) | Monte Carlo Tables 1–6 |
| [`04_replicate_figures.py`](https://github.com/merwanroudane/mllpanel/blob/main/examples/04_replicate_figures.py) | Monte Carlo Figures 1–6 |
| [`05_forecasting.py`](https://github.com/merwanroudane/mllpanel/blob/main/examples/05_forecasting.py) | out-of-sample and rolling one-step-ahead forecasts |
| [`06_inference_study.py`](https://github.com/merwanroudane/mllpanel/blob/main/examples/06_inference_study.py) | coverage of the Theorem 1 sandwich |

---

## Real economic data, bundled

The paper's own application uses licensed CSMAR data on Chinese A-share firms,
which cannot be redistributed. Four public panels ship inside the wheel.

| Loader | $N$ | $T$ | Balanced | Why it is here |
|---|---|---|---|---|
| `load_snmesp()` | 738 | 8 | yes | **Spanish manufacturing, 1983–1990.** Log output, labour, capital **and log intermediate consumption**, so it supports the full Table 8 comparison against OP, LP and ACF. The closest public analogue of the paper's Section 4, which used 459 Chinese firms over 20 years. |
| `load_produc()` | 48 | 17 | yes | **US states, 1970–1986.** Real gross state product, labour, private and public capital. Small $N$, long $T$. |
| `load_grunfeld()` | 10 | 20 | yes | **Firm investment, 1935–1954.** $N=10$ makes the $O(1/N)$ simultaneity of the CCE step visible. A cautionary example. |
| `load_empluk()` | 140 | 7–9 | **no** | **UK firms, 1976–1984.** For checking that your code path handles unequal $T_i$. |

Sources: Alonso-Borrego and Arellano (1999); Munnell (1990); Grunfeld (1958);
Arellano and Bond (1991) — all redistributed from the R package `plm` (GPL-2).

> `snmesp` has no investment series, so `log_inv` is imputed from
> $I_{it}=K_{it}-(1-\delta)K_{i,t-1}$. This is non-positive for **15.5%** of the
> panel — not a defect, but exactly the censoring the paper cites as the
> motivation for Levinsohn–Petrin over Olley–Pakes. Because it is imputed rather
> than measured, its column in a Table 8 replication is not numerically
> comparable to the paper's.

---

## Two results the paper does not report

### 1. The cross-sectional averages really do absorb the dependence

The paper's central claim is directly testable, and the paper does not test it.
`diagnostics.dependence_report` does. On the Spanish panel:

| Specification | Pesaran CD | p-value | mean ρ | mean \|ρ\| |
|---|---|---|---|---|
| Pooled | 276.60 | 0.000 | 0.188 | 0.483 |
| FE | 174.80 | 0.000 | 0.119 | 0.457 |
| CR | −1.05 | 0.293 | −0.001 | 0.429 |
| **CR + FE** | **−1.14** | **0.253** | **−0.001** | 0.430 |

The CD statistic collapses from **277 to −1.1**. Adding *additive* fixed effects
barely helps (175) — a factor structure is not an additive effect, which is the
paper's Panel C result appearing in real data.

**But read the last column.** Mean *absolute* pairwise correlation moves only
from 0.483 to 0.430. The CD test has no power against a factor structure whose
loadings average to zero, so together these columns say the block removes the
*average* co-movement while individual firm pairs remain strongly correlated.
`mllpanel` reports both, always.

![Cross-sectional dependence by specification](https://raw.githubusercontent.com/merwanroudane/mllpanel/main/docs/figures/fig_A1_dependence.png)

### 2. Which standard errors actually work

Theorem 1 gives $\Sigma^{-1}V\Sigma^{-1}$ but supplies **no sample analogue**,
and the paper reports no standard errors anywhere and runs no coverage
experiment. `examples/06_inference_study.py` runs one — Type 1 design, $N=20$,
$T=200$, $d=5$, 150 replications, against the closed-form $\beta$ from
Appendix A:

| `vcov` | CR block se_ratio | CR coverage | FE block se_ratio | FE coverage |
|---|---|---|---|---|
| `theorem1` (i.i.d. plug-in) | 6.63 | 1.000 | 16.37 | 1.000 |
| `cluster_unit` | 6.75 | 1.000 | **0.86** | **0.892** |
| `cluster_time` | **0.92** | **0.912** | 16.87 | 1.000 |
| `dk` (Driscoll–Kraay) | **0.92** | **0.908** | 16.66 | 1.000 |
| `twoway` (default) | 3.99 | 0.989 | 6.53 | 1.000 |

`se_ratio` = reported standard error ÷ true Monte Carlo standard deviation; one
is the target.

**Match the clustering dimension to the block you are testing.** The
cross-sectional-average block varies over $t$, so cluster by period; the
time-average block varies over $i$, so cluster by unit. The **literal i.i.d.
plug-in is 6–16× too wide**, so a test built on it would essentially never
reject — worth knowing before anyone implements Theorem 1 naively.

The sandwich is validated against `statsmodels`: `theorem1` and `cluster_unit`
reproduce `OLS(...).fit(cov_type="HC0")` and `cov_type="cluster"` on the
residualised variables to six decimal places.

---

## Reproducing the paper

```bash
# Tables 1-6.  Start small; --estimate-runtime before anything ambitious.
python examples/03_replicate_tables.py --reps 20 --learners lasso rf nn
python examples/03_replicate_tables.py --reps 1000 --full --jobs 16   # days

# Figures 1-6
python examples/04_replicate_figures.py --reps 200 --band --jobs 8

# Section 4 on real data
python examples/02_real_data_tfp.py --outdir output
```

| Paper | Command |
|---|---|
| Tables 1–2 ($MSE(f)$) | `03_replicate_tables.py`, `value="mse_f"` |
| Tables 3–4 (in-sample $R^2$) | `value="r2_is"` |
| Tables 5–6 (out-of-sample $R^2$) | `value="r2_oos"` |
| Figures 1–6 | `04_replicate_figures.py` |
| Tables 7–8 | `02_real_data_tfp.py` |

### Does it replicate?

$MSE(f)$ under Type 1 at $N=10$, $T=100$, $d=5$, LASSO — this package at **10**
replications against the paper's 1,000:

| Panel | | Pooled | FE | CR + FE |
|---|---|---|---|---|
| A ($c_1=c_2=0$) | paper | 0.0109 | 0.0116 | 0.0146 |
| | `mllpanel` | 0.0192 | 0.0181 | 0.0176 |
| B ($c_1=1, c_2=0$) | paper | 0.2269 | 0.0111 | 0.0145 |
| | `mllpanel` | 0.2634 | 0.0241 | 0.0307 |
| C ($c_1=c_2=1$) | paper | 0.5601 | 0.1182 | 0.0164 |
| | `mllpanel` | 0.4138 | 0.1302 | 0.0227 |

The pattern that carries the paper's argument reproduces cleanly: pooled
estimation collapses the moment either latent term appears, FE alone is not
enough once there are common factors, and CR + FE is an order of magnitude
better in Panel C. All three panels of Table 5 (out-of-sample $R^2$) reproduce
the paper's ordering exactly.

**Where it diverges, and why.** The random forest columns come out roughly four
times larger than the paper's. The paper does not report its forest
hyper-parameters, and scikit-learn's defaults grow trees fully, which overfits
300 training observations hard. Tuning on the validation block (`--tune`) closes
much of the gap. We keep the transparent default rather than reverse-engineering
values that happen to match the published numbers — see
[THEORY.md §6.1](https://github.com/merwanroudane/mllpanel/blob/main/docs/THEORY.md#61-choices-the-paper-leaves-open).

### Section 4 on real data

Returns to scale across all nine specifications on the Spanish panel:

| Model | Returns to scale | Table 8: `ln M` | `ln I` |
|---|---|---|---|
| LASSO | 0.938 | 0.145 | 0.030 |
| Random Forests | 0.994 | 0.147 | 0.018 |
| Neural Networks | 0.905 | 0.190 | 0.080 |
| Pooled | 0.951 | 0.475 | 0.069 |
| Fixed Effects | 0.936 | 0.139 | 0.023 |
| Olley and Pakes (1996) | 1.632 | 0.056 | — |
| Levinsohn and Petrin (2003) | 0.481 | — | 0.042 |
| ACF, $h(\ln M, k, n)$ | 1.092 | — | 0.048 |
| ACF, $h(\ln I, k, n)$ | 0.958 | 0.057 | — |

Reported as found, not as one might wish. Three things to note:

- The three machine learning specifications and the two simple benchmarks all
  land near constant returns; **OP (1.63) and LP (0.48) do not**. LP's labour
  elasticity is 0.20 against 0.62–0.67 everywhere else — the
  Gandhi–Navarro–Rivers (2020) critique biting, because this is a *gross output*
  measure and materials is a flexible input. A genuine property of applying LP
  here, not a bug.
- On `ln M`, the control-function estimators do **better** (0.056, 0.057) than
  the machine learning rows (0.145–0.190) — the reverse of what the paper found
  on Chinese data.
- The machine learning rows carry $2d+2$ extra regressors, so a smaller
  correlation is partly mechanical. Table 8 is not a like-for-like contest.

![Isoquants of the estimated production function](https://raw.githubusercontent.com/merwanroudane/mllpanel/main/docs/figures/fig_3_isoquants.png)

Parallel straight contours would mean Cobb–Douglas was adequate. This is the
single most informative picture of what the machine learning step bought.

---

## Specifications

| `spec` | $\dim z$ | Blocks | Paper's label |
|---|---|---|---|
| `"crfe"` | $2d+2$ | cross-sectional **and** time averages | CR + FE |
| `"fe"` | $d+1$ | time averages only | FE |
| `"cr"` | $d+1$ | cross-sectional averages only | *addition* |
| `"pooled"` | 0 | none | Pooled |

`"cr"` isolates the contribution of the CCE block, which the paper's three
specifications cannot. **Report more than one** — the comparison *is* the
evidence.

---

## Learners

The estimator is learner-agnostic: anything satisfying condition (C3),
$\hat\eta-\eta = o_p((NT)^{-1/4})$, will do.

| Key | Estimator | In the paper |
|---|---|---|
| `"lasso"` | standardised `LassoCV` | yes |
| `"lasso_fixed"` | `Lasso` with an explicit penalty, for validation-block tuning | — |
| `"rf"` | `RandomForestRegressor(n_estimators=500)` | yes |
| `"nn"` | `MLPRegressor(hidden_layer_sizes=(32,16,8))` — the Gu–Kelly–Xiu architecture | yes |
| `"ols"` | `LinearRegression` — makes the estimator analytically checkable | — |
| `"ridge"`, `"gbr"`, `"et"`, `"krr"` | ridge, boosting, extra trees, kernel ridge | boosting named as admissible |

Any scikit-learn regressor works, and the three stages can differ:

```python
MLLPanel(learner="rf",       # flexible g and m — these only need to predict well
         learner_f="krr")    # smooth f — this one gets differentiated
```

The paper's own simulation result is a clean rule: **LASSO when $f$ is sparse
and linear** (Tables 1, 5), **forests or networks when it is not** (Tables 2, 6,
where LASSO hits a bias floor no amount of data fixes).

> `"lasso"` is a `LassoCV` and has no `alpha` to tune. For the paper's
> validation-block protocol use `"lasso_fixed"` with
> `grid={"model__alpha": ...}`.

---

## Standard errors

| `vcov` | Meat matrix | Use when |
|---|---|---|
| `"theorem1"` / `"hc"` | $\sum\psi\psi'$ | you want the paper's literal formula (but see above) |
| `"cluster_unit"` | by unit | the time-average block carries your argument |
| `"cluster_time"` | by period | the cross-sectional-average block carries your argument |
| `"twoway"` | Cameron–Gelbach–Miller, PSD-projected | **default**; never anti-conservative |
| `"dk"` | Driscoll–Kraay with Bartlett HAC | as `cluster_time`, with serial correlation |

Two Wald tests run automatically:

```python
model.wald_["common_factors"]   # H0: no unobserved common factors
model.wald_["fixed_effects"]    # H0: Mundlak — effects uncorrelated with x
```

> When a block's switch is genuinely off, its coefficients are individually
> **unidentified** (the block is collinear up to an $O_p(T^{-1/2})$ term) even
> though its fitted *contribution* is pinned at zero. Test the contribution via
> `productivity()["factor_part"]`, not a single coefficient.

---

## Faithfulness to the paper

Every equation is mapped to the code that implements it in
[THEORY.md §5](https://github.com/merwanroudane/mllpanel/blob/main/docs/THEORY.md#5-fidelity-map-equation--code).
Where the paper leaves a choice open, `mllpanel` picks a documented default and
exposes the alternative:

| Choice | Paper | Default here | Alternative |
|---|---|---|---|
| Fold construction | unspecified | `fold_scheme="block"` | `"interleaved"`, `"random"` |
| Leftover periods when $T \nmid K$ | $T_0=[T/K]$ drops them | keep them | `drop_remainder=True` |
| $\hat\beta$ aggregation | eq. (6), pooled | `beta_agg="pooled"` | `"mean"` (footnote 1) |
| Second-stage $f$ | full sample | full sample | `f_crossfit=True` |
| Standard errors | none given | `vcov="twoway"` | `"theorem1"` and three others |
| $R^2$ | uncentred (Gu–Kelly–Xiu) | both reported; `r2_` is centred | `r2_uncentered_` |
| Forest hyper-parameters | unreported | 500 trees, sklearn defaults | tune on a validation block |

Two limitations of the method itself, documented because they bound what any
implementation can deliver:

- **$\beta$ is derived locally but used globally.** Appendix A makes the
  coefficients on $\bar x_t$ and $\bar x_i$ equal $-\theta(x)$, the gradient of
  $f$ *at the expansion point*. Under a nonlinear $f$ that varies with $x$, so a
  constant $\beta$ is an approximation whose error goes into $e_{it}$.
- **Homogeneous factor loadings are implicitly required.** The coefficient on
  $\bar y_t$ is $\gamma_i'(\bar\gamma\bar\gamma')^{-1}\bar\gamma$, which is
  unit-specific; a pooled $\beta_1$ needs it constant. The paper's simulations
  set $\gamma_i=0.5$ for every $i$, so the design never stresses this. Test it
  with `simulate_panel(gamma="random")`.

Also documented: a short
[errata list](https://github.com/merwanroudane/mllpanel/blob/main/docs/THEORY.md#66-errata-in-the-published-paper)
for the published paper — most visibly, every panel of Figures 1–6 is labelled
`"ISS R-Square"` on the vertical axis, including Figures 1 and 2, which plot
mean squared error. This package labels axes by the quantity plotted.

Run `print(mllpanel.deviations())` for the list at runtime.

---

## Tables and figures

Journal-quality output, ready to `\input` into a manuscript.

- Tables use `booktabs`: no vertical rules, `\multicolumn` groups with
  `\cmidrule`, full-width panel captions for Panels A/B/C, optional bolding of
  the best cell per row, and a notes block. LaTeX escaping is automatic but
  leaves any cell containing `$` alone, so `$d$` survives while `omega_it` gets
  its underscore escaped.
- Figures export as vector PDF with `pdf.fonttype 42`, so text stays editable,
  plus 400 dpi PNG for slides.
- Every figure uses colour **and** dash pattern **and** marker, so it survives
  greyscale printing.

Ten figure types: sensitivity panels (Figures 1–6), convergence, coefficient
plots, partial dependence, isoquant contours, 3-D surfaces, factor paths,
productivity distributions, nuisance-quality bars, dependence reports.

LaTeX preamble: `\usepackage{booktabs}`, plus `rotating` if you use
`landscape=True`.

---

## Cost and performance

Each fit trains $K(1+\dim z)$ nuisance models plus one second stage — with
`spec="crfe"`, $d=5$ and $K=5$ that is **61 model fits per estimation**. Plan
accordingly.

| Task | Rough cost |
|---|---|
| One `crfe` fit, LASSO, $N{=}738$, $T{=}8$, $d{=}2$ | ~1 s |
| One `crfe` fit, random forest, same | ~20 s |
| Section 4 script, 3 learners + 6 benchmarks + 8 figures | ~10 min |
| Tables 1–6, 6 designs × 10 reps × 2 learners | ~15 min |
| The published grid, 72 designs × 1000 reps × 3 learners × 3 specs | **days** |

Levers: `n_jobs=` parallelises folds and Monte Carlo replications;
`paper_designs(quick=True)` while developing; `montecarlo.estimate_runtime()`
times one replication per design and extrapolates before you commit;
`make_learner("rf", n_estimators=100)` for a cheaper forest.

---

## Troubleshooting

| Symptom | Cause | Fix |
|---|---|---|
| `RuntimeError: No candidate ... could be fitted` | `"lasso"` is a `LassoCV` with no `alpha` | use `"lasso_fixed"` |
| $R^2$ of 0.9999 for every model | uncentred $R^2$ on log levels | report `model.r2_`, not `r2_uncentered_` |
| Coefficient on `ybar_i` far from 1 | the panel does not match the model's structure, or $T$ is tiny | see [GUIDE.md §8](https://github.com/merwanroudane/mllpanel/blob/main/docs/GUIDE.md#8-fit-and-read-the-output-properly) |
| Negative out-of-fold $R^2$ in `nuisance_fit_` | a $z$ column varying over few periods, extrapolated by block folds | try `fold_scheme="interleaved"`; see [GUIDE.md §9](https://github.com/merwanroudane/mllpanel/blob/main/docs/GUIDE.md#9-check-condition-c3-before-you-believe-anything) |
| Wald test rejects but `factor_part` ≈ 0 | block coefficients unidentified when the switch is off | test the contribution, not a coefficient |
| `Duplicate (unit, time) pairs detected` | repeated cells | aggregate or de-duplicate first |
| Huge standard errors under `vcov="theorem1"` | the i.i.d. plug-in is 6–16× too wide | use `"twoway"`, or match the dimension |
| `singular ... using the pseudo-inverse` | a time-invariant regressor is collinear with its own time average | drop it from `x`, or use `spec="cr"` |
| Forest elasticities with enormous `Std. dev.` | piecewise-constant $\hat f$ | raise `relative_step`, or set `learner_f="krr"` |
| `predict` raises on unseen units | their time averages are unidentified | `unseen_units="grand_mean"`, and say so |
| `unbalanced` warning | Appendix A assumes balance | proceed but report it |

The full 12-row version, with explanations, is in
[GUIDE.md §16](https://github.com/merwanroudane/mllpanel/blob/main/docs/GUIDE.md#16-common-pitfalls).

---

## Testing

```bash
pytest -q      # 148 tests
```

The suite is not a smoke test. Under the paper's Type 1 design the whole
estimator has a closed-form answer, so
[`tests/test_oracle.py`](https://github.com/merwanroudane/mllpanel/blob/main/tests/test_oracle.py)
checks $\hat\beta$ against the algebra of Appendix A:

- $\hat\beta$ matches $(c_2, -c_2\theta', c_1, -c_1\theta')$ to within 0.03;
- the coefficients on `ybar_t` and `ybar_i` are 1, and on `xbar_i[x1..x3]` are
  $-0.2$, with the irrelevant regressors at 0;
- $MSE(f)$ orders `crfe` < `fe` < `pooled` under Panel C, with a factor of ten
  between the extremes, reproducing the paper's headline result;
- $MSE(f)$ halves as $T$ quadruples — Theorem 2 in practice;
- with no latent terms, `pooled` wins slightly, reproducing Panel A.

Also tested: every covariance estimator gives finite positive standard errors
and a PSD matrix; the score is orthogonal at $\hat\beta$; folds partition the
time index; the CD test detects a factor structure and accepts independence; all
ten figure types build; every LaTeX table renders.

---

## Citation

```bibtex
@unpublished{yang2024mlpanel,
  author = {Yang, Bingduo and Long, Wei and Cai, Zongwu},
  title  = {Machine Learning Based Panel Data Models},
  year   = {2024},
  note   = {Working paper}
}

@software{roudane2026mllpanel,
  author  = {Roudane, Merwan},
  title   = {mllpanel: Machine Learning Based Panel Data Models},
  version = {0.1.0},
  year    = {2026},
  url     = {https://github.com/merwanroudane/mllpanel}
}
```

Or `print(mllpanel.citation("bibtex"))`.

---

## License

MIT — see
[LICENSE](https://github.com/merwanroudane/mllpanel/blob/main/LICENSE).

The bundled datasets are redistributed from the R package `plm` under GPL-2;
their original sources are documented in `src/mllpanel/datasets/_loaders.py`.

---

## Author

**Dr Merwan Roudane** ·
[github.com/merwanroudane](https://github.com/merwanroudane)
