EconBase
← Back to paper

Estimation with Pairwise Observations

Extracted main text — title through conclusion, appendix excluded. This is what our citation measures are computed over, published so the extraction can be checked by eye.

79,537 characters · 26 sections · 15 citation commands

Rendered from LaTeX for readability, not typeset faithfully. Citation keys are highlighted; maths is left as source; figures, tables and equation environments are summarised rather than reproduced; unrecognised commands are greyed out so nothing is silently dropped. Email addresses are removed.

Estimation with Pairwise Observations

\vskip 2cm \vskip 0.5cm

abstractThe paper introduces a new estimation method for the standard linear regression model. The procedure is not driven by the optimisation of any objective function rather, it is a simple weighted average of slopes from observation pairs. The paper shows that such estimator is consistent for carefully selected weights. Other properties, such as asymptotic distributions, have also been derived to facilitate valid statistical inference. Unlike traditional methods, such as Least Squares and Maximum Likelihood, among others, the estimated residual of this estimator is not by construction orthogonal to the explanatory variables of the model. This property allows a wide range of practical applications, such as the testing of endogeneity, i.e.,the correlation between the explanatory variables and the disturbance terms. \begin{keywords}{Linear regression model, consistent estimation, endogeneity, testing for endogeneity.} \end{keywords} \begin{JEL} {C01, C13, C20, C26, C51} \end{JEL}

Introduction

To start with, the paper intents to answer a simple question relating to the estimation of a linear regression model. Let us assume that we have a very basic model with only one explanatory variable:

equation[equation omitted — 123 chars of source]

Assuming the $n$ observations the data can be divided in clusters in such a way that in each cluster there are only two observations. This clustering may be carried out in two different ways:

itemize• The $(x_{i}, y_{i})$ observations are ordered based on the values of $x_{i}$ and two adjacent observations are grouped into one cluster. Let us call this the {\it sorted adjacent} approach. There are $n-1$ such pairs. When we do not apply any sorting and just group each adjacent observations into a cluster, we call such approach {\it non-sorted adjacent}. • All feasible combinations \begin{equation*} [(x_{i}, y_{i}) ; (x_{j}, y_{j})],\,\, (i \neq j) \end{equation*} of the observations are considered to form clusters. There are in total $\frac{n(n-1)}{2}$ of such clusters. Let us call this the (non-sorted) {\it full-pairwise approach}. Here, again we may consider sorting the observations based on the values of $x_{i}$, which gives the {\it sorted full-pairwise approach}.

In each cluster there are only two observations so the “regression line" in a cluster is the line going through these two observations, and there is no estimation involved. In such a way we can get as many $\beta_0$ and $\beta_1$ parameters as the number of clusters (see Figures 1 and 2 for example).

figure[figure omitted — 117 chars of source]
figure[figure omitted — 132 chars of source]

Our main question of interest is how these “cluster-wise" or let us call them {\it pairwise} parameters relate to the “true" parameters of the data generation process formulated in ((ref)). Let us denote the pairwise parameters in the $i$th cluster with ${\beta_{i,0}}$ and ${\beta_{i, 1}}$ and consider the following two estimators:

itemize• The weighted average of $\beta_{i,1}$ \begin{equation} \hat{\beta}_1 = \frac{\sum_{i=1}^{n} w_i \beta_{i,1}}{\sum_{i=1}^{n} w_i } \end{equation} with some $w_i$ weights. • Minimise a loss function such as \begin{equation} \hat{\beta}_1 = \underset{\beta_1}{\arg\min} \sum_{i = 1}^{n} (w_{i}(\beta_{i,1}-\hat{\beta}_1))^2. \end{equation}

Estimators for $\beta_0$ can also be defined similarly and they will be discussed with more details in Section (ref). The question now is how $\hat{\beta}_1$ relates to the true parameter, $\beta_1$ formulated in Equation ((ref)). Let us call this estimation procedure {\it Estimation with Pairwise Observations} (EwPO).

Some Monte Carlo Simulation Results

First, we adopted a Petri-dish type approach and carried out a large number of Monte Carlo simulation. The setup and the main results can be found in Appendix (ref).\footnote{Additional simulation results can be found in the Online Supplement EwPO-Sup2023.} The results provide some interesting insights in terms of the consistency of the proposed estimators under different choices of $w_i$. The EwPO turned out to be consistent in the following cases:

enumerate• EwPO as defined in Equation ((ref)) with \begin{enumerate} • Sorted full-pairwise and $w_{ij} = x_i - x_j$ • Non-sorted adjacent and $w_i = |x_i-x_{i-1}|$ • Sorted full-pairwise and $w_{ij} = |x_i-x_j|$ • Non-sorted full-pairwise and $w_{ij} = |x_i - x_j|$ • Sorted full-pairwise and $w_{ij} = \sqrt{(x_i-x_j)^2 + (y_i-y_j)^2}$ • Non-sorted full-pairwise and $w_{ij} = \sqrt{(x_i-x_j)^2 + (y_i-y_j)^2}$ \end{enumerate} • Loss function defined in Equation ((ref)) \begin{enumerate} • Non-sorted adjacent and $w_i = x_i - x_{i-1}$ • Sorted full-pairwise and $w_{ij} = x_i-x_j$ • Non-sorted full-pairwise and $w_{ij} = x_i-x_j$ \end{enumerate} Note: In case of the loss function taking the absolute value of weights does not matter as weights are squared in the functional form!

Some Analytical Results

Next, let us presents some analytical explanations for the simulations results. We also present the notations and definitions used throughout this paper.

Notations and Definitions

In order to develop the theoretical foundation, let us start with introducing a few notations, definitions and concepts. Let $\Delta x_{ij}$ denotes the difference between $x_i$ and $x_j$ that is,

equation[equation omitted — 88 chars of source]

In the event that the subscript is suppressed, then $j=i-1$. That is $\Delta x_i = x_i - x_{i-1}$.

The pairwise parameter is defined to be

equation[equation omitted — 94 chars of source]

In the event that $j=i-1$ we simplify the notation by writing $\beta_{1,i}$. That is $\beta_{1,i}=\beta_{1,(i,j)}$.

