iop

# Monte Carlo: unit effects in the ordered probit

`fe_bias.R` compares four ways of handling unit heterogeneity in a panel
ordered probit when a covariate is correlated with the unit effect:

* `pooled` -- ordered probit ignoring the unit effect;
* `fe` -- unit dummies (`oprobit(..., fe = "unit")`);
* `re` -- a unit random intercept by adaptive Gauss--Hermite quadrature
  (`oprobit(..., re = "unit")`), which assumes the effect is uncorrelated with
  the covariates;
* `mundlak` -- the random intercept plus unit means of the covariates
  (`mundlak()` + `re = "unit"`), the correlated-random-effects device.

DGP: `y* = 0.8 x1 - 0.5 x2 + alpha_g + e`, `alpha_g ~ N(0, 1)`,
`x1 = N(0, 1) + 0.5 alpha_g`, four categories, G = 100 units, T in
{4, 8, 16, 32}, 100 replications per T. `fe_bias_results.csv` holds the mean,
percent bias, and RMSE of the `x1` coefficient (true value 0.8).

The qualitative pattern -- the pooled and random-intercept estimators are
biased by the omitted correlated effect, the dummies carry the
incidental-parameters bias that fades with T (Greene 2004), and the Mundlak
device is essentially unbiased at every T -- is why the documentation
recommends `mundlak()` (or `re` with Mundlak means) for short panels and
reserves `fe` for long ones.

## Audit (`fe_bias_verify.R`, log in `fe_bias_verify.log`)

* The pooled estimator's bias has a closed form under this DGP
  (coefficient 1.2/sqrt(1.8) = 0.894, +11.8 percent); the MC means are
  0.893-0.895 at every T.
* On identical draws, `oprobit(fe =)` equals `ordinal::clm()` with unit
  dummies to 1e-10 (same log-likelihood), `oprobit(re =)` equals
  `ordinal::clmm()` to 1e-6 (same sigma), and `mundlak()` + `re` equals
  `clmm()` with hand-computed unit means to 1e-6.
* The dummy bias is invariant to the number of units (T = 4: +28 percent at
  G = 100, +26 percent at G = 300) and a split-panel jackknife (Dhaene and
  Jochmans 2015) removes it (T = 8: +13 percent -> -3 percent): the
  incidental-parameters signature.
* Dropping units without identifying variation does not bias anything: the
  Mundlak estimator on the FE-retained sample (88 percent of observations at
  T = 4) is unbiased (-0.2 percent).

# Monte Carlo: size of `split_test()` under a common split

`split_test_size.R` draws 120 data sets from a common-split ZiOP (n = 1500,
J = 3, `beta = (0.8, -0.5)`, `tau = (-0.6, 0.7)`, `gamma = (0.3, 1)` on
`(1, z1)`, seeds `1000 + r`), fits the common-split model, and records the
Lagrange-multiplier statistic of `split_test()` with the outer-product (OPG,
the default) and observed-information (Hessian) forms and the likelihood-ratio
statistic from the category-specific refit (2 restrictions). It writes
`split_test_size_results.csv` (one row per replication) and prints the
rejection rates at the 10, 5, and 1 percent levels and the quantiles of the
LM/LR ratios; `split_test_size.log` is the output of the run reported in the
documentation. About six minutes.