diff --git a/examples/asset-pricing_examples.ipynb b/examples/asset-pricing_examples.ipynb
index d24f82ac8c..8a239f633c 100644
--- a/examples/asset-pricing_examples.ipynb
+++ b/examples/asset-pricing_examples.ipynb
@@ -24,7 +24,7 @@
    "source": [
     "## Import data\n",
     "\n",
-    "The data used comes from Ken French\"s [website](http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html) and includes 4 factor returns, the excess market, the size factor, the value factor and the momentum factor.  The available test portfolios include the 12 industry portfolios, a subset of the size-value two way sort, and a subset of the size-momentum two way sort. "
+    "The data used comes from Ken French's [website](http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html) and includes 4 factor returns, the excess market, the size factor, the value factor and the momentum factor.  The available test portfolios include the 12 industry portfolios, a subset of the size-value two way sort, and a subset of the size-momentum two way sort. "
    ]
   },
   {
@@ -93,7 +93,7 @@
    "source": [
     "The factor set is expanded to include both the size and the value factors.\n",
     "\n",
-    "While the extra factors lower the J statistic and increases the $R^2$, the model is still easily rejected."
+    "While the extra factors lower the J statistic and increase the $R^2$, the model is still easily rejected."
    ]
   },
   {
@@ -191,7 +191,7 @@
    "source": [
     "## 2-Step Estimation\n",
     "\n",
-    "When one of more factors are not returns on traded assets, it is necessary to use a two-step procedure (or GMM).  In the 2-step estimator, the first estimates factor loadings and the second uses the factor loadings to estimate the risk premia.  \n",
+    "When one or more factors are not returns on traded assets, it is necessary to use a two-step procedure (or GMM).  In the 2-step estimator, the first estimates factor loadings and the second uses the factor loadings to estimate the risk premia.  \n",
     "\n",
     "Here all four factors are used to attempt to price the size-momentum portfolios."
    ]
@@ -226,7 +226,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The identification requirements for this model are the the $\\beta$s have unique variation -- this requires some cross-sectional differences in exposures and that the loadings are not excessively correlated.  Since these test portfolios do not attempt to sort on value, it is relatively non-dispersed and also correlated with both the market and size. This might make the inference unreliable.  "
+    "The identification requirements for this model are that the $\\beta$s have unique variation -- this requires some cross-sectional differences in exposures and that the loadings are not excessively correlated.  Since these test portfolios do not attempt to sort on value, it is relatively non-dispersed and also correlated with both the market and size. This might make the inference unreliable.  "
    ]
   },
   {
@@ -265,7 +265,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The risk-free rate can be optionally estimated.  This is useful in case excess returns are not available of if the risk-free rate used to construct the excess returns might be misspecified.\n",
+    "The risk-free rate can be optionally estimated.  This is useful in case excess returns are not available or if the risk-free rate used to construct the excess returns might be misspecified.\n",
     "\n",
     "Here the estimated risk-free rate is small and insignificant and has little impact on the model J-statistic."
    ]
@@ -316,7 +316,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The final estimator is the GMM estimator which is similar to estimating the 2-step model in a single step.  In practice the GMM estimator is estimated at least twice, once to get an consistent estimate of the covariance of the moment conditions and the second time to efficiently estimate parameters. \n",
+    "The final estimator is the GMM estimator which is similar to estimating the 2-step model in a single step.  In practice the GMM estimator is estimated at least twice, once to get a consistent estimate of the covariance of the moment conditions and the second time to efficiently estimate parameters. \n",
     "\n",
     "The GMM estimator does not have a closed form and so a non-linear optimizer is required.  The default output prints the progress every 10 iterations.  Here the model is fit twice, which is the standard method to implement efficient GMM."
    ]