It is useful to also note that

equation[equation omitted — 95 chars of source]

Partial Sum

Some of the expressions for analysing the full-pairwise case contain several partial sums which are asymptotically related to a Brownian Motion or Wiener Process. So it is useful to introduce some additional notations and state the relevant results here. Define $\floor{x}$ as the largest integer that is less than or equal to $x$ and $\ceil{x}$ as the smallest integer that is larger than or equal to $x$. Let $r\in [0,1]$, $u_i \sim N(0,\sigma^2_u)$ and define

equation[equation omitted — 111 chars of source]

then under the conditional stated in billingsley:1999, the Functional Central Limit Theorem ensures $W_n(r) \overset{d}{\rightarrow} W(r)$, where $W(r)$ denotes a Wiener Process.

Adjacent Case

The estimator considered in this section is

equation[equation omitted — 149 chars of source]

where $w_i$ represents the weight for the $i^{th}$ pairwise parameter.

Case 1: \texorpdfstring{$w_i = \Delta x_i$}{weight equals to Delta x}

In this case, Equation ((ref)) reduces to

align*[align* omitted — 416 chars of source]

Now, in the case of adjacent pairwise parameter

align*[align* omitted — 129 chars of source]

The last line follows because the first term in every bracket is being cancelled out by the last term in the next bracket. The only two exceptions are $x_1$ and $x_n$. So this gives

equation[equation omitted — 100 chars of source]

Equation ((ref)) holds regardless whether the data was sorted or not. The only difference is that, in the sorted case $x_n = \max \{x_i\}^n_{i=1}$ and $x_1 = \min \{ x_i \}^n_{i=1}$. Equation ((ref)) does seem to suggest that $\hat{\beta}_1$ is not going to be consistent in this case, as $\displaystyle = \frac{\Delta u_{n1}}{\Delta x_{n1}}$ remains even when $n\rightarrow \infty$.

Case 2: \texorpdfstring{$w_i = |\Delta x_i|$}{Weight equals to absolute Delta x}

In this case, it would be useful to first define the $\sgn{x}$ function.

equation[equation omitted — 134 chars of source]

Note that $\sgn{x}x = |x|$, and that unlike the case when $w_i = \Delta x_i$,

equation*[equation* omitted — 95 chars of source]

and no simple cancellation occurs here when the observations are \textsl{not sorted}. However, when the observations are \textsl{sorted}, it is clear that $\sum^n_{i=2} |\Delta x_i| = \sum^n_{i=2} \Delta x_i$ since $\Delta x_i = |\Delta x_i|$. So the result for the sorted case is the same as the case when $w_i = \Delta x_i$.

For the non-sorted case,

align*[align* omitted — 678 chars of source]

Under the assumption that $\E \left (|\Delta x_i| \right ) \leq \infty$, it is clear that $\E \left ( |\Delta x_i| \right ) > 0$. Note that $\E (|\sgn{\Delta x_i}u_i|) = \E ( |u_i| )$, so under that assumption that $\E (u_i) = 0$, $\E (|u_i|) < \infty$ and $u_i \perp x_i$ then $\E \left [ \sgn{\Delta x_i}u_i \right ] = \E \left [ \sgn{\Delta x_i} \right ]\E [ u_i ] = 0$. Moreover, under Weak Law of Large Number (WLLN),

equation*[equation* omitted — 71 chars of source]

and hence

equation*[equation* omitted — 55 chars of source]

Full-pairwise

An important difference here is that the number of pairwise parameters is much larger. Let $N$ denotes the number of full-pairwise parameters then $N = \displaystyle \frac{n(n+1)}{2}$. It is also important to identify the exact pair from each pairwise parameter and redefine the estimator as

equation[equation omitted — 187 chars of source]

where $\beta_{1,(i,j)}$ is defined as in Equation ((ref)).

Case 1: \texorpdfstring{$w_{ij} = \Delta x_{ij}$}{Weight equals to Delta x}

We first deal with the non-sorted case. In this case, the estimator reduces to

align[align omitted — 447 chars of source]

Let us first focus on the denominator as the argument for the numerator will be similar.

align[align omitted — 283 chars of source]

For reason that will become clear later, we are going to assume $x_i\sim D(\mu_x, \sigma^2_x)$ with $\sigma^2<\infty$ and $\E |x_i|<\infty$. We analyse expression ((ref)) term-by-term. The first term in expression ((ref)) can be rewritten as

align*[align* omitted — 201 chars of source]

where $z_i = x_i - \mu_x$. Note that $z_i$ is the de-meaned version of $x_i$ and therefore $\E(z_i)=0$. To expand further the expression above,

align[align omitted — 615 chars of source]

Now for the second term in expression ((ref)),

align[align omitted — 373 chars of source]

Substitute expressions ((ref)) and ((ref)) to expression ((ref)) yields

equation[equation omitted — 178 chars of source]

Note that $N=O(n^2)$ so

equation*[equation* omitted — 81 chars of source]

The same argument applies to $\sum^n_{i=1} \sum^{i-1}_{j=1} \Delta u_{ij}$ and since both approach 0 at the same rate, $\hat{\beta}_1$ is going to be unstable as $n\rightarrow \infty$.

The case when the data is sorted means that $\sum^n_{i=1} \sum^{i-1}_{j=1} \Delta x_{ij}$ is not going to be 0 asymptotically since $\Delta x_{ij} > 0$ by construction. Equation ((ref)) stills hold when $x_i$ is sorted such that $x_i \geq x_j$ for all $i>j$. Since $u_i$ is independent from $x_i$ this means the numerator in Equation ((ref)) approaches 0 as $n\rightarrow \infty$ but the denominator either approaches a finite non-zero positive constant or diverges. In both cases, the $\hat{\beta}_1 - \beta_1 = o_p(1)$ under the Assumptions of Proposition (ref), which provides the asymptotic distribution for this case.

Case 2: \texorpdfstring{$w_{ij} = |\Delta x_{ij}|$}{Weight equals to Absolute Delta x}

Following similar derivation as above, it is straightforward to show that

