diff --git a/linearmodels/panel/model.py b/linearmodels/panel/model.py index a63ca58938..a377898bc8 100644 --- a/linearmodels/panel/model.py +++ b/linearmodels/panel/model.py @@ -3070,6 +3070,35 @@ def fit( ``debiased=True`` and ``group_debias=True`` give the clustered covariance of Stata's ``xtreg, re vce(cluster ...)``. See :ref:`panel-small-sample`. + + The variance of the idiosyncratic error is the residual variance of the + within regression, and the variance of the effects is estimated from the + residuals of the between regression. The degrees of freedom of the + variances are the numbers of coefficients that the regressions estimate. A + regressor that is constant within entities (the constant, or years of + education in a wage panel) or that is collinear with the others once the + means of the entities are removed (experience, if the model has year + dummies) is not in the within regression, and a regressor that is collinear + with the others in the means of the entities (the year dummies and the + constant in a balanced panel) is not in the between regression. A + ``ValueError`` is raised if either regression has no residual degrees of + freedom, since the variances cannot be estimated. + + .. versionchanged:: 7.1 + + The degrees of freedom of the variances count the coefficients that the + regressions estimate. Earlier versions counted every column of ``exog`` + (less one for the constant in the within regression), so that the + variances, theta and the estimates were slightly different for models with + time-invariant or collinear regressors, and for models without a constant. + With ``small_sample`` and an unbalanced panel, the residuals of the between + regression that the estimator of Baltagi and Chang (1994) uses are those of + the regression that weights the entities by their numbers of observations. + Earlier versions used the residuals of the unweighted regression of the + means, which gave a different variance of the effects. Models that do not + have residual degrees of freedom in the regressions raise a + ``ValueError``. Earlier versions raised an error from NumPy or SciPy, + or returned variances that were not defined. """ w = self.weights.values2d root_w = cast("linearmodels.typing.data.Float64Array", np.sqrt(w)) @@ -3079,6 +3108,15 @@ def fit( assert isinstance(demeaned_exog, PanelData) y = demeaned_dep.values2d x = demeaned_exog.values2d + # The regressors that do not vary within entities (a constant, a time-invariant + # regressor) or that are collinear once the means are removed have no + # coefficient in the within regression, and do not use a degree of freedom. + # Removing the means leaves rounding error in a regressor that does not vary, + # which is large if its units are large, so the rank is that of the regressors + # in the units that they had before the means were removed + x_norm = np.sqrt(((root_w * self.exog.values2d) ** 2).sum(0)) + x_norm[x_norm == 0.0] = 1.0 + rank_within = _lstsq(x / x_norm, y, rcond=None)[2] if self.has_constant: w_sum = w.sum() y_gm = (w * self.dependent.values2d).sum(0) / w_sum @@ -3090,28 +3128,59 @@ def fit( wybar = self.dependent.mean("entity", weights=self.weights) wxbar = self.exog.mean("entity", weights=self.weights) - params = _lstsq(np.asarray(wxbar), np.asarray(wybar), rcond=None)[0] + params, _, rank_between, _ = _lstsq( + np.asarray(wxbar), np.asarray(wybar), rcond=None + ) wu = np.asarray(wybar) - np.asarray(wxbar) @ params nobs = weps.shape[0] neffects = wu.shape[0] - nvar = x.shape[1] - sigma2_e = float(np.squeeze(weps.T @ weps)) / (nobs - nvar - neffects + 1) + df_within = nobs - neffects - rank_within + df_between = neffects - rank_between + if df_within <= 0: + raise ValueError( + "The variance of the idiosyncratic error cannot be estimated since " + "the within regression has no residual degrees of freedom. The model " + f"has {nobs} observations of {neffects} entities, and the regressors " + f"that vary within entities have a rank of {rank_within}, so that " + f"at least {neffects + rank_within + 1} observations are needed." + ) + if df_between <= 0: + raise ValueError( + "The variance of the effects cannot be estimated since the between " + "regression has no residual degrees of freedom. The model has " + f"{neffects} entities, and the means of the regressors have a rank " + f"of {rank_between}, so that at least {rank_between + 1} entities " + "are needed." + ) + sigma2_e = float(np.squeeze(weps.T @ weps)) / df_within ssr = float(np.squeeze(wu.T @ wu)) t = np.asarray(self.dependent.count("entity")) unbalanced = np.ptp(t) != 0 if small_sample and unbalanced: - ssr = float(np.squeeze((t * wu).T @ wu)) + # The between regression of the estimator of Baltagi and Chang (1994) is + # the regression of the data replaced by the means of the entities. It + # weights the means by the numbers of observations, as the trace below + # does, and the residuals of the regression of the means do not + root_t = np.sqrt(t) + xbar, ybar = np.asarray(wxbar), np.asarray(wybar) + params_t = _lstsq(root_t * xbar, root_t * ybar, rcond=None)[0] + wu_t = ybar - xbar @ params_t + ssr = float(np.squeeze((t * wu_t).T @ wu_t)) wx_df = cast("DataFrame", root_w * self.exog.dataframe) means = wx_df.groupby(level=0).transform("mean").values denom = means.T @ means sums = wx_df.groupby(level=0).sum().values num = sums.T @ sums - tr = np.trace(np.linalg.inv(denom) @ num) - sigma2_u = max(0, (ssr - (neffects - nvar) * sigma2_e) / (nobs - tr)) + # The means of the regressors can be collinear, e.g., the dummy variables + # of two periods that are observed for every entity have the same mean + # in all entities. The pseudo inverse uses the regressors that are not + # collinear + tr = np.trace(np.linalg.pinv(denom, hermitian=True) @ num) + sigma2_u = max(0, (ssr - df_between * sigma2_e) / (nobs - tr)) else: t_bar = neffects / ((1.0 / t).sum()) - sigma2_u = max(0, ssr / (neffects - nvar) - sigma2_e / t_bar) + sigma2_u = max(0, ssr / df_between - sigma2_e / t_bar) rho = sigma2_u / (sigma2_u + sigma2_e) theta = 1.0 - np.sqrt(sigma2_e / (t * sigma2_u + sigma2_e)) diff --git a/linearmodels/tests/panel/results/re-variance-components-reference.R b/linearmodels/tests/panel/results/re-variance-components-reference.R new file mode 100644 index 0000000000..644ff1dfb3 --- /dev/null +++ b/linearmodels/tests/panel/results/re-variance-components-reference.R @@ -0,0 +1,111 @@ +# Reference values for the variance components of RandomEffects, from the R +# package plm, which are read by ../test_random_effects_variance.py. +# +# Run from this directory: Rscript re-variance-components-reference.R +# Requires: plm (tested with 2.6.7 on R 4.6.1) +# Input: ../../../datasets/wage_panel/wage_panel.csv.bz2 +# Output: re-variance-components-reference.csv +# +# The model is the Swamy-Arora random effects model of plm, +# plm(..., model = "random"), whose estimator of the variance components uses the +# degrees of freedom of the regressions that estimate them: +# +# sigma2_idios = RSS(within) / (n - N - K_w) +# sigma2_id = RSS(between) / (N - K_b) - sigma2_idios / Tbar (balanced) +# +# where K_w is the number of coefficients of the within regression, which are the +# regressors that vary within the entities and are not collinear after the +# entities' means are removed, and K_b is the number of coefficients of the +# between regression, including the constant. Time-invariant regressors (black, +# hisp, educ) are not in the within regression, and the dummy variables of the +# years are collinear with the constant in the between regression of a balanced +# panel. The degrees of freedom do not count them. For an unbalanced panel plm +# uses the estimator of Baltagi and Chang (1994) (dfcor = 3), which is the +# small_sample option of RandomEffects.fit. That branch of ercomp counts the +# columns of the model matrix of the within regression for K_w, so that it counts +# a regressor that is aliased (exper with the year dummies), and this is not the +# df.residual of the within model of plm. The unbalanced cases have no aliased +# regressor. +# +# Stata's xtreg, re uses the same degrees of freedom: for the Cornwell and Rupert +# wage panel (A. C. Cameron's log of psidextract.dta), xtreg lwage exp exp2 wks ed, +# re reports sigma_e = .15220316, which is also the sigma_e of xtreg lwage exp exp2 +# wks, fe, where ed is omitted because of collinearity (n - N - 3 degrees of +# freedom). The same model has sigma_e = 0.1522032 in plm. +# +# Cases (the unbalanced cases drop observations by a fixed rule): +# balanced-notebook the model of the panel_examples notebook: year dummies, +# time-invariant regressors, and exper, which is +# collinear with the year dummies in the within regression +# balanced-time-invariant time-invariant regressors without year dummies +# balanced-time-varying regressors that all vary within the entities +# unbalanced-year-dummies year dummies and time-invariant regressors, without exper +# unbalanced-time-invariant as balanced-time-invariant + +suppressMessages(library(plm)) + +cat(R.version.string, "and plm", as.character(packageVersion("plm")), "\n") + +# CSV files with "\n" line endings on every platform +write_csv <- function(x, path, ...) { + con <- file(path, open = "wb") + on.exit(close(con)) + write.csv(x, con, row.names = FALSE, eol = "\n", ...) +} + +wages <- read.csv(bzfile("../../../datasets/wage_panel/wage_panel.csv.bz2")) + +# The observations that the unbalanced cases drop: 588 of 4,360, so that each +# entity has between 6 and 8 observations and the means of the year dummies differ +# between the entities +drop <- ((wages$nr * 7 + wages$year * 13) %% 5 == 0) & ((wages$nr + wages$year) %% 3 != 0) +cat("unbalanced cases drop", sum(drop), "of", nrow(wages), "observations\n") + +formulas <- list( + notebook = lwage ~ black + hisp + exper + expersq + married + educ + union + + factor(year), + `year-dummies` = lwage ~ black + hisp + expersq + married + educ + union + + factor(year), + `time-invariant` = lwage ~ black + hisp + exper + expersq + married + educ + union, + `time-varying` = lwage ~ expersq + married + union +) +cases <- list( + list(name = "balanced-notebook", data = wages, formula = "notebook"), + list(name = "balanced-time-invariant", data = wages, formula = "time-invariant"), + list(name = "balanced-time-varying", data = wages, formula = "time-varying"), + list(name = "unbalanced-year-dummies", data = wages[!drop, ], + formula = "year-dummies"), + list(name = "unbalanced-time-invariant", data = wages[!drop, ], + formula = "time-invariant") +) + +rows <- list() +for (case in cases) { + pd <- pdata.frame(case$data, index = c("nr", "year")) + fit <- plm(formulas[[case$formula]], data = pd, model = "random") + components <- ercomp(fit) + theta <- as.numeric(components$theta) + est <- coef(fit) + se <- sqrt(diag(vcov(fit))) + within <- plm(formulas[[case$formula]], data = pd, model = "within") + between <- plm(formulas[[case$formula]], data = pd, model = "between") + rows[[case$name]] <- rbind( + data.frame(case = case$name, quantity = "nobs", value = nobs(fit)), + data.frame(case = case$name, quantity = "df_within", value = df.residual(within)), + data.frame(case = case$name, quantity = "df_between", value = df.residual(between)), + data.frame(case = case$name, quantity = "sigma2_idios", + value = unname(components$sigma2["idios"])), + data.frame(case = case$name, quantity = "sigma2_id", + value = unname(components$sigma2["id"])), + data.frame(case = case$name, quantity = "theta_mean", value = mean(theta)), + data.frame(case = case$name, quantity = paste0("coef", seq_along(est) - 1), + value = unname(est)), + data.frame(case = case$name, quantity = paste0("se", seq_along(est) - 1), + value = unname(se)) + ) + cat(sprintf("%-27s nobs %5d df within %5d df between %4d idios %.7f id %.7f\n", + case$name, nobs(fit), df.residual(within), df.residual(between), + components$sigma2["idios"], components$sigma2["id"])) +} +write_csv(do.call(rbind, rows), "re-variance-components-reference.csv") +cat("wrote re-variance-components-reference.csv\n") diff --git a/linearmodels/tests/panel/results/re-variance-components-reference.csv b/linearmodels/tests/panel/results/re-variance-components-reference.csv new file mode 100644 index 0000000000..92c414c316 --- /dev/null +++ b/linearmodels/tests/panel/results/re-variance-components-reference.csv @@ -0,0 +1,129 @@ +"case","quantity","value" +"balanced-notebook","nobs",4360 +"balanced-notebook","df_within",3805 +"balanced-notebook","df_between",537 +"balanced-notebook","sigma2_idios",0.123193989960515 +"balanced-notebook","sigma2_id",0.105367200462246 +"balanced-notebook","theta_mean",0.642910879665865 +"balanced-notebook","coef0",0.0235864096542448 +"balanced-notebook","coef1",-0.139376734111937 +"balanced-notebook","coef2",0.0217317306549108 +"balanced-notebook","coef3",0.105754523245461 +"balanced-notebook","coef4",-0.00472394273755648 +"balanced-notebook","coef5",0.0639860359475824 +"balanced-notebook","coef6",0.091876272649444 +"balanced-notebook","coef7",0.106134422744568 +"balanced-notebook","coef8",0.0404619940706127 +"balanced-notebook","coef9",0.0309211619343808 +"balanced-notebook","coef10",0.0202806373503664 +"balanced-notebook","coef11",0.043118684873118 +"balanced-notebook","coef12",0.0578154064233405 +"balanced-notebook","coef13",0.0919475382668158 +"balanced-notebook","coef14",0.134928886286404 +"balanced-notebook","se0",0.150668257591444 +"balanced-notebook","se1",0.0477228164328931 +"balanced-notebook","se2",0.042606290027916 +"balanced-notebook","se3",0.0153668157460012 +"balanced-notebook","se4",0.00068949694503348 +"balanced-notebook","se5",0.0167742437355273 +"balanced-notebook","se6",0.0106597040958363 +"balanced-notebook","se7",0.0178538555251242 +"balanced-notebook","se8",0.024694610713068 +"balanced-notebook","se9",0.0323416128411864 +"balanced-notebook","se10",0.0415819882406935 +"balanced-notebook","se11",0.0513163474941151 +"balanced-notebook","se12",0.0612323120313417 +"balanced-notebook","se13",0.0712292614222332 +"balanced-notebook","se14",0.0813135284535609 +"balanced-time-invariant","nobs",4360 +"balanced-time-invariant","df_within",3811 +"balanced-time-invariant","df_between",537 +"balanced-time-invariant","sigma2_idios",0.123380320307742 +"balanced-time-invariant","sigma2_id",0.105343909168843 +"balanced-time-invariant","theta_mean",0.642640933868348 +"balanced-time-invariant","coef0",-0.107464203974044 +"balanced-time-invariant","coef1",-0.144130691111868 +"balanced-time-invariant","coef2",0.0201510730060844 +"balanced-time-invariant","coef3",0.112119493522384 +"balanced-time-invariant","coef4",-0.00406885475618792 +"balanced-time-invariant","coef5",0.0627951179727595 +"balanced-time-invariant","coef6",0.101224614698513 +"balanced-time-invariant","coef7",0.107378852592373 +"balanced-time-invariant","se0",0.110705725594299 +"balanced-time-invariant","se1",0.0476148274269489 +"balanced-time-invariant","se2",0.042601124174701 +"balanced-time-invariant","se3",0.00826087205583039 +"balanced-time-invariant","se4",0.000591825600027988 +"balanced-time-invariant","se5",0.0167728540561277 +"balanced-time-invariant","se6",0.00891328987417564 +"balanced-time-invariant","se7",0.0178300147703687 +"balanced-time-varying","nobs",4360 +"balanced-time-varying","df_within",3812 +"balanced-time-varying","df_between",541 +"balanced-time-varying","sigma2_idios",0.129581486663583 +"balanced-time-varying","sigma2_id",0.122493126638203 +"balanced-time-varying","theta_mean",0.658254629973685 +"balanced-time-varying","coef0",1.40520815095454 +"balanced-time-varying","coef1",0.00319171996193454 +"balanced-time-varying","coef2",0.131555169153797 +"balanced-time-varying","coef3",0.103450759968374 +"balanced-time-varying","se0",0.0192679695028832 +"balanced-time-varying","se1",0.000181336837300841 +"balanced-time-varying","se2",0.0168817492705327 +"balanced-time-varying","se3",0.0184859713861548 +"unbalanced-year-dummies","nobs",3772 +"unbalanced-year-dummies","df_within",3217 +"unbalanced-year-dummies","df_between",531 +"unbalanced-year-dummies","sigma2_idios",0.128302344184146 +"unbalanced-year-dummies","sigma2_id",0.108358317305489 +"unbalanced-year-dummies","theta_mean",0.618184470344336 +"unbalanced-year-dummies","coef0",0.605548790479954 +"unbalanced-year-dummies","coef1",-0.125821461978439 +"unbalanced-year-dummies","coef2",0.0174929636180146 +"unbalanced-year-dummies","coef3",-0.00130071500004915 +"unbalanced-year-dummies","coef4",0.0837165981326888 +"unbalanced-year-dummies","coef5",0.0646943293668789 +"unbalanced-year-dummies","coef6",0.127259889974332 +"unbalanced-year-dummies","coef7",0.112573187254709 +"unbalanced-year-dummies","coef8",0.188019604193127 +"unbalanced-year-dummies","coef9",0.242700240753059 +"unbalanced-year-dummies","coef10",0.333987251802706 +"unbalanced-year-dummies","coef11",0.394081555168255 +"unbalanced-year-dummies","coef12",0.470921208312212 +"unbalanced-year-dummies","coef13",0.562894669724611 +"unbalanced-year-dummies","se0",0.127130321119421 +"unbalanced-year-dummies","se1",0.0491356466494286 +"unbalanced-year-dummies","se2",0.0439467291801474 +"unbalanced-year-dummies","se3",0.00053586489757734 +"unbalanced-year-dummies","se4",0.0182603878294297 +"unbalanced-year-dummies","se5",0.0102145931419279 +"unbalanced-year-dummies","se6",0.0195234302249132 +"unbalanced-year-dummies","se7",0.0242011174360354 +"unbalanced-year-dummies","se8",0.0252961420112048 +"unbalanced-year-dummies","se9",0.0282487685579614 +"unbalanced-year-dummies","se10",0.0324318401001503 +"unbalanced-year-dummies","se11",0.0383367319727384 +"unbalanced-year-dummies","se12",0.0457883832019504 +"unbalanced-year-dummies","se13",0.0547496534149072 +"unbalanced-time-invariant","nobs",3772 +"unbalanced-time-invariant","df_within",3223 +"unbalanced-time-invariant","df_between",537 +"unbalanced-time-invariant","sigma2_idios",0.128427843646611 +"unbalanced-time-invariant","sigma2_id",0.107626857054882 +"unbalanced-time-invariant","theta_mean",0.6169198975766 +"unbalanced-time-invariant","coef0",-0.125812807588583 +"unbalanced-time-invariant","coef1",-0.151894553284766 +"unbalanced-time-invariant","coef2",0.0157023105776524 +"unbalanced-time-invariant","coef3",0.115990971886776 +"unbalanced-time-invariant","coef4",-0.00431687671630469 +"unbalanced-time-invariant","coef5",0.0683264347224204 +"unbalanced-time-invariant","coef6",0.100948952871547 +"unbalanced-time-invariant","coef7",0.124414893965101 +"unbalanced-time-invariant","se0",0.113929832936069 +"unbalanced-time-invariant","se1",0.0487220333669172 +"unbalanced-time-invariant","se2",0.0435788643259302 +"unbalanced-time-invariant","se3",0.00910200119006807 +"unbalanced-time-invariant","se4",0.000653075903819231 +"unbalanced-time-invariant","se5",0.0182689960010359 +"unbalanced-time-invariant","se6",0.00913644762476636 +"unbalanced-time-invariant","se7",0.0193869180875927 diff --git a/linearmodels/tests/panel/test_random_effects_no_df.py b/linearmodels/tests/panel/test_random_effects_no_df.py new file mode 100644 index 0000000000..3bcb10ad0b --- /dev/null +++ b/linearmodels/tests/panel/test_random_effects_no_df.py @@ -0,0 +1,150 @@ +""" +RandomEffects where the regressions of the variance components have no degrees of freedom + +The variance of the idiosyncratic error is the residual sum of squares of the within +regression divided by its degrees of freedom, nobs - N - rank, where rank is the rank +of the regressors once the means of the entities are removed. The variance of the +effects is estimated from the between regression, which has N - rank observations of +the means of the regressors left. The variances are not defined if there are no +degrees of freedom, e.g., if there is a single observation of each entity or fewer +entities than regressors. These models raised an error from the division, from NumPy or +from SciPy, or returned an infinite variance of the effects. +""" + +import itertools + +import numpy as np +from numpy.testing import assert_allclose +import pandas as pd +import pytest + +from linearmodels.panel import RandomEffects + + +def make_panel(counts, nvar, const, seed=0): + rg = np.random.default_rng(seed) + entities = np.repeat(np.arange(len(counts)), counts) + times = np.concatenate([np.arange(count) for count in counts]) + index = pd.MultiIndex.from_arrays([entities, times], names=["entity", "time"]) + x = pd.DataFrame( + rg.standard_normal((entities.shape[0], nvar)), + index=index, + columns=[f"x{i}" for i in range(nvar)], + ) + if const: + x.insert(0, "const", 1.0) + y = pd.Series(rg.standard_normal(entities.shape[0]), index=index, name="y") + return y, x + + +def degrees_of_freedom(x): + """The residual degrees of freedom of the within and the between regressions""" + nobs = x.shape[0] + means = x.groupby(level=0).mean() + demeaned = x - x.groupby(level=0).transform("mean") + within = nobs - means.shape[0] - np.linalg.matrix_rank(demeaned.to_numpy()) + between = means.shape[0] - np.linalg.matrix_rank(means.to_numpy()) + return within, between + + +def lsdv_variance(y, x): + """Residual variance of the regression on the dummy variables of the entities""" + codes = pd.factorize(x.index.get_level_values(0))[0] + dummies = np.eye(codes.max() + 1)[codes] + design = np.column_stack([x.to_numpy(), dummies]) + resid = y.to_numpy() - design @ np.linalg.lstsq(design, y.to_numpy(), rcond=None)[0] + return resid @ resid / (design.shape[0] - np.linalg.matrix_rank(design)) + + +WITHIN = "variance of the idiosyncratic error cannot be estimated" +BETWEEN = "variance of the effects cannot be estimated" +BALANCED = [[t] * n for n in (1, 2, 3, 4) for t in (1, 2, 3)] +UNBALANCED = [[1, 2, 1], [1, 1, 2], [2, 2, 1], [1, 3, 2, 1], [3, 1, 4, 2, 2]] + + +def full_rank_designs(): + """The designs where the regressors are not collinear, as the model requires""" + for counts, nvar, const in itertools.product( + BALANCED + UNBALANCED, [1, 2, 3], [False, True] + ): + x = make_panel(counts, nvar, const)[1] + if np.linalg.matrix_rank(x.to_numpy()) == x.shape[1]: + yield counts, nvar, const + + +DESIGNS = list(full_rank_designs()) +DESIGN_IDS = [f"{c}-{k}-{'const' if const else 'noconst'}" for c, k, const in DESIGNS] + + +@pytest.mark.parametrize("small_sample", [False, True]) +@pytest.mark.parametrize(("counts", "nvar", "const"), DESIGNS, ids=DESIGN_IDS) +def test_variances_need_degrees_of_freedom(counts, nvar, const, small_sample): + y, x = make_panel(counts, nvar, const) + within, between = degrees_of_freedom(x) + mod = RandomEffects(y, x) + if within <= 0: + with pytest.raises(ValueError, match=WITHIN): + mod.fit(small_sample=small_sample) + elif between <= 0: + with pytest.raises(ValueError, match=BETWEEN): + mod.fit(small_sample=small_sample) + else: + res = mod.fit(small_sample=small_sample) + variances = res.variance_decomposition + assert np.all(np.isfinite(variances)) + assert variances["Residual"] > 0 + assert variances["Effects"] >= 0 + assert np.all(np.isfinite(res.params)) + assert_allclose(variances["Residual"], lsdv_variance(y, x)) + + +def test_each_cell_of_the_grid_is_used(): + # Some designs of the grid fail for each reason, and some are estimated + reasons = [] + for counts, nvar, const in DESIGNS: + within, between = degrees_of_freedom(make_panel(counts, nvar, const)[1]) + reasons.append("within" if within <= 0 else "between" if between <= 0 else "ok") + assert {"within", "between", "ok"} == set(reasons) + assert all(reasons.count(reason) >= 10 for reason in set(reasons)) + + +@pytest.mark.parametrize("small_sample", [False, True]) +def test_one_observation_of_each_entity(small_sample): + y, x = make_panel([1] * 10, 1, True) + # The 10 means leave nothing of the 10 observations, and the regressor does not vary + with pytest.raises(ValueError, match=WITHIN) as exc_info: + RandomEffects(y, x).fit(small_sample=small_sample) + assert "at least 11 observations are needed" in str(exc_info.value) + assert "10 observations of 10 entities" in str(exc_info.value) + + +@pytest.mark.parametrize("small_sample", [False, True]) +def test_as_many_regressors_as_the_observations_of_the_entities_leave(small_sample): + # 3 entities with 2 observations: 3 regressors use the 3 degrees of freedom + y, x = make_panel([2, 2, 2], 3, False) + with pytest.raises(ValueError, match=WITHIN) as exc_info: + RandomEffects(y, x).fit(small_sample=small_sample) + assert "rank of 3" in str(exc_info.value) + assert "at least 7 observations are needed" in str(exc_info.value) + + +@pytest.mark.parametrize("small_sample", [False, True]) +@pytest.mark.parametrize("counts", [[10], [3, 3], [4, 4]]) +def test_not_more_entities_than_the_rank_of_the_means(counts, small_sample): + # The constant and a regressor have a rank of 2 in the means of 1 or 2 entities + y, x = make_panel(counts, 1, True) + with pytest.raises(ValueError, match=BETWEEN) as exc_info: + RandomEffects(y, x).fit(small_sample=small_sample) + rank = min(len(counts), 2) + assert f"{len(counts)} entities" in str(exc_info.value) + assert f"at least {rank + 1} entities are needed" in str(exc_info.value) + + +@pytest.mark.parametrize("small_sample", [False, True]) +def test_one_more_entity_than_the_rank_of_the_means(small_sample): + # The smallest panel that can be estimated: 3 entities, a constant and a regressor + y, x = make_panel([3, 3, 3], 1, True) + assert degrees_of_freedom(x) == (9 - 3 - 1, 1) + res = RandomEffects(y, x).fit(small_sample=small_sample) + assert np.all(np.isfinite(res.variance_decomposition)) + assert_allclose(res.variance_decomposition["Residual"], lsdv_variance(y, x)) diff --git a/linearmodels/tests/panel/test_random_effects_variance.py b/linearmodels/tests/panel/test_random_effects_variance.py new file mode 100644 index 0000000000..95b4dc40d4 --- /dev/null +++ b/linearmodels/tests/panel/test_random_effects_variance.py @@ -0,0 +1,333 @@ +""" +The variance components of RandomEffects + +The idiosyncratic variance is the residual variance of the within regression, and +the variance of the effects is estimated from the residuals of the between +regression. The degrees of freedom of the variances are those of the regressions: +the within regression has the coefficients of the regressors that vary within +entities and are not collinear once the means of the entities are removed, and the +between regression has the coefficients of the regressors that are not collinear +across entities. A time-invariant regressor, such as education in a wage panel, is +not in the within regression, and neither is a regressor that is collinear with the +dummy variables of the time periods within entities, such as experience. A +dummy variable of the time periods can be collinear with the constant in the between +regression. The degrees of freedom do not count them. + +The references are + +* the residual variance of the regression on the dummy variables of the entities, + which is the idiosyncratic variance of Swamy and Arora; +* the estimators written out in NumPy using the ranks of the regressions; +* plm, in results/re-variance-components-reference.csv, for the wage panel of Vella + and Verbeek. The R script that produced it is + results/re-variance-components-reference.R. Stata's xtreg, re counts the + coefficients in the same way: sigma_e of a model with a time-invariant regressor is + that of xtreg, fe, which omits the regressor. +""" + +import os + +import numpy as np +from numpy.testing import assert_allclose +import pandas as pd +import pytest +from scipy.linalg import qr + +from linearmodels.datasets import wage_panel +from linearmodels.panel import PanelOLS, RandomEffects + +RESULTS = os.path.join(os.path.split(os.path.abspath(__file__))[0], "results") +NPERIOD = 6 +DUMMIES = [f"d{p}" for p in range(1, NPERIOD)] +DESIGNS = { + "time-varying": ["const", "x1"], + "no-constant": ["x1"], + "time-invariant": ["const", "x1", "z"], + "time-invariant-no-constant": ["x1", "z"], + # exper is collinear with the dummies once the means of the entities are removed + "collinear-with-dummies": ["const", "x1", "z", "exper", *DUMMIES], +} + + +def make_data(balance, seed=0, nentity=40): + """ + A panel with a time-invariant regressor (z) and one that is collinear with the + dummies of the time periods within entities (exper), as the experience in a + wage panel. The unbalanced panels drop observations from the last three periods + ("collinear-means") or from all periods + """ + rg = np.random.default_rng(seed) + index = pd.MultiIndex.from_product( + [np.arange(nentity), np.arange(NPERIOD)], names=["entity", "time"] + ) + entity = index.codes[0] + period = index.codes[1] + nobs = nentity * NPERIOD + effect = rg.standard_normal(nentity)[entity] + x1 = 0.5 * effect + rg.standard_normal(nobs) + z = rg.standard_normal(nentity)[entity] + exper = period + rg.integers(0, 5, nentity)[entity] + shock = (1 + 0.3 * np.abs(x1)) * rg.standard_normal(nobs) + data = pd.DataFrame( + { + "y": 1 + 0.8 * x1 + 0.5 * z + 0.05 * exper + 1.2 * effect + shock, + "const": 1.0, + "x1": x1, + "z": z, + "exper": exper, + "w": rg.chisquare(4, nobs) / 4 + 0.25, + }, + index=index, + ) + for p in range(1, NPERIOD): + data[f"d{p}"] = (period == p).astype(float) + if balance != "balanced": + keep = rg.random((nentity, NPERIOD)) > 0.3 + if balance == "unbalanced-collinear-means": + # The dummies of the first periods have the same mean in all entities + keep[:, :3] = True + for i in range(nentity): + keep[i, rg.choice(NPERIOD, 3, replace=False)] = True + data = data.loc[keep.ravel()] + return data + + +@pytest.fixture( + params=["balanced", "unbalanced", "unbalanced-collinear-means"], scope="module" +) +def data(request): + return make_data(request.param) + + +def lsdv_variance(data, columns, weights): + """ + The residual variance of the regression on the dummy variables of the entities + + The regressors of the model that are not identified by the dummy variables + (the time-invariant ones, which they span) do not reduce the degrees of freedom. + """ + dummies = pd.get_dummies(data.index.get_level_values(0)).astype(float) + regressors = [c for c in columns if c != "const"] + exog = np.column_stack([data[regressors].to_numpy(), dummies.to_numpy()]) + root_w = np.ones(data.shape[0]) + if weights is not None: + root_w = np.sqrt((weights / weights.mean()).to_numpy()) + wexog, wy = root_w[:, None] * exog, root_w * data["y"].to_numpy() + beta, _, rank, _ = np.linalg.lstsq(wexog, wy, rcond=None) + return float(np.sum((wy - wexog @ beta) ** 2)) / (data.shape[0] - rank) + + +def swar_reference(y, x, weights=None, small_sample=False): + """ + The estimators of Swamy and Arora written out, with the ranks of the regressions + + Returns the idiosyncratic variance, the variance of the effects and the numbers + of coefficients of the within and the between regressions + """ + entity = y.index.get_level_values(0) + weights = pd.Series(1.0, index=y.index) if weights is None else weights + weights = weights / weights.mean() + total = weights.groupby(entity).transform("sum") + + def weighted_means(frame): + weighted = frame.mul(weights, axis=0) + return weighted.groupby(entity).transform("sum").div(total, axis=0) + + root_w = np.sqrt(weights.to_numpy()) + xd = root_w[:, None] * (x - weighted_means(x)).to_numpy() + yd = root_w * (y - weighted_means(y.to_frame()).iloc[:, 0]).to_numpy() + beta, _, rank_within, _ = np.linalg.lstsq(xd, yd, rcond=None) + nobs, nentity = x.shape[0], entity.nunique() + sigma2_eps = float(np.sum((yd - xd @ beta) ** 2)) / (nobs - nentity - rank_within) + # The regression of the means of the entities + xb = weighted_means(x).groupby(entity).first().to_numpy() + yb = weighted_means(y.to_frame()).iloc[:, 0].groupby(entity).first().to_numpy() + count = y.groupby(entity).size().to_numpy().astype(float) + rank_between = int(np.linalg.matrix_rank(xb)) + if not small_sample: + resid = yb - xb @ np.linalg.lstsq(xb, yb, rcond=None)[0] + tbar = nentity / np.sum(1 / count) + sigma2_effects = float(np.sum(resid**2)) / (nentity - rank_between) + sigma2_effects = max(0.0, sigma2_effects - sigma2_eps / tbar) + return sigma2_eps, sigma2_effects, rank_within, rank_between + # Baltagi and Chang (1994) for unbalanced panels: the means are weighted by the + # number of observations. This is the regression of the data replaced by the + # means of the entities. The trace uses columns that are not collinear + root_t = np.sqrt(count) + coef = np.linalg.lstsq(root_t[:, None] * xb, root_t * yb, rcond=None)[0] + rss = float(np.sum(count * (yb - xb @ coef) ** 2)) + basis = qr(xb, pivoting=True, mode="economic")[2][:rank_between] + xr = xb[:, basis] + trace = np.trace( + np.linalg.inv((count[:, None] * xr).T @ xr) + @ ((count[:, None] * xr).T @ (count[:, None] * xr)) + ) + sigma2_effects = (rss - (nentity - rank_between) * sigma2_eps) / (nobs - trace) + return sigma2_eps, max(0.0, float(sigma2_effects)), rank_within, rank_between + + +@pytest.mark.parametrize("weighted", [False, True], ids=["unweighted", "weighted"]) +@pytest.mark.parametrize("design", list(DESIGNS)) +def test_idiosyncratic_variance_against_dummy_regression(data, design, weighted): + columns = DESIGNS[design] + weights = data["w"] if weighted else None + res = RandomEffects(data["y"], data[columns], weights=weights).fit() + expected = lsdv_variance(data, columns, weights) + assert_allclose(res.variance_decomposition["Residual"], expected) + + +@pytest.mark.parametrize("design", ["time-invariant", "time-invariant-no-constant"]) +def test_idiosyncratic_variance_is_that_of_the_fixed_effects(data, design): + # Stata: sigma_e of xtreg, re with a time-invariant regressor is the sigma_e + # of xtreg, fe, which omits the regressor + columns = DESIGNS[design] + re_res = RandomEffects(data["y"], data[columns]).fit() + fe_columns = [c for c in columns if c not in ("const", "z")] + fe_res = PanelOLS(data["y"], data[fe_columns], entity_effects=True).fit() + assert_allclose(re_res.variance_decomposition["Residual"], fe_res.s2) + + +@pytest.mark.parametrize("scale", [1e3, 1e6, 1e8]) +def test_variance_does_not_depend_on_the_units_of_a_regressor_that_does_not_vary( + data, scale +): + # Removing the means of a regressor that does not vary within entities leaves + # rounding error, which is large next to the other regressors if its units are + # large, and it is not a coefficient of the within regression + columns = DESIGNS["time-invariant"] + scaled = data.assign(z=data["z"] * np.pi * scale) + res = RandomEffects(scaled["y"], scaled[columns]).fit() + expected = lsdv_variance(scaled, columns, None) + assert_allclose(res.variance_decomposition["Residual"], expected) + + +@pytest.mark.parametrize("small_sample", [False, True]) +@pytest.mark.parametrize("design", list(DESIGNS)) +def test_variance_components_against_estimators_written_out(data, design, small_sample): + columns = DESIGNS[design] + res = RandomEffects(data["y"], data[columns]).fit(small_sample=small_sample) + sigma2_eps, sigma2_effects, *_ = swar_reference( + data["y"], data[columns], small_sample=small_sample + ) + decomposition = res.variance_decomposition + assert_allclose(decomposition["Residual"], sigma2_eps) + assert_allclose(decomposition["Effects"], sigma2_effects) + rho = sigma2_effects / (sigma2_effects + sigma2_eps) + assert_allclose(decomposition["Percent due to Effects"], rho) + theta = res.theta["theta"] + count = data.groupby(level=0).size().reindex(theta.index) + expected_theta = 1 - np.sqrt(sigma2_eps / (count * sigma2_effects + sigma2_eps)) + assert_allclose(theta, expected_theta) + + +@pytest.mark.parametrize("design", list(DESIGNS)) +def test_weighted_variance_components_against_estimators_written_out(data, design): + columns = DESIGNS[design] + res = RandomEffects(data["y"], data[columns], weights=data["w"]).fit() + sigma2_eps, sigma2_effects, *_ = swar_reference( + data["y"], data[columns], weights=data["w"] + ) + assert_allclose(res.variance_decomposition["Residual"], sigma2_eps) + assert_allclose(res.variance_decomposition["Effects"], sigma2_effects) + + +def test_the_designs_have_regressors_that_the_regressions_do_not_estimate(): + # The tests are only informative if the numbers of coefficients of the + # regressions differ from the number of regressors. The within regression + # has x1 and the dummies: z and the constant do not vary within entities, and + # exper is collinear with the dummies. In the means of a balanced panel the + # dummies are the constant 1/T. If periods 0, 1 and 2 are observed for every + # entity, the means of their dummies are equal (1/T_i), and so the constant is + # 3 times the mean of the dummy of period 1 plus the means of the other dummies + columns = DESIGNS["collinear-with-dummies"] + for balance, between in [ + ("balanced", 4), + ("unbalanced", len(columns)), + ("unbalanced-collinear-means", len(columns) - 2), + ]: + sample = make_data(balance) + ranks = swar_reference(sample["y"], sample[columns])[2:] + assert ranks == (1 + len(DUMMIES), between), balance + + +def wage_data(unbalanced): + """The wage panel of the notebook, and its unbalanced version in the R script""" + data = wage_panel.load() + if unbalanced: + year, nr = data["year"], data["nr"] + drop = ((nr * 7 + year * 13) % 5 == 0) & ((nr + year) % 3 != 0) + data = data.loc[~drop] + year = pd.Categorical(data["year"]) + data = data.set_index(["nr", "year"]) + data["year"] = year + data["const"] = 1.0 + return data + + +TIME_INVARIANT = ["black", "hisp", "exper", "expersq", "married", "educ", "union"] +PLM_COLUMNS = { + "notebook": ["const", *TIME_INVARIANT, "year"], + "year-dummies": [ + "const", + "black", + "hisp", + "expersq", + "married", + "educ", + "union", + "year", + ], + "time-invariant": ["const", *TIME_INVARIANT], + "time-varying": ["const", "expersq", "married", "union"], +} +PLM_CASES = [ + ("balanced", "notebook"), + ("balanced", "time-invariant"), + ("balanced", "time-varying"), + ("unbalanced", "year-dummies"), + ("unbalanced", "time-invariant"), +] + + +@pytest.fixture(scope="module") +def plm_reference(): + path = os.path.join(RESULTS, "re-variance-components-reference.csv") + reference = pd.read_csv(path) + return { + case: group.set_index("quantity")["value"] + for case, group in reference.groupby("case") + } + + +@pytest.mark.parametrize(("balance", "model"), PLM_CASES) +def test_random_effects_against_plm(plm_reference, balance, model): + expected = plm_reference[f"{balance}-{model}"] + unbalanced = balance == "unbalanced" + data = wage_data(unbalanced) + mod = RandomEffects(data["lwage"], data[PLM_COLUMNS[model]]) + # plm uses the estimator of Baltagi and Chang (1994) for an unbalanced panel + res = mod.fit(small_sample=unbalanced) + assert res.nobs == expected["nobs"] + decomposition = res.variance_decomposition + assert_allclose(decomposition["Residual"], expected["sigma2_idios"], rtol=1e-8) + assert_allclose(decomposition["Effects"], expected["sigma2_id"], rtol=1e-8) + # plm reports the theta of each observation + entities = mod.dependent.index.get_level_values(0) + theta = res.theta["theta"].reindex(entities) + assert_allclose(theta.mean(), expected["theta_mean"], rtol=1e-8) + nparam = res.params.shape[0] + coefficients = expected[[f"coef{i}" for i in range(nparam)]].to_numpy() + std_errors = expected[[f"se{i}" for i in range(nparam)]].to_numpy() + assert_allclose(res.params.to_numpy(), coefficients, rtol=1e-7, atol=1e-10) + assert_allclose(res.std_errors.to_numpy(), std_errors, rtol=1e-7) + # The estimators written out use the numbers of coefficients that plm does + sigma2_eps, sigma2_effects, rank_within, rank_between = swar_reference( + mod.dependent.dataframe["lwage"], + mod.exog.dataframe, + small_sample=unbalanced, + ) + nentity = entities.nunique() + assert expected["df_within"] == res.nobs - nentity - rank_within + assert expected["df_between"] == nentity - rank_between + assert_allclose(sigma2_eps, expected["sigma2_idios"], rtol=1e-8) + assert_allclose(sigma2_effects, expected["sigma2_id"], rtol=1e-8)