diff --git a/examples/asset-pricing_formulas.ipynb b/examples/asset-pricing_formulas.ipynb
index 6d8d5e34d8..7b1c679401 100644
--- a/examples/asset-pricing_formulas.ipynb
+++ b/examples/asset-pricing_formulas.ipynb
@@ -6,7 +6,7 @@
    "source": [
     "# Using Formulas\n",
     "\n",
-    "The asset pricing model estimators all all formulas.  Since the models have multiple dependent variables (test portfolios) as well as multiple independent variables (factors), the standard [formulaic](https://github.com/matthewwardrop/formulaic/) syntax needs to be modified. There are two methods to use formulas.  The first specified both the test portfolio and the factors.  The second specifies only the factors and the test portfolios are passed using an optional keyword argument. The second syntax exists since in many models the number of test portfolios might be large and interest is usually in modifying the factors.\n",
+    "The asset pricing model estimators all support formulas.  Since the models have multiple dependent variables (test portfolios) as well as multiple independent variables (factors), the standard [formulaic](https://github.com/matthewwardrop/formulaic/) syntax needs to be modified. There are two methods to use formulas.  The first specifies both the test portfolio and the factors.  The second specifies only the factors and the test portfolios are passed using an optional keyword argument. The second syntax exists since in many models the number of test portfolios might be large and interest is usually in modifying the factors.\n",
     "\n",
     "## Available Syntax \n",
     "\n",
@@ -36,7 +36,7 @@
    "source": [
     "## Import data and transform to be excess returns\n",
     "\n",
-    "The data used comes from Ken French\"s [website](http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html) and includes 4 factor returns, the excess market, the size factor, the value factor and the momentum factor.  The available test portfolios include the 12 industry portfolios, a subset of the size-value two way sort, and a subset of the size-momentum two way sort. "
+    "The data used comes from Ken French's [website](http://mba.tuck.dartmouth.edu/pages/faculty/ken.french/data_library.html) and includes 4 factor returns, the excess market, the size factor, the value factor and the momentum factor.  The available test portfolios include the 12 industry portfolios, a subset of the size-value two way sort, and a subset of the size-momentum two way sort. "
    ]
   },
   {
diff --git a/examples/panel_data-formats.ipynb b/examples/panel_data-formats.ipynb
index d95d60314f..df68d07fc5 100644
--- a/examples/panel_data-formats.ipynb
+++ b/examples/panel_data-formats.ipynb
@@ -9,7 +9,7 @@
     "There are two primary methods to express data:\n",
     "\n",
     "  * MultiIndex DataFrames where the outer index is the entity and the inner is the time index.  This requires using pandas.\n",
-    "  * 3D structures were dimension 0 (outer) is variable, dimension 1 is time index and dimension 2 is the entity index.  It is also possible to use a 2D data structure with dimensions (t, n) which is treated as a 3D data structure having dimensions (1, t, n). These 3D data structures can be pandas, NumPy or xarray."
+    "  * 3D structures where dimension 0 (outer) is variable, dimension 1 is time index and dimension 2 is the entity index.  It is also possible to use a 2D data structure with dimensions (t, n) which is treated as a 3D data structure having dimensions (1, t, n). These 3D data structures can be pandas, NumPy or xarray."
    ]
   },
   {
@@ -59,7 +59,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The `MultiIndex` `DataFrame` can be used to initialized the model.  When only referencing a single series, the `MultiIndex` `Series` representation can be used."
+    "The `MultiIndex` `DataFrame` can be used to initialize the model.  When only referencing a single series, the `MultiIndex` `Series` representation can be used."
    ]
   },
   {
@@ -79,7 +79,7 @@
    "metadata": {},
    "source": [
     "## NumPy arrays\n",
-    "3D NumPy arrays can be used to hand panel data where the three axes are `0`, items, `1`, time and `2`, entity. NumPy arrays are not usually the best format for data since the results all use generic variable names. \n",
+    "3D NumPy arrays can be used to handle panel data where the three axes are `0`, items, `1`, time and `2`, entity. NumPy arrays are not usually the best format for data since the results all use generic variable names. \n",
     "\n",
     "Pandas dropped support for their Panel in 0.25."
    ]
@@ -99,7 +99,7 @@
     "ntime = mi_data.index.levels[1].shape[0]\n",
     "np_lscrap = np_lscrap.reshape((nentity, ntime)).T\n",
     "np_hrsemp = np_hrsemp.reshape((nentity, ntime)).T\n",
-    "np_hrsemp.shape = (1, ntime, nentity)"
+    "np_hrsemp = np_hrsemp.reshape((1, ntime, nentity))"
    ]
   },
   {
@@ -120,7 +120,7 @@
    "source": [
     "## xarray DataArrays\n",
     "\n",
-    "xarray is a relatively new entrant into the set of packages used for data structures.  The data structures provided by ``xarray`` are relevant in the context of panel models since pandas `Panel` is scheduled for removal in the futures, and so the only 3d data format that will remain viable is an `xarray` `DataArray`. `DataArray`s are similar to pandas `Panel` although `DataArrays` use some difference notation.  In principle it is possible to express the same information in a `DataArray` as one can in a `Panel`"
+    "xarray is a relatively new entrant into the set of packages used for data structures.  The data structures provided by ``xarray`` are relevant in the context of panel models since pandas `Panel` has been removed, and so the only 3d data format that remains viable is an `xarray` `DataArray`. `DataArray`s are similar to pandas `Panel` although `DataArrays` use some different notation.  In principle it is possible to express the same information in a `DataArray` as one can in a `Panel`"
    ]
   },
   {
@@ -212,6 +212,33 @@
     ").fit()\n",
     "print(res)"
    ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Aligning inputs\n",
+    "\n",
+    "The rows of `dependent` and `exog` are paired by position. An input that has labels, such as a pandas `Series` or `DataFrame`, must therefore use the same labels, in the same order, as `dependent`, and a `ValueError` is raised if it does not, since the values would be attached to the wrong observations. NumPy arrays do not have labels and are always paired by position.\n",
+    "\n",
+    "Here the regressor is sorted by year and then firm, while the dependent variable is sorted by firm and then year. `reindex` puts the rows of the regressor in the order of the dependent variable."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "hrsemp = mi_data.hrsemp.sort_index(level=\"year\")\n",
+    "try:\n",
+    "    PanelOLS(mi_data.lscrap, hrsemp, entity_effects=True)\n",
+    "except ValueError as exc:\n",
+    "    print(exc)\n",
+    "\n",
+    "hrsemp = hrsemp.reindex(mi_data.index)\n",
+    "print(PanelOLS(mi_data.lscrap, hrsemp, entity_effects=True).fit().params)"
+   ]
   }
  ],
  "metadata": {
diff --git a/examples/panel_examples.ipynb b/examples/panel_examples.ipynb
index 3f3f398dae..4d0afc689d 100644
--- a/examples/panel_examples.ipynb
+++ b/examples/panel_examples.ipynb
@@ -17,7 +17,7 @@
     "\n",
     "F. Vella and M. Verbeek (1998), \"Whose Wages Do Unions Raise? A Dynamic Model of Unionism and Wage Rate Determination for Young Men,\" _Journal of Applied Econometrics_ 13, 163-183.\n",
     "\n",
-    "The data set consists of wages and characteristics for men during the 1980s.  The entity identifier is ``nr`` and the time identified is ``year``.  This data is used extensively in Chapter 14 of _Introduction to Econometrics_ by Jeffrey Wooldridge.\n",
+    "The data set consists of wages and characteristics for men during the 1980s.  The entity identifier is ``nr`` and the time identifier is ``year``.  This data is used extensively in Chapter 14 of _Introduction to Econometrics_ by Jeffrey Wooldridge.\n",
     "\n",
     "Here a ``MultiIndex`` ``DataFrame`` is used to hold the data in a format that can be understood as a panel. Before setting the index, a year ``Categorical`` is created which facilitated making dummies."
    ]