align[align omitted — 1,011 chars of source]

Given that $x_i\perp u_j$ for all $i,j=1,\ldots, n$ and $n^{-1}\sum^n_{i=1} |\sgn{\Delta x_{ij}}| \leq 1$, the last two terms in Equation ((ref)) therefore converges to 0 following a similar argument as in the case $w_i = |\Delta x_{ij}|$. Using the same argument as above, $\sum^n_{i=1}\sum^{i-1}_{j=1} |\Delta x_{ij} | = O_p(n^2)$ with a non-zero positive bound. Thus $\hat{\beta}_1 - \beta_1 = o_p(1)$.

An interesting observation from the theoretical discussion above is that the intercept, $\beta_0$, does not play a role in the consistency of EwPO estimator. The assumption of exogeneity, however, remains essential for the identification and consistency of $\beta_1$ as well as $\beta_0$, as shown in Section (ref) below.

Confidence Intervals

The asymptotic distribution of $\hat{\beta}$ in the full-pairwise case with $w_{ij} = \Delta x_{ij}$ can be found in Proposition (ref). Analytical results for the asymptotic behaviour of the full-pairwise case with $w_{ij} = |\Delta x_{ij}|$ is more difficult as the presence of absolute value creates some technical challenges. In practice, however, hypothesis testing in this case can still be conducted via the jackknife procedure.

Regardless on the choice of $w_{ij}$, the distribution of $\hat{\beta}$ in the full-pairwise case is likely to be non-standard and involve stochastic integrals. Therefore, the critical values is still required to be simulated even if the analytical distributions can be obtained. From a practical viewpoint, the jackknife procedure does not necessarily induce more computational cost and it is generally much simpler to conduct as the procedure is the same across different choices of $w_{ij}$.

An example of the procedure to estimate critical values can be found as follows:

enumerate[{Step} 1.] • Set the level of significant for the test $\alpha$. • Set the number of observations to be removed from the main sample, call it $d$ where $\sqrt{n} < d < n$. • Set the number of replication $R < \begin{pmatrix} n \\ n-d \end{pmatrix}$ • set $i=1$ • Create a jackknife sample, $J_i$, by removing $d$ randomly selected observations from the main sample. • Compute $\hat{\beta}$ using the jackknife sample $J_i$, call it $\hat{\beta}_i$ • Set $i\leftarrow i+1$. • Repeat Steps (ref) to (ref) until $i=R$. • Sort the set $\{\hat{\beta}_i\}_{i=1}^R$ from lowest to highest. • The $\displaystyle \floor{\frac{\alpha}{2} R}^{th}$ and $\floor{\frac{(1-\alpha)}{2}R}^{th}$ elements in the set gives the lower and upper bounds of the confidence interval at the $\alpha$ significant level.

Table (ref) provides an example on the performance of the jackknife algorithm. As shown in the table, the algorithm performs well in obtaining the confidence intervals of the full-pairwise estimator.

table[table omitted — 1,742 chars of source]

What about the Intercept?

The discussion above focuses on the slope. This section investigates the relation between pairwise parameter and the intercept in Equation ((ref)). The most obvious approach is to consider

equation[equation omitted — 98 chars of source]

where $\hat{\beta}_1$ is any consistent estimator of $\beta_1$ with $\bar{y}_n$ and $\bar{x}_n$ denoting the averages of $y_i$ and $x_i$ over a sample of $n$ observations, respectively. That is, $\bar{y}_n = n^{-1} \sum^n_{i=1} y_i$ and similarly for $\bar{x}_n$.

It is straightforward to show that $\hat{\beta}_0$ is consistent since

align*[align* omitted — 140 chars of source]

and as $n\rightarrow \infty$, $\beta_1 - \hat{\beta}_1 = o_p(1)$, $\bar{x}_n - \E(x_i) = o_p(1)$ an $\bar{u}_n = o_p(1)$. Thus, by the Continuous Mapping Theorem, $\hat{\beta}_0 - \beta_0 = o_p(1)$.

While this provides a consistent estimator, we would also like to examine if there is a connection between the pairwise parameter of $\beta_0$ and $\beta_0$. Define the pairwise parameter of the intercept as

equation[equation omitted — 89 chars of source]

The adjacent case

In the adjacent case we consider the following estimator of $\beta_0$

equation[equation omitted — 123 chars of source]

Since with $w_i= |\Delta x_i|$ the non-sorted observations is the only case that produces consistent estimate of $\beta_1$, we will focus on this particular case for $\beta_0$.

align[align omitted — 388 chars of source]

Expanding the following term further gives

align*[align* omitted — 201 chars of source]

Substitute the last expression above into Equation ((ref)) yields

equation[equation omitted — 190 chars of source]

Under the assumption that $x_i \perp u_i$ and the usual regularity assumptions on $x_i$, then by the WLLN, $\sum^n_{i=2} |\Delta x_i| u_i = o_p(1)$ and $\sum^n_{i=2} \sgn{\Delta x_i} x_i \Delta u_i$. The former holds, because the sum is converging to $\E \left ( |\Delta x_i| u_i \right ) = \E \left (|\Delta x_i| \right ) \E (u_i) = 0$ and the second expression holds because the sum is converging to $\E \left [ \sgn{\Delta x_i} x_i \Delta u_i \right ] = \E \left [\sgn{\Delta x_i} x_i\right ] \E \left (\Delta u_i \right )=0$.

Loss Function

Next, let us analyse the properties of estimator ((ref)):

equation[equation omitted — 157 chars of source]

where

equation*[equation* omitted — 61 chars of source]

\texorpdfstring{$w_i = |\Delta x_i|$}{Weight equals to Absolute Delta x}

Given the quadratic loss, here show that the case when $w_i = \Delta x_i$ is identical to the case when $w_i = |\Delta x_i|$.

The solution to Equation ((ref)) implies

equation*[equation* omitted — 83 chars of source]

and therefore

equation*[equation* omitted — 118 chars of source]

Thus when $w_i = |\Delta x_i|$ or when $w_i = \Delta x_i$ we get

equation*[equation* omitted — 138 chars of source]