@@ -45,7 +45,7 @@
    "metadata": {},
    "source": [
     "## Basic regression on panel data\n",
-    "``PooledOLS`` is just plain OLS that understands that various panel data structures. It is useful as a base model. Here the log wage is modeled using all of the variables and time dummies."
+    "``PooledOLS`` is just plain OLS that understands the various panel data structures. It is useful as a base model. Here the log wage is modeled using all of the variables and time dummies."
    ]
   },
   {
@@ -82,7 +82,7 @@
    "metadata": {},
    "source": [
     "### Random effects\n",
-    "The random effects model is virtually identical to the pooled OLS model except that is accounts for the structure of the model and so is more efficient. Random effects uses a quasi-demeaning strategy which subtracts the time average of the within entity values to account for the common shock."
+    "The random effects model is virtually identical to the pooled OLS model except that it accounts for the structure of the model and so is more efficient. Random effects uses a quasi-demeaning strategy which subtracts the time average of the within entity values to account for the common shock."
    ]
   },
   {
@@ -118,7 +118,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The coefficient $\\theta_i$ determines how much demeaning takes place.  When this value is 1, the RE model reduces to the pooled model since this occurs when there is no variance in the effects.  When panels are unbalanced it will vary across entities, but in this balanced panel all values are the same."
+    "The coefficient $\\theta_i$ determines how much demeaning takes place.  When this value is 0, the RE model reduces to the pooled model since this occurs when there is no variance in the effects, and when it is 1, the RE model reduces to the fixed effects (within) estimator.  When panels are unbalanced it will vary across entities, but in this balanced panel all values are the same."
    ]
   },
   {
@@ -135,7 +135,7 @@
    "metadata": {},
    "source": [
     "### The between estimator\n",
-    "The between estimator is an alternative, usually less efficient estimator, can can be used to estimate model parameters.  It is particular simple since it first computes the time averages of $y$ and $x$ and then runs a simple regression using these averages.\n",
+    "The between estimator is an alternative, usually less efficient, estimator that can be used to estimate model parameters.  It is particularly simple since it first computes the time averages of $y$ and $x$ and then runs a simple regression using these averages.\n",
     "\n",
     "The year dummies are dropped since the averaging removes differences due to the year.  ``expersq`` was also dropped since it is fairly co-linear with ``exper``. These results are broadly similar to the previous models."
    ]
@@ -160,7 +160,7 @@
    "metadata": {},
    "source": [
     "#### Other options\n",
-    "The ``fit`` method of ``BetweenOLS`` takes an optional input ``reweight`` which uses GLS-like weighting when panels are unbalanced.  In this dataset the panel is balanced and so using this option would have no effect. By default this is ``False`` and so the averages are directly used in OLS even when there are difference numbers of observations across entities."
+    "The ``fit`` method of ``BetweenOLS`` takes an optional input ``reweight`` which uses GLS-like weighting when panels are unbalanced.  In this dataset the panel is balanced and so using this option would have no effect. By default this is ``False`` and so the averages are directly used in OLS even when there are different numbers of observations across entities."
    ]
   },
   {
@@ -203,10 +203,10 @@
    "metadata": {},
    "source": [
     "### Time Effects\n",
-    "Time effect can be added using ``time_effects=True``.  Here the time dummies are removed. Note that the core coefficients are identical. The only change is in the test statistic for poolability since not the \"effects\" include both entity and time, whereas before only entity were included.\n",
+    "Time effect can be added using ``time_effects=True``.  Here the time dummies are removed. Note that the core coefficients are identical. The only change is in the test statistic for poolability since now the \"effects\" include both entity and time, whereas before only entity were included.\n",
     "\n",
     "#### Effects vs Dummies\n",
-    "For variable which can be consistently estimated, such as time effects in the usual large N, fixed T panel, it is not necessary to include these as effects.  They can instead be implemented as dummy variables.  "
+    "For variables which can be consistently estimated, such as time effects in the usual large N, fixed T panel, it is not necessary to include these as effects.  They can instead be implemented as dummy variables.  "
    ]
   },
   {
@@ -227,7 +227,7 @@
    "metadata": {},
    "source": [
     "#### Other Options\n",
-    "Some additional options are available when using ``PanelOLS`` and effects.  The ``use_lsdv`` can be used to force the model to be estimated using dummy variables.  This is not necessary and will be slower in most circumstances.  This options is primarily included for testing.  ``auto_df`` instructs ``PanelOLS`` to determine the degree of freedom adjustment to make.  Unlike most other estimators, the treatment of effects depends on the covariance estimator.  For example, when using a classic covariance estimator, effects count as lost degrees of freedom.  When using entity effects and clustering by entity, they do not.  Setting ``auto_df=False`` will allow the entities to be counted or not, depending on the value of ``count_effects``."
+    "Some additional options are available when using ``PanelOLS`` and effects.  The ``use_lsdv`` can be used to force the model to be estimated using dummy variables.  This is not necessary and will be slower in most circumstances.  This option is primarily included for testing.  ``auto_df`` instructs ``PanelOLS`` to determine the degree of freedom adjustment to make.  Unlike most other estimators, the treatment of effects depends on the covariance estimator.  For example, when using a classic covariance estimator, effects count as lost degrees of freedom.  When using entity effects and clustering by entity, they do not.  Setting ``auto_df=False`` will allow the entities to be counted or not, depending on the value of ``count_effects``."
    ]
   },
   {
@@ -236,13 +236,13 @@
    "source": [
     "### Other Effects\n",
     "\n",
-    "`PanelOLS` can be used to estimated models with up to 2-effects fixed effects. These can include any combination of\n",
+    "`PanelOLS` can be used to estimate models with up to 2-effects fixed effects. These can include any combination of\n",
     "\n",
     "* Entity effects\n",
     "* Time effects\n",
     "* Other effects\n",
     "\n",
-    "Other effects are specified using the `other_effects` input of `PanelOLS`. The input should have the same number of observations as the `engog` and can have 1 or two columns.  Below we reproduce the model fitted with time effects only using `other_effects` to set the effect."
+    "Other effects are specified using the `other_effects` input of `PanelOLS`. The input should have the same number of observations as the `exog` and can have 1 or two columns.  Below we reproduce the model fitted with entity and time effects using `other_effects` to set the time effect."
    ]
   },
   {
@@ -270,7 +270,7 @@
    "metadata": {},
    "source": [
     "### Using first differences\n",
-    "First differencing is an alternative to using fixed effects when there might be correlation.  When using first differences, time-invariant variables must be excluded.  Additionally, only one linear time-trending variable can be included since this will look like a constant.  This variable will soak up all time-trends in the data, and so interpretations of these variable can be challenging."
+    "First differencing is an alternative to using fixed effects when there might be correlation.  When using first differences, time-invariant variables must be excluded.  Additionally, only one linear time-trending variable can be included since this will look like a constant.  This variable will soak up all time-trends in the data, and so interpretations of these variables can be challenging."
    ]
   },
   {
@@ -307,6 +307,84 @@
     "print(compare({\"BE\": be_res, \"RE\": re_res, \"Pooled\": pooled_res}))"
    ]
   },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "## Specification tests\n",