which is in fact the OLS estimator for the `first difference' equation

equation[equation omitted — 66 chars of source]

Full-pairwise

Reconsider the estimator in the full-pairwise case

equation[equation omitted — 155 chars of source]

which gives

equation*[equation* omitted — 168 chars of source]

Now, if we set $w_{ij}=|\Delta x_{ij}|$ then

equation[equation omitted — 215 chars of source]

which under the assumption that $x_i\perp u_i$, shows that the estimator is consistent. The same applies to $w_{ij} = \Delta x_{ij}$. The equivalence is mostly due to the quadratic loss. An interesting case is $w_{ij} = \sqrt{|\Delta x_{ij} |}$, then

equation*[equation* omitted — 187 chars of source]

which is in fact our full-pairwise estimator.

Multivariate Extension

One way to extend the approach to multiple explanatory variables is to utilise the residual matrix. Consider the following Data Generating Process

equation[equation omitted — 142 chars of source]

with the matrix representation

equation[equation omitted — 71 chars of source]

where $\bfy = \left (y_1, \ldots, y_n\right )'$ is a $n\times 1$ vector containing the dependent variable, $\bfX = \left (\bfx_1, \ldots, \bfx_K \right )$ is a $n\times K$ matrix containing the $K$ explanatory variables such that $\bfX$ is $\bfx_k = \left ( x_{k1}, \ldots, x_{kn} \right )'$ for $k=1,\ldots, K$ and $\bfu = \left (u_1, \ldots, u_K \right )'$ is a $n\times 1$ vector containing the residuals.

It is convenient to introduce some special matrices here. These matrices allow us to represent the pairwise parameter as a sequence of matrix operations. Define $\bfd_{i,j}$ as a $1\times n$ vector such that the $i^{th}$ and $j^{th}$ elements are 1 and -1, respectively while all other elements are 0. Let $\bfD_F = \left ( \bfd_{n,n-1}', \bfd_{n,n-2}', \ldots, \bfd_{n,1}', \bfd_{n-1,n-2}', \ldots, \bfd_{n-1,1}', \ldots, \bfd_{2,1}' \right )'$ and $\bfD_a = \left ( \bfd_{2,1}', \bfd_{3,2}', \bfd_{4,3}'\ldots, \bfd_{n,n-1}'\right )'$. Given these matrices, it is straightforward to show that $\bfD_F \bfx = \{ \Delta x_{ij} \}_{i=2,j=1}^{n,i-1}$ is a $n(n-1)/2 \times 1$ vector and $\bfD_a\bfx = \{ \Delta x_i \}_{i=2}^{n}$ is a $(n-1) \times 1$ vector. The former gives the vector of Full-pairwise Difference while the latter yields the Adjacent Difference.

Let $\bfS:\R^n\rightarrow \R^n$ be a function that returns a $n\times 1$ zero-one selection matrix such that if $\bfz = \bfS(\bfx)\bfx$ then $\bfz$ is a $n\times 1$ vector containing the same elements as $\bfx$ such that $z_i \geq z_j$ for $j>i$. Thus, $\bfS(\bfx)\bfx$ produces the sorted version of $\bfx$ in descending order. Finally, define the residual maker

equation[equation omitted — 106 chars of source]

where $\bfX_{-k}$ denote the $n\times K-1$ matrix which contains all the columns of $\bfX$ except the $k^{th}$ column.

The first step is to develop the matrix representation of the univariate case, then the multivariate case can be derived by repeat applications of the residual maker. That is, we construct the multivariate version by combining the univariate estimator for each of the regressors.

equation*[equation* omitted — 65 chars of source]

and note that $\bfD \bfy = \bfD\bfx \beta_1 + \bfD\bfu$ for $\bfD = \{\bfD_a, \bfD_F\}$. The pairwise parameters for the sample can be written as

equation[equation omitted — 90 chars of source]

where $\bfdiag{\bfx}= {\left [ \bfI \otimes \left (\bfD \bfx \right )' \right ]\bfS_n }$ is a diagonal matrix with the elements of $\bfx$ in the diagonal. $\bfS_n$ is a $n^2\times n$ zero-one matrix such that all the $( (i-1)n+1, i)$ elements are 1 for $i=1,\ldots, n$ and all the other elements are 0. Note that Equation ((ref)) is reasonably general and covers all the previously discussed case by defining the transformation matrix $\bfD$ differently. A summary can be found in Table (ref).

table[table omitted — 393 chars of source]

Given the vector of pairwise parameters $\bfbeta_1$, the estimator of $\beta_1$ can be written as

equation[equation omitted — 121 chars of source]

where $\bfw_1$ is the vector of weights associated with each $\beta_{1,(i,j)}$ in $\bfbeta_1$. Following from the analysis above, some of the choices include $\bfw_1 = \left ( |\Delta x_1|, \ldots, |\Delta x_n| \right )'$ and $\bfw_1 = \left ( \Delta x_1, \ldots, \Delta x_n \right )'$. It is straightforward to generalise Equation ((ref)), and subsequently Equation ((ref)), to the case when $\bfX$ is a $n\times K$ matrix with $K$ regressors. Consider

align[align omitted — 177 chars of source]

where $\bfbeta_{-k}$ contains the same elements as $\bfbeta$ with the $k^{th}$ element removed. Following Equation ((ref)) the pairwise parameter of $\beta_k$ is

equation[equation omitted — 105 chars of source]

Similarly, Equation ((ref)) can be generalised as

equation[equation omitted — 125 chars of source]

and

equation[equation omitted — 526 chars of source]

As shown in previous results, suitable choice of $\bfw_k$ leads to consistency of $\hat{\beta}_k$ and therefore $\hat{\bfbeta} - \bfbeta = o_p(1)$ under the same conditions as those discussed above for all $k=1,\ldots, K$.

Testing for Endogeneity

There are multiple endogeneity tests available for empirical use, but all rely on some kind of additional or external information, mostly in the form of instrumental variables (see e.g., Hausmann1978, or Wooldridge2002, pp. 118-122).

The aim in this section is to develop two tests for endogeneity based on the estimator of the form

equation*[equation* omitted — 159 chars of source]

which solely relies on data already in use for the estimation.

Residuals Test

The residuals test is particularly useful when $\beta_0=0$. Consider $w_{ij} = |\Delta x_{ij}|$ which implies

equation[equation omitted — 76 chars of source]

where

equation[equation omitted — 196 chars of source]

and under $x_i \perp u_i$, $\delta_n = o_p(1)$. Now given the following specification

equation[equation omitted — 66 chars of source]

the estimated residual from the EwPO estimator is

align*[align* omitted — 124 chars of source]

It is straightforward to show that

align*[align* omitted — 161 chars of source]

Now as $n\rightarrow \infty$ the last line above is

equation[equation omitted — 92 chars of source]

under the assumption that $\E(u_i) = 0$. There are two cases when $n^{-1}\sum^n_{i=1} \hat{u}_i$ is $o_p(1)$ namely, $x_i \perp u_i$ or $\mu_x = 0$. Assuming $\mu_x\neq 0$, which can be easily verified in practice, it is possible to directly test $x_i \perp u_i$ by testing $H_0: \E(u_i) = 0$. This provides the foundation for testing endogeneity by examining the mean of the estimated residuals using standard testing procedure, such as the $t$-test.\footnote{The argument here also applies to the usual Ordinary Least Squares (OLS) estimator. That is, when $\beta_0=0$, the estimated residuals do not have 0 mean. Thus, it also provides a test of endogeneity in this special case.}

Covariance Test

Of course $\beta_0=0$ is often too restrictive in practice. One way to alleviate this is to remove the intercept by considering the model in difference. Consider $\Delta \hat{u}_{pq} = \Delta y_{pq} - \Delta x_{pq} \hat{\beta}$ i.e., the `estimated' residuals in the form