+    "\n",
+    "The results of the panel models include two specification tests. The Hausman test examines whether the effects are correlated with the regressors, since random effects is only consistent if they are not. The Pesaran CD test examines whether the idiosyncratic shocks of different entities are correlated.\n",
+    "\n",
+    "### Hausman test\n",
+    "\n",
+    "The Hausman test compares the random effects estimates with those of a model that is consistent whether or not the effects are correlated with the regressors, usually a fixed effects model. If the effects are uncorrelated with the regressors, both estimators are consistent and random effects is more efficient. If they are correlated, only the fixed effects estimator is consistent. The test is a method of the random effects results, and its argument is the results of the other model. Both models must use the default unadjusted covariance estimator.\n",
+    "\n",
+    "The random effects model above includes the year dummies, and so it is compared with the fixed effects model that has entity and time effects, `fe_te_res`. The test compares the coefficients that are in both models, which are `expersq`, `married` and `union`. The time effects are not coefficients, and so only the entity effects are tested."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "re_res.hausman(fe_te_res)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "The p-value is close to 0, and so the null is rejected. There is evidence that the effects are correlated with the regressors, so that the random effects estimates are not consistent and the fixed effects estimates should be used. The random effects coefficients on `married` and `union` are larger than the fixed effects ones, which is consistent with unobserved characteristics that raise wages being more common among married men and union members.\n",
+    "\n",
+    "The result is a `HausmanTestResults` named tuple. Its fields are `statistic`, the test statistic with its p-value, and `estimates`, the coefficients that were compared and the standard error of their difference, and so `statistic, estimates = re_res.hausman(fe_te_res)` also works."
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "### Pesaran CD test\n",
+    "\n",
+    "The Pesaran CD test examines the null that the idiosyncratic shocks of different entities are uncorrelated. Its statistic is a scaled sum of the correlations over time of the shocks of each pair of entities, and it is asymptotically standard normal under the null. The test is the property `pesaran_cd` of the results, and it is computed from the idiosyncratic shocks, `idiosyncratic`.\n",
+    "\n",
+    "The test is applied here to fixed effects and random effects models with the same regressors and without year effects, which cannot be used with this test, as explained below."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "exog_vars = [\"expersq\", \"union\", \"married\"]\n",
+    "exog = sm.add_constant(data[exog_vars])\n",
+    "fe_cd_res = PanelOLS(data.lwage, exog, entity_effects=True).fit()\n",
+    "re_cd_res = RandomEffects(data.lwage, exog).fit()\n",
+    "print(\"Fixed effects\")\n",
+    "print(fe_cd_res.pesaran_cd)\n",
+    "print(\"\\nRandom effects\")\n",
+    "print(re_cd_res.pesaran_cd)"
+   ]
+  },
+  {
+   "cell_type": "markdown",
+   "metadata": {},
+   "source": [
+    "Both statistics are far in the tail of the standard normal distribution, and so independence is rejected: the wages of different men are correlated after the individual effects have been removed. Some of this dependence might be due to the common changes in wages over time, which these models do not include.\n",
+    "\n",
+    "`pesaran_cd` returns an `InapplicableTestStatistic` instead of a test statistic for models with time effects, since time effects demean the residuals across the entities that are observed in the same periods, which makes the residuals negatively correlated across entities and biases the statistic downwards. Year dummies that are included as regressors have the same effect on the residuals, which is why they are left out above. The test is also not applicable to `BetweenOLS` and `FamaMacBeth`. Printing the result shows why."
+   ]
+  },
+  {
+   "cell_type": "code",
+   "execution_count": null,
+   "metadata": {},
+   "outputs": [],
+   "source": [
+    "print(fe_te_res.pesaran_cd)"
+   ]
+  },
   {
    "cell_type": "markdown",
    "metadata": {},
@@ -319,7 +397,7 @@
    "metadata": {},
    "source": [
     "### Heteroskedasticity Robust Covariance\n",
-    "White\"s robust covariance can be used by setting ``cov_type=\"robust``. This estimator adds some robustness against certain types of specification issues but should not be used when using fixed effects (entity effects) since it is no longer robust.  Instead a clustered covariance is required."
+    "White's robust covariance can be used by setting ``cov_type=\"robust\"``. This estimator adds some robustness against certain types of specification issues but should not be used when using fixed effects (entity effects) since it is no longer robust.  Instead a clustered covariance is required."
    ]
   },
   {
@@ -339,7 +417,7 @@
    "metadata": {},
    "source": [
     "### Clustered by Entity\n",
-    "The usual variable to cluster are are entity or entity and time.  The can be implemented using ``cov_type=\"clustered\"`` and the additional keyword arguments ``cluster_entity=True`` and/or ``cluster_time=True``. "
+    "The usual variables to cluster on are entity or entity and time.  These can be implemented using ``cov_type=\"clustered\"`` and the additional keyword arguments ``cluster_entity=True`` and/or ``cluster_time=True``. "
    ]
   },
   {
@@ -373,9 +451,9 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "An ``OrderedDict`` is used to hold the results for comparing models.  This allows the models to be named as well as for the order of the models to be specified.  A standard ``dict`` will produce effectively random order.\n",
+    "An ``OrderedDict`` is used to hold the results for comparing models.  This allows the models to be named as well as for the order of the models to be specified.  A standard ``dict`` also preserves the order in which the models are added.\n",
     "\n",
-    "Clustering on entity reduced the t-stats across the board.  This suggests there is important correlation in the residuals per entity. Clustering by both also decreases the t-stats which suggests that there is cross-sectional dependence in the data. *Note*: clustering by entity addresses correlation across time and clustering by time controls for correlation between entities in a time period."
+    "Clustering on entity reduced the t-stats across the board.  This suggests there is important correlation in the residuals per entity. Clustering by both entity and time gives t-stats that are close to those from clustering by entity alone, some larger and some smaller, and so most of the reduction relative to the robust covariance comes from the correlation within an entity. *Note*: clustering by entity addresses correlation across time and clustering by time controls for correlation between entities in a time period."
    ]
   },
   {
diff --git a/examples/panel_using-formulas.ipynb b/examples/panel_using-formulas.ipynb
index 6e20bc0283..be3d8c1abd 100644
--- a/examples/panel_using-formulas.ipynb
+++ b/examples/panel_using-formulas.ipynb
@@ -6,7 +6,7 @@
    "source": [
     "# Using formulas to specify models\n",
     "\n",
-    "All of the models can be specified using formulas. The formulas used here utilize [formulaic](https://matthewwardrop.github.io/formulaic/) are similar to those in [statsmodels](http://www.statsmodels.org).  The basis formula syntax for a single variable regression would be\n",
+    "All of the models can be specified using formulas. The formulas used here utilize [formulaic](https://matthewwardrop.github.io/formulaic/) and are similar to those in [statsmodels](http://www.statsmodels.org).  The basic formula syntax for a single variable regression would be\n",
     "\n",
     "```\n",
     "y ~ 1 + x\n",
@@ -55,7 +55,7 @@
    "source": [
     "## PanelOLS with Entity Effects and a constant\n",
     "\n",
-    "The constant can be explicitly included using the `1 + ` notation.  When a constant is included in the model, and additional constraint is imposed that the number of the effects is 0.  This allows the constant to be identified using the grand mean of the dependent and the regressors."
+    "The constant can be explicitly included using the `1 + ` notation.  When a constant is included in the model, an additional constraint is imposed that the sum of the effects is 0.  This allows the constant to be identified using the grand mean of the dependent and the regressors."
    ]
   },
   {
diff --git a/examples/system_examples.ipynb b/examples/system_examples.ipynb
index 03bafd369b..e061762664 100644
--- a/examples/system_examples.ipynb
+++ b/examples/system_examples.ipynb
@@ -29,7 +29,7 @@
     "\n",
     "$$\\hat{\\beta}_{GLS} = \\left(X^\\prime \\Omega^{-1} X\\right)^{-1}X^\\prime  \\Omega^{-1} Y$$\n",
     "\n",
-    "where $\\Omega^{-1} = \\Sigma^{-1} \\otimes I_{T}$, $\\Sigma_{ij}$ is the covariance between $\\epsilon_i$ and $\\epsilon_j$ and $T$ is the number of observations. The GLS estimator is only beneficial when the regressors in different models differ  and when residuals are correlated. There GLS estimates are identical to the multivariate OLS estimates when all regressors are common."
+    "where $\\Omega^{-1} = \\Sigma^{-1} \\otimes I_{T}$, $\\Sigma_{ij}$ is the covariance between $\\epsilon_i$ and $\\epsilon_j$ and $T$ is the number of observations. The GLS estimator is only beneficial when the regressors in different models differ  and when residuals are correlated. The GLS estimates are identical to the multivariate OLS estimates when all regressors are common."
    ]
   },
   {
@@ -50,7 +50,7 @@
    "metadata": {},
    "source": [
     "## Data\n",
-    "Two data sets will be used.  The first is from Munnell which looks at the effect of capital on state GDP.  This example follows the example in Chapter 10 in recent editions of Greene\"s _Econometric Analysis_.\n",
+    "Two data sets will be used.  The first is from Munnell which looks at the effect of capital on state GDP.  This example follows the example in Chapter 10 in recent editions of Greene's _Econometric Analysis_.\n",
     "\n",
     "The data is state-level but the model is estimated in region.  The first step is to aggregate the data by region.  All capital measures are summed and the unemployment rate is averaged using weights proportional to the total employment in each state."
    ]
@@ -117,7 +117,7 @@
    "source": [
     "## Basic Usage\n",
     "\n",
-    "Seemingly Unrelated Models are fairly complex and each equation could have a different number of regressors.  As a result, it is not possibly to use standard `pandas` or `numpy` data structures, and so dictionaries (or technically dictionary-like objects) are used.  In practice, it is strongly recommended to use a `OrderedDictionary` from the `collections` module.  This ensures that equation order will be preserved. In addition, the dictionary must have the following structure:\n",
+    "Seemingly Unrelated Models are fairly complex and each equation could have a different number of regressors.  As a result, it is not possible to use standard `pandas` or `numpy` data structures, and so dictionaries (or technically dictionary-like objects) are used.  In practice, it is strongly recommended to use an `OrderedDict` from the `collections` module.  This ensures that equation order will be preserved. In addition, the dictionary must have the following structure:\n",
     "\n",
     "* `keys` **must be strings** and will be used as equation labels\n",
     "* The value associated with each key must be either a dictionary or a tuple.\n",
@@ -125,7 +125,7 @@
     "  * When a dictionary is used, it must have two keys, `dependent` and `exog`.  It can optionally have a third key `weights` which provides weights to use in the regression.\n",
     "  * When a tuple is used, it must have two elements and takes the form `(dependent, exog)`.  It can optionally contains weights in which case it takes the form `(dependent, exog, weights)`.\n",
     "\n",
-    "This example uses the dictionary syntax to contain the data for each region and uses the region identified as the equation label."
+    "This example uses the dictionary syntax to contain the data for each region and uses the region identifier as the equation label."
    ]
   },
   {
@@ -152,7 +152,7 @@
    "source": [
     "Fitting the model is virtually identical to  fitting any other model with the exception of the special data structure required. \n",
     "\n",
-    "The fitting options here ensure that the homoskedastic covariance estimator are used (`cov_type=\"unadjusted\"`) and that a small sample adjustment is applied. By default, GLS is used (this can be overridden using `method=\"ols\"`."
+    "The fitting options here ensure that the homoskedastic covariance estimator is used (`cov_type=\"unadjusted\"`). No small sample adjustment is applied, which is the default (`debiased=False`). By default, GLS is used (this can be overridden using `method=\"ols\"`)."
    ]
   },
   {
@@ -231,7 +231,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Individual results are contained in a dictionary located at the attribute `equations` and can be accessed using equation labels (available using the attribute `equation_labels`).  Additional information about the model is presented in this view. The West Coast results are show.  "
+    "Individual results are contained in a dictionary located at the attribute `equations` and can be accessed using equation labels (available using the attribute `equation_labels`).  Additional information about the model is presented in this view. The West Coast results are shown.  "
    ]
   },
   {
@@ -267,7 +267,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "These results can be compared to the results in Greene -- they are unsurprisingly identical."
+    "These results can be compared to the FGLS estimates in Table 10.1 of Greene, which is shown below. They agree with the table, up to differences in the last printed digit, except for three printed entries that look like misprints: the entries for `UNEMP` in SW and for `lnUTIL` in WC have a misplaced decimal point, and the entry for `lnPC` in NE has the opposite sign. The estimates also agree with an independent computation of the FGLS estimator."
    ]
   },
   {
@@ -285,7 +285,7 @@
    "source": [
     "The GLS estimation method requires stronger assumptions for parameter estimates to be consistent.  If these are violated then it might be the case that OLS is still consistent (in some sense) and so OLS can be used by passing `method=\"ols\"` when calling `fit`. \n",
     "\n",
-    "These results can be compared to Greene\"s table -- they are identical except the final value which seems to have a small typo."
+    "These results can be compared to the OLS estimates in Table 10.1, which is shown above. They agree with the table, up to differences in the last printed digit, except for five printed entries: `UNEMP` for CN, which is the final value, has a misplaced decimal point, `lnEMP` for MA has two digits transposed, and the intercepts for GF and SO and `lnUTIL` for SO differ by about 0.003 to 0.005."
    ]
   },
   {
@@ -311,7 +311,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The parameter estimates for one coefficient -- unemployment -- can be compared across the two estimation methods.  "
+    "The parameter estimates for one coefficient -- private capital (`lnPC`) -- can be compared across the two estimation methods.  "
    ]
   },
   {
@@ -338,7 +338,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The parameters and their standard errors match those reported in Greene."
+    "The parameters and their standard errors match those in Table 10.3 of Greene, which is shown below, except for the FGLS estimate for NE, which has the opposite sign in the printed table."
    ]
   },
   {
@@ -420,7 +420,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The estimator can be forced to use GLS by setting `method=\"gls\"`.  As can be seen below, the parameter estimates and standard errors do not change even though two-stage GLS is used.  The $R^2$ **does change** since the left-hand side variable is transformed by the GLS weighting before estimation. "
+    "The estimator can be forced to use GLS by setting `method=\"gls\"`.  As can be seen below, the parameter estimates, standard errors and $R^2$ do not change even though two-stage GLS is used. "
    ]
   },
   {
@@ -501,11 +501,11 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The estimator supports heteroskedasticity robust covariance estimation b setting `cov_type=\"robust\"`.  The estimator allows for arbitrary correlation across series with the same time index but not correlation across time periods, which is the same assumption as in the unadjusted estimator. The main difference is that this estimator will correct standard errors for dependence between regressors (or squared regressors) and squared residuals.\n",
+    "The estimator supports heteroskedasticity robust covariance estimation by setting `cov_type=\"robust\"`.  The estimator allows for arbitrary correlation across series with the same time index but not correlation across time periods, which is the same assumption as in the unadjusted estimator. The main difference is that this estimator will correct standard errors for dependence between regressors (or squared regressors) and squared residuals.\n",
     "\n",
     "#### Heteroskedasticity Robust Covariance Estimation\n",
     "\n",
-    "In the fringe benefit model there are some large differences between standard errors computed using the the homoskedastic covariance estimator and the heteroskedasticity robust covariance estimator (e.g., `exper`)."
+    "In the fringe benefit model there are some large differences between standard errors computed using the homoskedastic covariance estimator and the heteroskedasticity robust covariance estimator (e.g., `exper`)."
    ]
   },
   {
@@ -522,7 +522,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Other supported covariance estimators include `\"kernel`\" which implements a HAC and `\"clustered\"` which supports 1 and 2-way clustering. \n",
+    "Other supported covariance estimators include `\"kernel\"` which implements a HAC and `\"clustered\"` which supports 1 and 2-way clustering. \n",
     "\n",
     "#### Kernel (HAC)\n",
     "\n",
@@ -654,7 +654,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "The constraint can be added to an existing model using the `add_constraint` method.  This method requires one input, `r` and optionally `q`.  if `q` is not provided, it is set to 0. \n",
+    "The constraint can be added to an existing model using the `add_constraints` method.  This method requires one input, `r` and optionally `q`.  if `q` is not provided, it is set to 0. \n",
     "\n",
     "Here the constraint is added, the model is estimated, and the parameters for unemployment are displayed.  They all have the same value."
    ]