equation[equation omitted — 205 chars of source]

In the case of a consistent estimator $\hat{\beta}$ can be expressed in the form of $\hat{\beta} = \beta + \delta(\bfw, \bfu)$ where $\bfw$ and $\bfu$ denote the vectors of the weights i.e., $w_{ij}$ and the residuals $u_i$, respectively, for $i,j=1,\ldots, n$ and $B$ vanishes to 0 asymptotically under ideal conditions, Therefore

equation*[equation* omitted — 93 chars of source]

Hence, one can test $H_0:\E [\delta (\bfw, \bfu)] = 0$ by considering the test statistics

equation[equation omitted — 165 chars of source]

Case 1: \texorpdfstring{$\bfw = \{\Delta x_{ij} \}_{i,j=1}^n$}{Delta x }

In this case,

equation*[equation* omitted — 201 chars of source]

and the test statistics

equation[equation omitted — 363 chars of source]

Under the null that $x_i \perp u_i$, $S(\Delta x_{ij})$ has mean 0. The asymptotic behaviour of EwPO estimator as defined in Equation ((ref)) with $w_{ij} = \Delta x_{ij}$ and that of the test statistics as defined above can be found in Propositions (ref) and (ref), respectively. Before presenting the theoretical results, consider the following assumptions

enumerate[{Assumption} 1.] • $\E (x_i) = \mu_x <\infty$ $\forall i=1,\ldots, n$ and $\exists \sigma^2_x < \infty$ such that $n^{-1} \sum^n_{i=1} (x_i-\mu_x)^2 - \sigma_x^2 = o_p(1)$. • $\E (u_i) = 0$ $\forall i=1,\ldots, n$ and $\exists \sigma^2_u < \infty$ such that $n^{-1} \sum^n_{i=1} u^2_i - \sigma^2_u = o_p(1)$. • $\E (u_i| x_j) =0$ for $i,j=1,\ldots, n$. • $\Delta x_{ij}$ and $\Delta u_{ij}$ are mixingales as defined in mcleish_invariance_1975 and satisfy all the conditions for Theorem 3.8 in mcleish_invariance_1975.

Assumptions (ref) and (ref) are standard in the sense that the existence of second moments are required for both the regressor and the residuals. The exogeneity assumption as presented in Assumption (ref) is required for the consistency of the estimator. It also stipulates the null underlying the distribution of the test statistics as asserted in Proposition (ref). Assumption (ref) is required to ensure that the various partial sums converge to the Brownian Motion process using the Functional Central Limit Theorem. It is required due to the partial sums from the full-pair construction of the estimator. See also wooldridge_invariance_1988 for potentially more general assumptions on $x_i$ and $u_i$.

The asymptotic distribution of the full-pairwise EwPO with $w_{ij} = \Delta x_{ij}$ is presented in Proposition (ref) below

propositionUnder Assumptions (ref) to (ref) with $\hat{\beta}$ as defined in Equation ({(ref)}) and $\bfw = \{\Delta x_{ij} \}^n_{i,j=1}$: \begin{equation} \hat{\beta} - \beta \overset{d}{\rightarrow} \frac{ W(1) - 2\int^1_0 W(\lambda) d\lambda}{B(1) - 2\int^1_0 B(\lambda) d\lambda} . \end{equation}
proofSee Appendix (ref).

The proposition below gives the asymptotic distribution of the covariance test statistics generated by the estimator considered in Proposition (ref) above.

propositionUnder the assumptions of Proposition (ref) \begin{equation} \begin{split} (\sigma_x\sigma_u)^{-1} S(\bfw) \overset{d}{\rightarrow}& B(1)W(1) + \int^1_0 B(\lambda) W(\lambda) d\lambda - \int^1_0B(\lambda)dW(\lambda) \&\quad - \int^1_0W(\lambda) dB(\lambda) - B^2(1) \frac{ W(1) - 2\int^1_0 W(\lambda) d\lambda}{B(1) - 2\int^1_0 B(\lambda) d\lambda}, \end{split} \end{equation} where $B(\lambda)$ and $W(\lambda)$ are independent standard Brownian motion processes.
proofSee Appendix (ref)
remarkIt can be shown that both terms on the right hand side of Equation ((ref)) are negative on expectation. Under the alternative that $u_i$ and $x_i$ are correlated, the values of the last two stochastic integrals in Equation ((ref)) increase as the correlation increases.
remarkThe asymptotic distribution as stated in Equation ((ref)) is obviously non-standard and must be simulated in order to obtain the critical value(s) for the purpose of inference. Since $B(\lambda)$ and $W(\lambda)$ are independent Brownian processes, each of the stochastic integrals can be simulated as \begin{align} \int^1_0 B(\lambda) d\lambda \approx& \sum^n_{p=1} z_{1p} \\ \int^1_0 W(\lambda) d\lambda \approx& \sum^n_{p=1} z_{3p} \\ \end{align} where $z_{ip} = z_{ip-1} + \varepsilon_{ip}$, $\bfepsilon_p = (\varepsilon_{1p}, \varepsilon_{2p})' \sim NID(0, \bfI)$ and \begin{align} \int^1_0 W(\lambda) dB(\lambda) \approx & \sum^n_{p=1} z_{2p}\varepsilon_{1p} \\ \int^1_0 B(\lambda) dW(\lambda) \approx & \sum^n_{p=1} z_{1p}\varepsilon_{2p} . \end{align} Under the null, other quantities in the distribution namely, $\sigma_x$ and $\sigma_u$ can be estimated consistently from the data and estimated residuals, respectively. Thus the asymptotic distribution can be simulated by repeated computation of Equations ((ref)) - ((ref)) and the substitution their values into Equation ((ref)) for desired number of replications. For further details on the simulation of stochastic integrals, see for example, johansen:1995. By way of demonstration, the critical values from the simulated sample of Equation ((ref)) using Equations ((ref)) - ((ref)) with $\sigma_x=\sigma_u = 1$ and $n=10000$ can be found in Table (ref) and the respective histogram in Figure (ref). Note that this is simulated with $w_{ij} = \Delta x_{ij}$ and may not be appropriate when $w_{ij} = |\Delta x_{ij}|$. In this case, for example, the jackknife procedure may provide more accurate approximation.

{\color{blue}

figure[figure omitted — 157 chars of source]
table[table omitted — 359 chars of source]

}

Case 2: \texorpdfstring{$w_{ij} = |\Delta x_{ij}|$}{Absolute Delta x}

Similar derivation as above yields

equation[equation omitted — 385 chars of source]

The asymptotic distribution of EwPO with $w_{ij}=|\Delta x_{ij}|$ and the associated test statistics are presented in Propositions (ref) and (ref).

propositionUnder Assumptions (ref) to (ref) with $\hat{\beta}$ as defined in Equation ({(ref)}) and $\bfw = \{|\Delta x_{ij}| \}^n_{i,j=1}$: \begin{equation} \sqrt{n} \left (\hat{\beta} - \beta \right ) \overset{d}{\rightarrow} 2\frac{\sigma_u}{\mu_{\Delta x_s}} \left ( W(1) - 2\int^1_0 W(\lambda) d\lambda \right ) . \end{equation} where $\mu_{\Delta x_s} = \E (\Delta x_s)$ where $x_s$ denotes the sorted version of the regressor $x$ i.e., $\Delta x_{s,ij} = x_i - x_j$ where $x_i > x_j$.
proofSee Appendix (ref).
propositionUnder Assumptions (ref) to (ref) with $\hat{\beta}$ as defined in Equation ({(ref)}) and $\bfw = \{|\Delta x_{ij}| \}^n_{i,j=1}$: \begin{equation} \begin{split} (\sigma_x\sigma_u)^{-1} S(\bfw) \overset{d}{\rightarrow} B(1)W(1) +& \int^1_0 B(\lambda)W(\lambda) d\lambda\\ & - \int^1_0 B(\lambda)dW(\lambda) - \int^1_0 B(\lambda) dW(\lambda). \end{split} \end{equation}
proofSee Appendix (ref).

Appendix (ref) contains some Monte Carlo simulations of the test statistics with the two different weighting schemes under the null and selected alternatives. As shown in Appendix (ref), the mean test statistics deviates from zeros as the correlation between the regressor and the error term grows. Further investigations suggest that the covariance test has the right size with relatively good power, even in small sample, so could really be useful in practice.

Implications

Basically, the results can be used in two ways. First, we could calculate the appropriate critical values for a given data set and perform the testing procedure formally. However, this requires the simulation of the critical values which is not always convenient in practice.

Another way, however, to look at this is based on the fact that the value of the test statistic is a monotonically increasing function (in absolute value) of the correlation between $x_i$ and $u_i$, i.e., the degree of endogeneity. Therefore, the following selection procedure can be performed to reduce or minimise the risk of estimation bias due to endogeneity:

enumerate[a)] • Estimate the model, e.g., $y_i = \beta x_i +u_i$ with EwPO; • Estimate the model with EwPO after being transformed by the potential instrumental variable (IV) candidate, $z_i = \beta w_i + v_i$ where $z_i = g_iy_i$, $w_i = g_ix_i$ and $v_i = g_iu_i$ with $g_i$ an IV candidate that satisfies the same assumptions as $x_i$; • Select the model for inference with the lowest statistics value. This minimises the risk of endogeneity bias.

Note that the procedure above does not require the simulation of any distribution but a repeat application of the EwPO estimator. Thus, it is potentially a simpler but effective procedure.

Hausman Test

A third test is inspired by the well known Hausman Test proposed in Hausmann1978. The general idea is that the difference between OLS and EwPO with a suitably chosen weight should be zero asymptotically under the assumption of exogeneity. While EwPO may not be consistent in the presence of endogeneity, its asymptotic bias is generally not the same as that of the OLS depending on the choice of weight. This means their difference is non-zero under the alternative of endogeneity. This can be easily demonstrated by comparing the OLS estimator with an EwPO estimator. Recall the OLS estimator for a simple linear regression model $y_i = \beta_0 + \beta_1 x_i + u_i$ can be expressed as

align*[align* omitted — 143 chars of source]

Let $\delta_n(w_{ij}) = \left ( \sum^n_{i=2}\sum^{i-1}_{j=1} w_{ij} \right )^{-1} \left ( \sum^n_{i=2}\sum^{i-1}_{j=1} w_{ij} u_{ij} \right )$ then