@@ -682,7 +682,7 @@
     "\n",
     "so that the restriction is identical to the previous, only applied to all variables excluding the constant.\n",
     "\n",
-    "Here the estimated from the first two equations are shown. All coefficients except the intercept are identical across equations.\n",
+    "Here the estimates from the first two equations are shown. All coefficients except the intercept are identical across equations.\n",
     "\n",
     "**Note**: When linear constraints are imposed, the parameter covariance matrix will be singular.  Caution is needed to ensure test statistics are meaningful."
    ]
@@ -747,7 +747,7 @@
    "source": [
     "### Using GLS with common regressors\n",
     "\n",
-    "As noted previously, forcing GLS has no effect (except on changing the $R^2$).  This can be seen below."
+    "As noted previously, forcing GLS has no effect.  This can be seen below."
    ]
   },
   {
diff --git a/examples/system_formulas.ipynb b/examples/system_formulas.ipynb
index 16f63e20e8..fd02c1cee0 100644
--- a/examples/system_formulas.ipynb
+++ b/examples/system_formulas.ipynb
@@ -12,7 +12,7 @@
     "* Single formula separated using {}\n",
     "\n",
     "These examples use data on fringe benefits from F. Vella (1993), \"A Simple Estimator for Simultaneous Models with Censored\n",
-    "Endogenous Regressors\" which appears in Wooldridge (2002).  The model consists of two equations, one for hourly wage and the other for hourly benefits.  The initial model uses the same regressors in both equations. "
+    "Endogenous Regressors\" which appears in Wooldridge (2002).  The model consists of two equations, one for hourly wage and the other for hourly benefits.  The two equations have different regressors. "
    ]
   },
   {
@@ -138,7 +138,7 @@
     "\n",
     "$$ V(y_i \\times w_i) = \\sigma^2 \\,\\,\\forall i.$$\n",
     "\n",
-    "Weights can be passed using a `DataFrame` where each column.  \n",
+    "Weights can be passed using a `DataFrame` with a column for each equation.  \n",
     "\n",
     "Here the results are printed to ensure that the estimates are different from those in the standard GLS model."
    ]
@@ -149,7 +149,8 @@
    "metadata": {},
    "outputs": [],
    "source": [
-    "random_weights = np.random.chisquare(5, size=(616, 2))\n",
+    "rng = np.random.default_rng(20261012)\n",
+    "random_weights = rng.chisquare(5, size=(616, 2))\n",
     "random_weights = pd.DataFrame(random_weights, columns=[\"benefits\", \"earnings\"])\n",
     "weighted_mod = SUR.from_formula(formula, data, weights=random_weights)\n",
     "print(weighted_mod.fit())"
diff --git a/examples/system_three-stage-ls.ipynb b/examples/system_three-stage-ls.ipynb
index 2b30ae18e4..8284686b88 100644
--- a/examples/system_three-stage-ls.ipynb
+++ b/examples/system_three-stage-ls.ipynb
@@ -12,11 +12,11 @@
     "\\begin{eqnarray}\n",
     "hours & = & \\beta_0 + \\beta_1 \\ln(wage) + \\beta_2 educ + \\beta_3 age + \\beta_4 kidslt6 + \\beta_5 nwifeinc + \\epsilon^h_i \n",
     "\\\\\n",
-    "\\ln(wage) & = & \\gamma_0 + \\gamma_1 hours + \\gamma_2 educ + \\gamma_3 educ^2 + \\gamma_4 exper + \\epsilon^w_i \n",
+    "\\ln(wage) & = & \\gamma_0 + \\gamma_1 hours + \\gamma_2 educ + \\gamma_3 exper + \\gamma_4 exper^2 + \\epsilon^w_i \n",
     "\\end{eqnarray}\n",
     "$$\n",
     "\n",
-    "Each equation has a single exogenous variables.  The instruments for the endogenous variables are the regressors that appear in one equation but not the other. "
+    "Each equation has a single endogenous variable.  The instruments for the endogenous variables are the regressors that appear in one equation but not the other. "
    ]
   },
   {
@@ -25,7 +25,7 @@
    "source": [
     "## Data\n",
     "\n",
-    "The data set is the MORZ data set from Wooldridge (2002)."
+    "The data set is the MROZ data set from Wooldridge (2002)."
    ]
   },
   {
@@ -66,7 +66,7 @@
     "\n",
     "* `IV2SLS` - single equation 2-stage least squares\n",
     "* `IV3SLS` - system estimation using instrumental variables\n",
-    "* `SUR` - system estimation without endogenous variables\n"
+    "* `IVSystemGMM` - system estimation using instrumental variables and GMM\n"
    ]
   },
   {
@@ -145,7 +145,7 @@
    "cell_type": "markdown",
    "metadata": {},
    "source": [
-    "Using `method=\"gls\"` will use GLS estimates which can be more efficient than the usual estimates. Here only the first equation changes.  This is due to the structure of the problem."
+    "Using `method=\"gls\"` will use GLS estimates which can be more efficient than the usual estimates. Here the estimates of both equations change."
    ]
   },
   {
@@ -261,7 +261,7 @@
     "$$\n",
     "\\begin{align}\n",
     "y_1 = y_2 + x_1 + \\epsilon_1 \\\\\n",
-    "y_2 = -y_1 + x_2 + \\epsilon_1 \\\\\n",
+    "y_2 = -y_1 + x_2 + \\epsilon_2 \\\\\n",
     "\\end{align}\n",
     "$$\n",
     "\n",
@@ -269,7 +269,7 @@
     "$$\n",
     "\\begin{align}\n",
     "y_1 = \\frac{1}{2}\\left(x_1 + x_2 + \\epsilon_1 + \\epsilon_2\\right) \\\\\n",
-    "y_1 = \\frac{1}{2}\\left(x_2 - x_1 + \\epsilon_2 - \\epsilon_1\\right) \\\\\n",
+    "y_2 = \\frac{1}{2}\\left(x_2 - x_1 + \\epsilon_2 - \\epsilon_1\\right) \\\\\n",
     "\\end{align}\n",
     "$$\n",
     "\n",
@@ -294,7 +294,7 @@
     "in_sample = df.iloc[:-10000]\n",
     "oos = df.iloc[-10000:]\n",
     "mod = IV3SLS.from_formula(\n",
-    "    {\"y1\": \"y1 ~  x1 + [y2 ~ x2]\", \"y2\": \"y2 ~ x2 + [y1 ~ x1]\"}, data=df\n",
+    "    {\"y1\": \"y1 ~  x1 + [y2 ~ x2]\", \"y2\": \"y2 ~ x2 + [y1 ~ x1]\"}, data=in_sample\n",
     ")\n",
     "res = mod.fit()\n",
     "print(res)"