equation*[equation* omitted — 100 chars of source]

where $\hat{\beta}_{1(i,j)} (w_{ij})$ denotes an EwPO estimator with the weight $w_{ij}$. It is clear that under the null of exogeneity, $\epsilon_n=o_p(1)$ and $\delta_n(w_{ij}) = o_p(1)$ with the appropriate choice of $w_{ij}$. Under the alternative that $u_i$ and $x_i$ are correlated however, $|\epsilon_n - \delta_n(w_{ij})| = C + o_p(1)$ where $C$ is a positive (non-zero) real number. Note that if EwPO is consistent under the null and the alternative then $\delta_n(w_{ij}) = o_p(1)$ but $\epsilon_n$ does not vanish asymptotically. This was the basis of the Hausman test, but $\delta_n(w_{ij}) = o_p(1)$ is not required for the test to be valid since the test statistics still has non-zero expectation under the alternative if the asymptotic bias are different between the two estimators. This means as long as both estimators are consistent under the null and their asymptotic bias are different under the alternative, it is still possible to construct a Hausman-type test even though both estimators are inconsistent under the alternative.

To demonstrate the general idea, consider the EwPO estimator with $w_{ij} = |\Delta x_{ij}|$, $\hat{\beta}_{1,(i,j)}$, Equation ((ref)) implies that

equation*[equation* omitted — 246 chars of source]

In the case that $u_i = \rho (x_i - \mu_x) + v_i$ with $x_i \perp v_i$, it is straightforward to show that $\epsilon_n - \delta_n(|\Delta x_{ij}|) = o_p(1)$ since both bias terms are $\rho + o_p(1)$. Thus, the test will not work with $\hat{\beta}_{1,(i,j)}$ in this case, which rule out the case when $x_i$ and $u_i$ follows a joint normal distribution.

However, when $u_i$ cannot be expressed as a linear function of $x_i$ then Hausman type test may be possible. One example is to consider $u_i = \rho x_i^3 + v_i$ where $x_i$ follows a symmetrical distribution with $\E(x_i) = 0$, $\E(x_i^4) = \kappa < \infty$ and $x_i \perp v_i$, then it is possible to show that $\epsilon_n = \kappa\rho + o_p(1)$ and $\delta_n (|\Delta x_{ij}|) = o_p(1)$ \footnote{EwPO with $w_{ij} = |\Delta x_{ij}|$ is consistent in this special case, despite correlation between $u_i$ and $x_i$}.

One possible way to construct the test statistics is to consider

equation[equation omitted — 275 chars of source]

where $\Delta \bfSigma \left [ \hat{\bfbeta}_1 - \hat{\bfbeta}_{1,(i,j)} \right ] = \bfSigma(\hat{\bfbeta}_1) - \bfSigma\left(\hat{\bfbeta}_{1,(i,j)}\right )$ with $\bfSigma(\hat{\bfbeta}_1)$ and $\bfSigma\left(\hat{\bfbeta}_{1,(i,j)} \right )$ denoting the variance-covariance matrices of $\hat{\bfbeta}_1$ and $\hat{\bfbeta}_{1, (i,j)}$, respectively.

The asymptotic distribution of the test statistics in the case $K=1$ with $w_{ij} = |\Delta x_{ij}|$ can be readily obtained using Proposition (ref) as demonstrated in Proposition (ref).

propositionUnder Assumptions (ref), (ref) and (ref) along with the assumption that $|\epsilon_n - \delta_n(|\Delta x_{ij}|)| = C + o_p(1)$ where $C>0$, then \begin{equation} H(|\Delta x_{ij}|) \overset{d}{\sim} \frac{\sigma_u^2}{d\mu^2_{\Delta x_s}} \left ( \frac{\mu_{\Delta x_s} - 2\sigma_x}{\sigma_x} W(1) - 4 \int^1_0 W(\lambda) d\lambda \right )^2 \end{equation} where $d = Var(\hat{\beta}_1) - Var(\hat{\beta}_{1,(i,j)})$.
proofSee Appendix (ref).

The result from Proposition (ref) is that the distribution of the test is non-standard and therefore its critical values must be simulated. Specifically, there are two components that would require simulation. First, the expression within the bracket of Equation ((ref)) and second, $Var \left ( \hat{\beta}_{1,(i,j)} \right )$, following the result from Proposition (ref). In both cases, Equations ((ref)) to ((ref)) can be used to carry out the simulation. Other quantities, such as $\sigma_u$, $\sigma_x$ and $\mu_{\Delta x_S}$, can be replaced by their sample estimates. Table (ref) provides simulated critical values for the case $\sigma_u = 1$, $\sigma_x=1$ and $\mu_{\Delta x_S} = \sqrt{2/\pi}$ with $n=5000$ and 10,000 replications.

table[table omitted — 399 chars of source]

The power of the test depends on the exact relation between $x_i$ and $u_i$ as well as the choice of $w_{ij}$ and is an interesting issue for further research.

Conclusion

Endogeneity has been an enduring problem in econometrics. The practice to compensate for the information loss due to the correlation between the explanatory variable(s) and the disturbance terms has been to use additional information in the form of instrumental variables or moment conditions for testing and consistent estimation. The problem with this approach is that the results hinge on the “quality” of this information. The common wisdom to the day is that testing for endogeneity is not possible based exclusively on the information provided by the observations of the variables a model. This paper through the introduction of a new estimation method, the so-called estimation with pairwise observations (EwPO), demonstrates that this belief is misconceived and introduces three different direct endogeneity tests, ready to be used in applications. The EwPO new approach may also open the door for further new research paths in econometric and statistical estimation theory.

\vskip 3cm \Appendix{Monte Carlo Simulations of EwPO Estimation with Selected Weights} {

This Appendix presents some Monte Carlo (MC) simulation results to examine the properties of some EwPO estimators with selected weights. Additional results can be found in the Online Supplement (see EwPO-Sup2023).

The data generating process for the MC simulations is based on the model

equation*[equation* omitted — 53 chars of source]

and the MC experiments consider two possible distributions for $u_i$, namely

enumerate\(u_{i} \sim N(0,1)\), • \(u_{i} \sim \) skewed normal distribution,

where the skewed normal distribution is generated as

equation*[equation* omitted — 53 chars of source]

with $\xi=-\lambda\sqrt{\frac{2}{\pi}}$, $v_i\sim N(0,1)$ and $z_i\sim N(0, \sigma^2)$ such that $v_i$ and $w_i$ are independently distributed.

The MC experiments consider uniform distribution U(-10,10) for the regressor $x_i$.

Two parameter vectors have been considered $(\beta_0, \beta_1) = (1,0.5)$ and $(\beta_0, \beta_1) = (1, 1.5)$ for purposes of robustness checking. The number of MC replications was 1000. For the sake of brevity, only the results for $(\beta_0, \beta_1) = (1,0.5)$ are reported here. The other results can be found in the Online Supplement (EwPO-Sup2023).

table[table omitted — 1,240 chars of source]
table[table omitted — 1,240 chars of source]

Note: No mistake, the sorted and non-sorted full-pairwise estimation results given in Tables (ref) and (ref) are exactly the same.

table[table omitted — 1,257 chars of source]
table[table omitted — 1,261 chars of source]

\Appendix{Monte Carlo Simulations for the Covariance Test}

The MC setup considers again sample sizes $n=50$, $500$ and $5000$ with $1000$ replications.

enumerate[{Step} 1.] • Generated the same model as above with one explanatory variable: $$ y_i = \beta_0 + \beta_1x_i + u_i $$ with $\beta_0=1$, $\beta_1=0.5$ to star with, and $u_i$ was generated as $N(0, 1)$. Finally, $x$ was generated as $N(0, 5)$ and also as $U(-5, 5)$. The simulations of $x_i$ and $u_i$ were carried out under four different correlations, namely $\rho=0$ (benchmark exogeneity), $\rho=0.2$ (small), $\rho=0.5$ (medium), and $\rho=0.8$ (large). • Estimate the model with EwPO with $w_{ij} = \Delta x_{ij}$ and $w_{ij} =|\Delta x_{ij}|$. In each case, calculate the test statistics as defined in Equation ((ref)).

For further simulation results please refer to the Online Supplement (EwPO-Sup2023).

table[table omitted — 1,322 chars of source]
figure[figure omitted — 265 chars of source]
figure[figure omitted — 269 chars of source]
figure[figure omitted — 270 chars of source]
table[table omitted — 1,335 chars of source]

\vskip 1cm

table[table omitted — 1,346 chars of source]

\Appendix{Monte Carlo Results of the Residual Test and the Consistent Estimation} The DGP in this case is

equation*[equation* omitted — 39 chars of source]

where $x_i \sim N(5,2)$ and $u_i \sim N(0,1)$. The MC experiments considered no correlation and $\rho = 0.2, 0.5$ and 0.8 with sample sizes ranging from 50 to 5000 for each case.

table[table omitted — 1,407 chars of source]
table[table omitted — 1,214 chars of source]
table[table omitted — 1,148 chars of source]

\eject \Appendix{Mathematical Appendix} Proof of Proposition (ref). Define

align[align omitted — 188 chars of source]

where $1/n \leq \lambda \leq 1$ and $\floor{x}$ denotes the largest integer contains in $x$. Under Assumptions ((ref)) to ((ref)) $B_n(\lambda) \overset{d}{\rightarrow} B(\lambda)$ and $W_n(\lambda) \overset{d}{\rightarrow} W(\lambda)$ following the arguments in mcleish_invariance_1975. The proof first identifies the asymptotic distribution of each term in Equation ((ref)) and the result follow by an application of the Continuous Mapping Theorem as shown in billingsley:1999. From Equation ((ref)), it is straightforward to show that

equation[equation omitted — 175 chars of source]

and follow the same argument, it can also be shown that

equation[equation omitted — 175 chars of source]

Now, note that $\Delta x_{pq}^2 = \left [ (x_p - \mu_x) - (x_q-\mu_x) \right ]^2$ and

align*[align* omitted — 538 chars of source]

The last line follows from a similar argument as Equations ((ref)) and ((ref)). Therefore,

equation[equation omitted — 187 chars of source]

For the last term, repeat the argument above and note that $u_p = W_n(\frac{p}{n})-W_n(\frac{p-1}{n})$

align*[align* omitted — 322 chars of source]

Therefore

equation[equation omitted — 311 chars of source]

Combining Equations ((ref)), ((ref)), ((ref)) and ((ref)) and apply the Continuus Mapping Theorem gives the result. This completes the proof. $\blacksquare$

Proof of Proposition (ref).

The key to the proof is the fact that $\sum^n_{i=2} \sum^{i-1}_{j=1} |\Delta x_{ij}| = \sum^n_{i=2} \sum^{i-1}_{j=1} \Delta x_{s,ij}$. Under Assumption (ref), the result as stated in Equation ((ref)) still holds since reordering $x_i$ does not affect the distribution of $u_i$. Thus

equation[equation omitted — 236 chars of source]

by the Law of Large Number for triangular arrays. Recall

equation*[equation* omitted — 169 chars of source]

Hence, combining Equations ((ref)) and ((ref)) and apply the Continuous Mapping Theorem gives the result. This completes the proof. $\blacksquare$

Proof of Proposition (ref) The proof follows the same argument as Proposition (ref). Proposition (ref) shows that the rate of convergence for EwPO with $w_{ij} = |\Delta x_{ij}|$ is $\sqrt{n}$ and therefore,

equation*[equation* omitted — 96 chars of source]

since EwPO is a consistent estimator under Assumption (ref) i.e., $\Delta \hat{u}_{pq} = \Delta u_{pq} + o_p(1)$. Therefore the result follows directly from Equation ((ref)). This completes the proof. $\blacksquare$

Proof of Proposition (ref) Using the results from Proposition (ref), the proof follows the same argument as Theorem 1 in phillips_understanding_1986. This completes the proof. $\blacksquare$

\vskip 3cm