# Estimation with Pairwise Observations

Felix Chan  Curtin University, Perth<EMAIL_ADDRESS>László Mátyás  Central European University, Budapest, Hungary and Vienna,
<EMAIL_ADDRESS>
###### Abstract

The 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, and potentially several others.

Acknowledgement: The authors would like to thank Tom Wansbeek for insightful
discussion. Contribution by Gyorgy Ruzicska in the early stages of this
project and Kristof Reizinger’s invaluable coding work for the simulations are
also kindly acknowledged.

##  1 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:

$y_{i}=\beta_{0}+\beta_{1}x_{i}+u_{i},\qquad u_{i}\sim
D(0,\sigma^{2}_{u}),\quad i=1,\ldots,n.$ (1)

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:

  * •

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 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 non-sorted adjacent.

  * •

All feasible combinations

$[(x_{i},y_{i});(x_{j},y_{j})],\,\,(i\neq j)$

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) full-
pairwise approach. Here, again we may consider sorting the observations based
on the values of $x_{i}$, which gives the 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 1: Sorted adjacent Figure 2: Non-sorted full-pairwise

Our main question of interest is how these “cluster-wise” or let us call them
pairwise parameters relate to the “true” parameters of the data generation
process formulated in (1). 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:

  * •

The weighted average of $\beta_{i,1}$

$\hat{\beta}_{1}=\frac{\sum_{i=1}^{n}w_{i}\beta_{i,1}}{\sum_{i=1}^{n}w_{i}}$
(2)

with some $w_{i}$ weights.

  * •

Minimise a loss function such as

$\hat{\beta}_{1}=\underset{\beta_{1}}{\arg\min}\sum_{i=1}^{n}(w_{i}(\beta_{i,1}-\hat{\beta}_{1}))^{2}.$
(3)

Estimators for $\beta_{0}$ can also be defined similarly and they will be
discussed with more details in Section 3.5.

The question now is how $\hat{\beta}_{1}$ relates to the true parameter,
$\beta_{1}$ formulated in Equation (1). Let us call this estimation procedure
Estimation with Pairwise Observations (EwPO).

##  2 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 A.333Additional simulation results can be found in the Online
Supplement Chan et al. (2023). 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:

  1. 1.

EwPO as defined in Equation (2) with

    1. (a)

Sorted full-pairwise and $w_{ij}=x_{i}-x_{j}$

    2. (b)

Non-sorted adjacent and $w_{i}=|x_{i}-x_{i-1}|$

    3. (c)

Sorted full-pairwise and $w_{ij}=|x_{i}-x_{j}|$

    4. (d)

Non-sorted full-pairwise and $w_{ij}=|x_{i}-x_{j}|$

    5. (e)

Sorted full-pairwise and $w_{ij}=\sqrt{(x_{i}-x_{j})^{2}+(y_{i}-y_{j})^{2}}$

    6. (f)

Non-sorted full-pairwise and
$w_{ij}=\sqrt{(x_{i}-x_{j})^{2}+(y_{i}-y_{j})^{2}}$

  2. 2.

Loss function defined in Equation (3)

    1. (a)

Non-sorted adjacent and $w_{i}=x_{i}-x_{i-1}$

    2. (b)

Sorted full-pairwise and $w_{ij}=x_{i}-x_{j}$

    3. (c)

Non-sorted full-pairwise and $w_{ij}=x_{i}-x_{j}$

Note: In case of the loss function taking the absolute value of weights does
not matter as weights are squared in the functional form!

##  3 Some Analytical Results

This section presents some analytical explanations on the simulations results.
It will also present the notations and definitions used throughout this paper.

###  3.1 Notations and Definitions

In order to develop the theoretical foundation, we introduce a few notations,
definitions and concepts. Let $\Delta x_{ij}$ denotes the difference between
$x_{i}$ and $x_{j}$ that is,

$\Delta x_{ij}=x_{i}-x_{j},\qquad i,j=1,\ldots,n.$ (4)

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

$\beta_{1,(i,j)}=\frac{\Delta y_{ij}}{\Delta x_{ij}}.$ (5)

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

$\Delta y_{ij}=\beta_{1}\Delta x_{ij}+\Delta u_{ij}.$ (6)

###  3.2 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 $\lfloor x\rfloor$ as the largest integer that
is less than or equal to $x$ and $\lceil x\rceil$ 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

$W_{n}(r)=\frac{1}{\sqrt{n}}\sum^{\lfloor
rn\rfloor}_{i=1}\frac{u_{i}}{\sigma_{u}},$ (7)

then under the conditional stated in (Billingsley, 1999, Theorem 8.2, P.90),
the Functional Central Limit Theorem ensures
$W_{n}(r)\overset{d}{\rightarrow}W(r)$, where $W(r)$ denotes a Wiener Process.

###  3.3 Adjacent Case

The estimator considered in this section is

$\hat{\beta}_{1}=\left(\sum^{n}_{i=2}w_{i}\right)^{-1}\left(\sum^{n}_{i=2}w_{i}\beta_{1,i}\right)$
(8)

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

####  3.3.1 Case 1: $w_{i}=\Delta x_{i}$

In this case, Equation (8) reduces to

$\displaystyle\hat{\beta}_{1}=$
$\displaystyle\left(\sum^{n}_{i=2}w_{i}\right)^{-1}\left(\sum^{n}_{i=2}w_{i}\beta_{1,i}\right)$
$\displaystyle=$ $\displaystyle\left(\sum^{n}_{i=2}\Delta
x_{i}\right)^{-1}\left(\sum^{n}_{i=2}\Delta y_{i}\right)$ $\displaystyle=$
$\displaystyle\left(\sum^{n}_{i=2}\Delta
x_{i}\right)^{-1}\left(\sum^{n}_{i=2}\beta_{1}\Delta x_{i}+\Delta
u_{i}\right)$ $\displaystyle=$
$\displaystyle\beta_{1}+\frac{\sum^{n}_{i=2}\Delta u_{i}}{\sum^{n}_{i=2}\Delta
x_{i}}.$

Now, in the case of adjacent pairwise parameter

$\displaystyle\sum^{n}_{i=2}\Delta x_{i}=$
$\displaystyle(x_{2}-x_{1})+(x_{3}-x_{2})+(x_{4}-x_{3})+\ldots+(x_{n}-x_{n-1})$
$\displaystyle=$ $\displaystyle x_{n}-x_{1}.$

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

$\hat{\beta}_{1}=\beta_{1}+\frac{u_{n}-u_{1}}{x_{n}-x_{1}}.$ (9)

Equation (9) 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 (9) 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}}$ will remain even when
$n\rightarrow\infty$.

####  3.3.2 Case 2: $w_{i}=|\Delta x_{i}|$

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

$\text{sgn}\left(x\right)=\begin{cases}1&\text{ if }x\geq 0\\\ -1&\text{ if
}x<0.\end{cases}$ (10)

Note that $\text{sgn}\left(x\right)x=|x|$, and that unlike the case when
$w_{i}=\Delta x_{i}$,

$\sum^{n}_{i=1}|\Delta
x_{i}|=|x_{2}-x_{1}|+|x_{3}-x_{2}|+\ldots+|x_{n}-x_{n-1}|$

and no simple cancellation occurs here when the observations are not sorted.
However, when the observations are 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,

$\displaystyle\hat{\beta}_{1}=$
$\displaystyle\left(\sum^{n}_{i=2}w_{i}\right)^{-1}\left(\sum^{n}_{i=2}w_{i}\beta_{1,i}\right)$
$\displaystyle=$ $\displaystyle\left(\sum^{n}_{i=2}|\Delta
x_{i}|\right)^{-1}\left(\sum^{n}_{i=2}\text{sgn}\left(\Delta
x_{i}\right)\Delta y_{i}\right)$ $\displaystyle=$
$\displaystyle\left(\sum^{n}_{i=2}|\Delta
x_{i}|\right)^{-1}\left(\sum^{n}_{i=2}\beta_{1}\text{sgn}\left(\Delta
x_{i}\right)\Delta x_{i}+\text{sgn}\left(\Delta x_{i}\right)\Delta
u_{i}\right)$ $\displaystyle=$ $\displaystyle\left(\sum^{n}_{i=2}|\Delta
x_{i}|\right)^{-1}\left(\sum^{n}_{i=2}\beta_{1}|\Delta
x_{i}|+\text{sgn}\left(\Delta x_{i}\right)\Delta
u_{i}\right)\quad\because\text{sgn}\left(\Delta x_{i}\right)\Delta
x_{i}=|\Delta x_{i}|$ $\displaystyle=$
$\displaystyle\beta_{1}+\frac{\sum^{n}_{i=2}\text{sgn}\left(\Delta
x_{i}\right)u_{i}}{\sum^{n}_{i=1}|\Delta x_{i}|}.$

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

$(n-1)^{-1}\sum^{n}_{i=2}\text{sgn}\left(\Delta x_{i}\right)u_{i}=o_{p}(1)$

and hence

$\hat{\beta}_{1}-\beta_{1}=o_{p}(1).$

###  3.4 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

$\hat{\beta}_{1}=\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}\beta_{1,(i,j)}\right)$
(11)

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

####  3.4.1 Case 1: $w_{ij}=\Delta x_{ij}$

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

$\displaystyle\hat{\beta}_{1}=$
$\displaystyle\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}\beta_{1,(i,j)}\right)$
$\displaystyle=$ $\displaystyle\left(\sum^{n}_{i=2}\sum^{i-1}_{j=1}\Delta
x_{ij}\right)^{-1}\left(\sum^{n}_{i=2}\sum^{i-1}_{j=1}\Delta y_{ij}\right)$
$\displaystyle=$
$\displaystyle\beta_{1}+\frac{\sum^{n}_{i=2}\sum^{i-1}_{j=1}\Delta
u_{ij}}{\sum^{n}_{i=2}\sum^{i-1}_{j=1}\Delta x_{ij}}.$ (12)

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

$\displaystyle\sum^{n}_{i=2}\sum^{i-1}_{j=1}\Delta x_{ij}=$
$\displaystyle\sum^{n}_{i=2}\sum^{i-1}_{j=1}(x_{i}-x_{j})$ $\displaystyle=$
$\displaystyle\sum^{n}_{i=2}\sum^{i-1}_{j=1}x_{i}-\sum^{n}_{i=2}\sum^{i-1}_{j=1}x_{j}$
$\displaystyle=$
$\displaystyle\sum^{n}_{i=2}(i-1)x_{i}-\sum^{n}_{i=2}\sum^{i-1}_{j=1}x_{j}.$
(13)

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
$\mathbb{E}|x_{i}|<\infty$. We will analyse expression (13) term-by-term. The
first term in expression (13) can be rewritten as

$\displaystyle\sum^{n}_{i=2}(i-1)x_{i}=$ $\displaystyle
N\mu_{x}+\sum^{n}_{i=2}(i-1)(x_{i}-\mu_{x})$ $\displaystyle=$ $\displaystyle
N\mu_{x}+\sum^{n}_{i=2}iz_{i}-\sum^{n}_{i=2}z_{i}$ $\displaystyle=$
$\displaystyle N\mu_{x}-\sum^{n}_{i=2}z_{i}+\sum^{n}_{i=2}iz_{i}\,,$

where $z_{i}=x_{i}-\mu_{x}$. Note that $z_{i}$ is the de-meaned version of
$x_{i}$ and therefore $\mathbb{E}(z_{i})=0$. To expand further the expression
above,

$\displaystyle N\mu_{x}-\sum^{n}_{i=2}z_{i}+\sum^{n}_{i=2}iz_{i}=$
$\displaystyle N\mu_{x}-\sum^{n}_{i=2}z_{i}+\sum^{n}_{i=1}\sum^{n}_{j=i}z_{j}$
$\displaystyle=$ $\displaystyle
N\mu_{x}-\sum^{n}_{i=2}z_{i}+\sum^{n}_{i=1}\sum^{n}_{j=\lceil
rN\rceil}z_{j}\text{ where }r=i/N$ $\displaystyle=$ $\displaystyle
N\mu_{x}-\sum^{n}_{i=2}z_{i}+\sqrt{n}\sigma_{x}\sum^{n}_{i=1}\left(W_{n}(1)-W_{n}(r)\right)$
$\displaystyle=$ $\displaystyle
N\mu_{x}-\sum^{n}_{i=2}z_{i}+n\sqrt{n}\sigma_{x}\sum^{n}_{i=1}\left(\int^{i/n}_{(i-1)/n}W_{n}(1)-W_{n}(s)ds\right)$
$\displaystyle=$ $\displaystyle
N\mu_{x}-\sum^{n}_{i=2}z_{i}+n\sqrt{n}\sigma_{x}\left(W_{n}(1)-\int^{1}_{0}W_{n}(s)ds\right).$
(14)

Now for the second term in expression (13),

$\displaystyle\sum^{n}_{i=2}\sum^{i-1}_{j=1}x_{j}=$ $\displaystyle
N\mu_{x}+\sum^{n}_{i=2}\sum^{i-1}_{j=1}z_{j}$ $\displaystyle=$ $\displaystyle
N\mu_{x}+\sqrt{n}\sigma_{x}\sum^{n}_{i=1}W_{n}(r)\quad r=\frac{i-1}{n}$
$\displaystyle=$ $\displaystyle
N\mu_{x}+n\sqrt{n}\sigma_{x}\sum^{n}_{i=2}\int^{i/n}_{(i-1)/n}W_{n}(s)ds$
$\displaystyle=$ $\displaystyle
N\mu_{x}+n\sqrt{n}\sigma_{x}\int^{1}_{0}W_{n}(s)ds..$ (15)

Substitute expressions (14) and (15) to expression (13) yields

$\sum^{n}_{i=2}\sum^{i-1}_{j=1}\Delta
x_{ij}=n\sqrt{n}\sigma_{x}\left(W(1)-2\int^{1}_{0}W(s)ds\right)-\sum^{n}_{i=1}z_{i}.$
(16)

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

$N^{-1}\sum^{n}_{i=2}\sum^{i-1}_{j=1}\Delta x_{ij}=o_{p}(1).$

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 (12) 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 (12)
approaches 0 as $n\rightarrow\infty$ but the denominator will either approach
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 1,
which provides the asymptotic distribution for this case.

####  3.4.2 Case 2: $w_{ij}=|\Delta x_{ij}|$

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

$\displaystyle\hat{\beta}_{1,(i,j)}=$
$\displaystyle{\left(\sum^{n}_{i=2}\sum^{i-1}_{j=1}|\Delta
x_{ij}|\right)}^{-1}\left(\sum^{n}_{i=2}\sum^{i-1}_{j=1}\frac{\Delta
y_{ij}}{\Delta x_{ij}}|\Delta x_{ij}|\right)$ $\displaystyle=$
$\displaystyle{\left(\sum^{n}_{i=2}\sum^{i-1}_{j=1}|\Delta
x_{ij}|\right)}^{-1}\left(\sum^{n}_{i=2}\sum^{i-1}_{j=1}\Delta
y_{ij}\text{sgn}\left(\Delta x_{ij}\right)\right)$ $\displaystyle=$
$\displaystyle{\left(\sum^{n}_{i=2}\sum^{i-1}_{j=1}|\Delta
x_{ij}|\right)}^{-1}\left(\sum^{n}_{i=2}\sum^{i-1}_{j=1}\beta_{1}|\Delta
x_{ij}|+\Delta u_{ij}\text{sgn}\left(\Delta x_{ij}\right)\right)$
$\displaystyle=$
$\displaystyle{\beta_{1}+\left(\sum^{n}_{i=2}\sum^{i-1}_{j=1}|\Delta
x_{ij}|\right)}^{-1}\left(\sum^{n}_{i=2}\sum^{i-1}_{j=1}\Delta
u_{ij}\text{sgn}\left(\Delta x_{ij}\right)\right)$ $\displaystyle=$
$\displaystyle\beta_{1}+\left(\sum^{n}_{i=2}\sum^{i-1}_{j=1}|\Delta
x_{ij}|\right)^{-1}\left(\sum^{n}_{i=2}\sum^{i-1}_{j=1}u_{i}\text{sgn}\left(\Delta
x_{ij}\right)+\sum^{n}_{i=2}\sum^{i-1}_{j=1}u_{j}\text{sgn}\left(\Delta
x_{ij}\right)\right)\,.$ (17)

Given that $x_{i}\perp u_{j}$ for all $i,j=1,\ldots,n$ and
$n^{-1}\sum^{n}_{i=1}|\text{sgn}\left(\Delta x_{ij}\right)|\leq 1$, the last
two terms in Equation (17) will therefore converge 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)$.

####  3.4.3 Confidence Intervals

The asymptotic distribution of $\hat{\beta}$ in the full-pairwise case with
$w_{ij}=\Delta x_{ij}$ can be found in Proposition 2. 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 will 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 will be the same across different choices of $w_{ij}$.

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

  1. Step 1.

Set the level of significant for the test $\alpha$.

  2. Step 2.

Set the number of observations to be removed from the main sample, call it $d$
where $\sqrt{n}<d<n$.

  3. Step 3.

Set the number of replication $R<\begin{pmatrix}n\\\ n-d\end{pmatrix}$

  4. Step 4.

set $i=1$

  5. Step 5.

Create a jackknife sample, $J_{i}$, by removing $d$ randomly selected
observations from the main sample.

  6. Step 6.

Compute $\hat{\beta}$ using the jackknife sample $J_{i}$, call it
$\hat{\beta}_{i}$

  7. Step 7.

Set $i\leftarrow i+1$.

  8. Step 8.

Repeat Steps 5 to 7 until $i=R$.

  9. Step 9.

Sort the set $\\{\hat{\beta}_{i}\\}_{i=1}^{R}$ from lowest to highest.

  10. Step 10.

The $\displaystyle\lfloor\frac{\alpha}{2}R\rfloor^{th}$ and
$\lfloor\frac{(1-\alpha)}{2}R\rfloor^{th}$ elements in the set gives the lower
and upper bounds of the confidence interval at the $\alpha$ significant level.

Table 1 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.

Full-pairwise MC | Jackknife confidence intervals lower and upper bounds
---|---
|  Parameter |  Full-pairwise one single sample MC estimates |  Lower bound ($\beta_{0.025}$) |  Upper bound ($\beta_{0.975}$)
n=50 |  Exogen |  0.4479 |  0.3142 |  0.5618
$\rho=0.2$ |  0.6085 |  0.4763 |  0.8134
$\rho=0.5$ |  0.5925 |  0.4515 |  0.7357
$\rho=0.8$ |  0.5971 |  0.4415 |  0.7929
n=500 |  Exogen |  0.4725 |  0.4541 |  0.5352
$\rho=0.2$ |  0.5347 |  0.4994 |  0.5809
$\rho=0.5$ |  0.6051 |  0.5881 |  0.6683
$\rho=0.8$ |  0.6516 |  0.6188 |  0.6971
n=1000 |  Exogen |  0.5027 |  0.4388 |  0.5003
$\rho=0.2$ |  0.5558 |  0.5106 |  0.5660
$\rho=0.5$ |  0.6101 |  0.5626 |  0.6195
$\rho=0.8$ |  0.6527 |  0.6404 |  0.6912
n=5000 |  Exogen |  0.5117 |  0.4899 |  0.5155
$\rho=0.2$ |  0.5415 |  0.5291 |  0.5554
$\rho=0.5$ |  0.6007 |  0.5904 |  0.6155
$\rho=0.8$ |  0.6644 |  0.6561 |  0.6801
Table 1: Jakckknife simulation of $\beta$ quantiles for one given sample of
the Monte Carlo (MC) simulation results reported in the $2^{\text{nd}}$
column, DGP $\sim$ N(0,5), $|\Delta x|$ weighted full-pairwise estimator,
number of repetitions: 10000 and jackknife sample size ($d$): 2535

###  3.5 What about the Intercept?

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

$\hat{\beta}_{0}=\bar{y}_{n}-\hat{\beta}_{1}\bar{x}_{n}\,,$ (18)

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

$\displaystyle\hat{\beta}_{0}=$
$\displaystyle\bar{y}_{n}-\hat{\beta}_{1}\bar{x}_{n}$ $\displaystyle=$
$\displaystyle\left(\beta_{1}-\hat{\beta}_{1}\right)\bar{x}_{n}+\bar{u}_{n}$

and as $n\rightarrow\infty$, $\beta_{1}-\hat{\beta}_{1}=o_{p}(1)$,
$\bar{x}_{n}-\mathbb{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

$\beta_{0,(i,j)}=y_{i}-\beta_{1,(i,j)}x_{i}.$ (19)

####  3.5.1 The adjacent case

In the adjacent case we consider the following estimator of $\beta_{0}$

$\hat{\beta}_{0}={\left(\sum^{n}_{i=1}w_{i}\right)}^{-1}\left(\sum^{n}_{i=1}w_{i}\beta_{0,i}\right).$
(20)

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}$.

$\displaystyle\hat{\beta}_{0}=$ $\displaystyle{\left(\sum^{n}_{i=2}|\Delta
x_{i}|\right)}^{-1}\left(\sum^{n}_{i=2}|\Delta x_{i}|\beta_{0,i}\right)$
$\displaystyle=$ $\displaystyle\beta_{0}+{\left(\sum^{n}_{i=2}|\Delta
x_{i}|\right)}^{-1}\left[\beta_{1}\sum^{n}_{i=2}|\Delta
x_{i}|x_{i}-\sum^{n}_{i=2}|\Delta x_{i}|x_{i}\beta_{1,i}+\sum^{n}_{i=2}|\Delta
x_{i}|u_{i}\right].$ (21)

Expanding the following term further gives

$\displaystyle\sum^{n}_{i=2}|\Delta x_{i}|x_{i}\beta_{1,i}=$
$\displaystyle\sum^{n}_{i=2}\text{sgn}\left(\Delta x_{i}\right)x_{i}\Delta
y_{i}$ $\displaystyle=$ $\displaystyle\beta_{1}\sum^{n}_{i=2}|\Delta
x_{i}|x_{i}+\sum^{n}_{i=2}\text{sgn}\left(\Delta x_{i}\right)x_{i}\Delta
u_{i}.$

Substitute the last expression above into Equation (21) yields

$\hat{\beta}_{0}=\beta_{0}+{\left(\sum^{n}_{i=2}|\Delta
x_{i}|\right)}^{-1}\left(\sum^{n}_{i=2}|\Delta
x_{i}|u_{i}-\sum^{n}_{i=2}\text{sgn}\left(\Delta x_{i}\right)x_{i}\Delta
u_{i}\right).$ (22)

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}\text{sgn}\left(\Delta
x_{i}\right)x_{i}\Delta u_{i}$. The former holds, because the sum is
converging to $\mathbb{E}\left(|\Delta
x_{i}|u_{i}\right)=\mathbb{E}\left(|\Delta x_{i}|\right)\mathbb{E}(u_{i})=0$
and the second expression holds because the sum is converging to
$\mathbb{E}\left[\text{sgn}\left(\Delta x_{i}\right)x_{i}\Delta
u_{i}\right]=\mathbb{E}\left[\text{sgn}\left(\Delta
x_{i}\right)x_{i}\right]\mathbb{E}\left(\Delta u_{i}\right)=0$.

###  3.6 Loss Function

Next, we analyse the properties of the estimator (3):

$\hat{\beta}=\underset{\beta}{\arg\min}\sum^{n}_{i=1}\left[w_{i}\left(\beta_{i}-\beta\right)\right]^{2}\,,$
(23)

where

$\beta_{i}=\frac{\Delta y_{i}}{\Delta x_{i}}.$

####  3.6.1 $w_{i}=|\Delta x_{i}|$

Given the quadratic loss, this section will 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 (23) implies

$\sum^{n}_{i=1}w_{i}^{2}\left(\beta_{i}-\hat{\beta}\right)=0,$

and therefore

$\hat{\beta}=\left(\sum^{n}_{i=1}w_{i}^{2}\right)^{-1}\left(\sum^{n}_{i=1}w_{i}^{2}\beta_{i}\right).$

Thus when $w_{i}=|\Delta x_{i}|$ or when $w_{i}=\Delta x_{i}$ we get

$\hat{\beta}=\left(\sum^{n}_{i=1}|\Delta
x_{i}|^{2}\right)^{-1}\left(\sum^{n}_{i=1}\Delta x_{i}\Delta y_{i}\right)\,,$

which is in fact the OLS estimator for the ‘first difference’ equation

$\Delta y_{i}=\Delta x_{i}\beta+\Delta u_{i}\,.$ (24)

####  3.6.2 Full-pairwise

Reconsider the estimator in the full-pairwise case

$\hat{\beta}=\underset{\beta}{\arg\min}\sum^{n}_{i=1}\sum^{i-1}_{j=1}\left[w_{ij}\left(\beta_{ij}-\beta\right)\right]^{2}\,,$
(25)

which gives

$\hat{\beta}=\left(\sum^{n}_{i=1}\sum^{i-1}_{j=1}w_{ij}^{2}\right)^{-1}\left(\sum^{n}_{i=1}\sum^{i-1}_{j=1}w_{ij}^{2}\Delta\beta_{ij}\right).$

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

$\hat{\beta}=\beta+\left(\sum^{n}_{i=1}\sum^{i-1}_{j=1}|\Delta
x_{ij}|^{2}\right)^{-1}\left(\sum^{n}_{i=1}\sum^{i-1}_{j=1}\Delta x_{ij}\Delta
u_{ij}\right)\,,$ (26)

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

$\hat{\beta}=\left(\sum^{n}_{i=1}\sum^{i-1}_{j=1}|\Delta
x_{ij}|\right)^{-1}\left(\sum^{n}_{i=1}\sum^{i-1}_{j=1}|\Delta
x_{ij}|\hat{\beta}_{ij}\right)\,,$

which is our full-pairwise estimator.

##  4 Multivariate Extension

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

$y_{i}=\beta_{0}+\sum^{K}_{k=1}\beta_{k}x_{ki}+u_{i}\qquad u_{i}\sim
D(0,\sigma^{2}_{u})\quad i=1,\ldots,n$ (27)

with the matrix representation

$\mathbf{y}=\mathbf{X}\beta+u\,,$ (28)

where $\mathbf{y}=\left(y_{1},\ldots,y_{n}\right)^{\prime}$ is a $n\times 1$
vector containing the dependent variable,
$\mathbf{X}=\left(\mathbf{x}_{1},\ldots,\mathbf{x}_{K}\right)$ is a $n\times
K$ matrix containing the $K$ explanatory variables such that $\mathbf{X}$ is
$\mathbf{x}_{k}=\left(x_{k1},\ldots,x_{kn}\right)^{\prime}$ for $k=1,\ldots,K$
and $\mathbf{u}=\left(u_{1},\ldots,u_{K}\right)^{\prime}$ is a $n\times 1$
vector containing the residuals.

It will be convenient to introduce some special matrices here. These matrices
allow us to represent the pairwise parameter as a sequence of matrix
operations. Define $\mathbf{d}_{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
$\mathbf{D}_{F}=\left(\mathbf{d}_{n,n-1}^{\prime},\mathbf{d}_{n,n-2}^{\prime},\ldots,\mathbf{d}_{n,1}^{\prime},\mathbf{d}_{n-1,n-2}^{\prime},\ldots,\mathbf{d}_{n-1,1}^{\prime},\ldots,\mathbf{d}_{2,1}^{\prime}\right)^{\prime}$
and
$\mathbf{D}_{a}=\left(\mathbf{d}_{2,1}^{\prime},\mathbf{d}_{3,2}^{\prime},\mathbf{d}_{4,3}^{\prime}\ldots,\mathbf{d}_{n,n-1}^{\prime}\right)^{\prime}$.
Given these matrices, it is straightforward to show that
$\mathbf{D}_{F}\mathbf{x}=\\{\Delta x_{ij}\\}_{i=2,j=1}^{n,i-1}$ is a
$n(n-1)/2\times 1$ vector and $\mathbf{D}_{a}\mathbf{x}=\\{\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 $\mathbf{S}:\mathbb{R}^{n}\rightarrow\mathbb{R}^{n}$ be a function that
returns a $n\times 1$ zero-one selection matrix such that if
$\mathbf{z}=\mathbf{S}(\mathbf{x})\mathbf{x}$ then $\mathbf{z}$ will be a
$n\times 1$ vector containing the same elements as $\mathbf{x}$ such that
$z_{i}\geq z_{j}$ for $j>i$. Thus, $\mathbf{S}(\mathbf{x})\mathbf{x}$ produces
the sorted version of $\mathbf{x}$ in descending order. Finally, define the
residual maker

$\mathbf{M}_{k}=\mathbf{I}_{n}-\mathbf{X}_{-k}{\left(\mathbf{X}_{-k}^{\prime}\mathbf{X}_{-k}\right)}^{-1}\mathbf{X}_{-k}\,,$
(29)

where $\mathbf{X}_{-k}$ denote the $n\times K-1$ matrix which contains all the
columns of $\mathbf{X}$ 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 will construct the multivariate version by
combining the univariate estimator for each of the regressors.

$\mathbf{y}=\mathbf{1}\beta_{0}+\mathbf{x}\boldsymbol{\beta}_{1}+\mathbf{u}$

and note that
$\mathbf{D}\mathbf{y}=\mathbf{D}\mathbf{x}\beta_{1}+\mathbf{D}\mathbf{u}$ for
$\mathbf{D}=\\{\mathbf{D}_{a},\mathbf{D}_{F}\\}$. The pairwise parameters for
the sample can be written as

$\boldsymbol{\beta}_{1}=\text{{diag}}\left(\mathbf{x}\right)^{-1}\mathbf{D}\mathbf{y}\,,$
(30)

where
$\text{{diag}}\left(\mathbf{x}\right)={\left[\mathbf{I}\otimes\left(\mathbf{D}\mathbf{x}\right)^{\prime}\right]\mathbf{S}_{n}}$
is a diagonal matrix with the elements of $\mathbf{x}$ in the diagonal.
$\mathbf{S}_{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 (30) is reasonably general and covers all the
previously discussed case by defining the transformation matrix $\mathbf{D}$
differently. A summary can be found in Table 2.

Cases | $\mathbf{D}$
---|---
Adjacent | $\mathbf{D}_{a}$
Adjacent Sorted | $\mathbf{D}_{a}\mathbf{S}(\mathbf{x})$
Full-pairwise | $\mathbf{D}_{F}$
Full-pairwise Sorted | $\mathbf{D}_{F}\mathbf{S}(\mathbf{x})$
Table 2: Definitions of $\mathbf{D}$ based on different construction of
pairwise parameters

Given the vector of pairwise parameters $\boldsymbol{\beta}_{1}$, the
estimator of $\beta_{1}$ can be written as

$\hat{\beta}_{1}={\left(\mathbf{w}_{1}\mathbf{i}^{\prime}\right)}^{-1}\mathbf{w}_{1}^{\prime}\boldsymbol{\beta}_{1}\,,$
(31)

where $\mathbf{w}_{1}$ is the vector of weights associated with each
$\beta_{1,(i,j)}$ in $\boldsymbol{\beta}_{1}$. Following from the analysis
above, some of the choices include $\mathbf{w}_{1}=\left(|\Delta
x_{1}|,\ldots,|\Delta x_{n}|\right)^{\prime}$ and $\mathbf{w}_{1}=\left(\Delta
x_{1},\ldots,\Delta x_{n}\right)^{\prime}$.

It is straightforward to generalise Equation (30), and subsequently Equation
(31), to the case when $\mathbf{X}$ is a $n\times K$ matrix with $K$
regressors. Consider

$\displaystyle\mathbf{y}=$
$\displaystyle\mathbf{X}\boldsymbol{\beta}+\mathbf{u}$ (32) $\displaystyle=$
$\displaystyle\mathbf{x}_{k}\beta_{k}+\mathbf{X}_{-k}\boldsymbol{\beta}_{-k}+\mathbf{u}$
(33) $\displaystyle\mathbf{M}_{k}\mathbf{y}=$
$\displaystyle\mathbf{M}_{k}\mathbf{x}_{k}\beta_{k}+\mathbf{u}\,,$ (34)

where $\boldsymbol{\beta}_{-k}$ contains the same elements as
$\boldsymbol{\beta}$ with the $k^{th}$ element removed. Following Equation
(30) the pairwise parameter of $\beta_{k}$ is

$\boldsymbol{\beta}_{k}=\text{{diag}}\left(\mathbf{M}_{k}\mathbf{x}_{k}\right)^{-1}\mathbf{D}\mathbf{M}_{k}\mathbf{y}.$
(35)

Similarly, Equation (31) can be generalised as

$\hat{\beta}_{k}={\left(\mathbf{w}_{k}\mathbf{i}^{\prime}\right)}^{-1}\mathbf{w}_{k}^{\prime}\boldsymbol{\beta}_{k}$
(36)

and

$\hat{\boldsymbol{\beta}}=\begin{bmatrix}{\left(\mathbf{w}_{1}\mathbf{i}^{\prime}\right)}^{-1}\mathbf{w}_{1}^{\prime}&\mathbf{0}&\ldots&\mathbf{0}&\mathbf{0}\\\
\mathbf{0}&\ddots&\vdots&\ldots&\mathbf{0}\\\
\mathbf{0}&\ldots&{\left(\mathbf{w}_{k}\mathbf{i}^{\prime}\right)}^{-1}\mathbf{w}_{k}^{\prime}&\ldots&\mathbf{0}\\\
\vdots&\ldots&\vdots&\ddots&\ldots\\\
\mathbf{0}&\ldots&\ldots&\ldots&{\left(\mathbf{w}_{K}\mathbf{i}^{\prime}\right)}^{-1}\mathbf{w}_{K}^{\prime}\end{bmatrix}\begin{bmatrix}\boldsymbol{\beta}_{1}\\\
\vdots\\\ \boldsymbol{\beta}_{k}\\\ \vdots\\\
\boldsymbol{\beta}_{K}\end{bmatrix}.$ (37)

As shown in previous results, suitable choice of $\mathbf{w}_{k}$ leads to
consistency of $\hat{\beta}_{k}$ and therefore
$\hat{\boldsymbol{\beta}}-\boldsymbol{\beta}=o_{p}(1)$ under the same
conditions as those discussed above for all $k=1,\ldots,K$.

##  5 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., Hausman (1978), or Wooldridge (2002), pp.
118-122).

The aim here is to develop two tests for endogeneity based on the estimator of
the form

$\hat{\beta}=\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}\beta_{ij}\right)$

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

###  5.1 Residuals Test

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

$\hat{\beta}=\beta+\delta_{n}\,,$ (38)

where

$\delta_{n}=\left(\sum^{n}_{i=2}\sum^{i-1}_{j=1}|\Delta
x_{ij}|\right)^{-1}\left(\sum^{n}_{i=2}\sum^{i-1}_{j=1}\text{sgn}\left(\Delta
x_{ij}\right)\Delta u_{ij}\right)$ (39)

and under $x_{i}\perp u_{i}$, $\delta_{n}=o_{p}(1)$. Now given the following
specification

$y_{i}=x_{i}\beta+u_{i}$ (40)

the estimated residual from the EwPO estimator is

$\displaystyle\hat{u}_{i}=$ $\displaystyle y_{i}-x_{i}\hat{\beta}$
$\displaystyle=x_{i}(\beta-\hat{\beta})+u_{i}$
$\displaystyle=-x_{i}\delta_{n}+u_{i}.$

It is straightforward to show that

$\displaystyle n^{-1}\sum^{n}_{i=1}\hat{u}_{i}=$
$\displaystyle-n^{-1}\sum^{n}_{i=1}x_{i}\delta_{n}+u_{i}$ $\displaystyle=$
$\displaystyle-\delta_{n}n^{-1}\sum^{n}_{i=1}x_{i}+n^{-1}\sum^{n}_{i=1}u_{i}\,.$

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

$n^{-1}\sum^{n}_{i=1}\hat{u}_{i}=-\delta_{n}\mu_{x}+o_{p}(1)$ (41)

under the assumption that $\mathbb{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 is always possible by appropriate
transformations on the regressors, it is therefore possible to directly test
for endogeneity by testing if the mean of the estimated residuals is
statistically different from 0 using standard testing procedure, such as the
t-test.444The 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.

###  5.2 Covariance Test

Note that $\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

$\begin{split}\Delta\hat{u}_{pq}=&\Delta y_{pq}-\Delta x_{pq}\hat{\beta}\\\
=&\Delta x_{pq}\left(\beta-\hat{\beta}\right)+\Delta u_{pq}.\end{split}$ (42)

In the case of a consistent estimator $\hat{\beta}$ can be expressed in the
form of $\hat{\beta}=\beta+\delta(\mathbf{w},\mathbf{u})$ where $\mathbf{w}$
and $\mathbf{u}$ 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

$\Delta\hat{u}_{pq}=\Delta x_{pq}\delta(\mathbf{w},\mathbf{u})+\Delta u_{pq}.$

So consider the test statistics

$\begin{split}S(\mathbf{w})=n^{-2}\sum^{n}_{p=2}\sum^{p-1}_{q=1}\Delta
x_{pq}\Delta\hat{u}_{pq}.\end{split}$ (43)

####  5.2.1 Case 1: $\mathbf{w}=\\{\Delta x_{ij}\\}_{i,j=1}^{n}$

In this case,

$\Delta\hat{u}_{pq}=\Delta x_{pq}\left(\sum^{n}_{i=2}\sum^{i-1}_{j=1}\Delta
x_{ij}\right)^{-1}\left(\sum^{n}_{i=2}\sum^{i-1}_{j=1}\Delta
u_{ij}\right)+\Delta u_{pq}$

and the test statistics

$\begin{split}S(\mathbf{w})=&n^{-2}\left(\sum^{n}_{p=1}\sum^{p-1}_{q=1}\Delta
x_{pq}^{2}\right)\left(\sum^{n}_{i=2}\sum^{i-1}_{j=1}\Delta
x_{ij}\right)^{-1}\left(\sum^{n}_{i=2}\sum^{i-1}_{j=1}\Delta u_{ij}\right)\\\
&\qquad+n^{-2}\left(\sum^{n}_{p=2}\sum^{p-1}_{q=1}\Delta x_{pq}\Delta
u_{pq}\right)\,.\end{split}$ (44)

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 (11) with
$w_{ij}=\Delta x_{ij}$ and that of the test statistics as defined above can be
found in Propositions 1 and 2, respectively.

Before presenting the theoretical results, consider the following assumptions

  1. Assumption 1.

$\mathbb{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)$.

  2. Assumption 2.

$\mathbb{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)$.

  3. Assumption 3.

$\mathbb{E}(u_{i}|x_{j})=0$ for $i,j=1,\ldots,n$.

  4. Assumption 4.

$\Delta x_{ij}$ and $\Delta u_{ij}$ are mixingales as defined in McLeish
(1975) and satisfy all the conditions for Theorem 3.8 in McLeish (1975).

Assumptions 1 and 2 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 3 is required for the consistency of the
estimator. It also stipulates the null underlying the distribution of the test
statistics as asserted in Proposition 2. Assumption 4 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 and
White (1988) for potentially more general assumptions on $x_{i}$ and $u_{i}$.

Asymptotic distribution of the full-pairwise EwPO with $w_{ij}=\Delta x_{ij}$
can be found in Proposition 1 below

######  Proposition 1.

Under Assumptions 1 to 4 with $\hat{\beta}$ as defined in Equation (11) and
$\mathbf{w}=\\{\Delta x_{ij}\\}^{n}_{i,j=1}$:

$\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}.$
(45)

###### Proof.

See Appendix. ∎

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

######  Proposition 2.

Under the assumptions of Proposition 1

$\begin{split}(\sigma_{x}\sigma_{u})^{-1}S(\mathbf{w})\overset{d}{\rightarrow}&B(1)W(1)+\int^{1}_{0}B(\lambda)W(\lambda)d\lambda-\int^{1}_{0}B(\lambda)dW(\lambda)\\\
&\quad-\int^{1}_{0}W(\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}$
(46)

where $B(\lambda)$ and $W(\lambda)$ are independent standard Brownian motion
processes.

###### Proof.

See Appendix C ∎

Remark 1. It can be shown that both terms on the right hand side of Equation
(46) are negative on expectation. Under the alternative that $u_{i}$ and
$x_{i}$ are correlated, then the values of the last two stochastic integrals
in Equation (46) increase as the correlation increases.

Remark 2. The asymptotic distribution as stated in Equation (46) 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

$\displaystyle\int^{1}_{0}B(\lambda)d\lambda\approx$
$\displaystyle\sum^{n}_{p=1}z_{1p}$ (47)
$\displaystyle\int^{1}_{0}W(\lambda)d\lambda\approx$
$\displaystyle\sum^{n}_{p=1}z_{3p}$ (48)

where $z_{ip}=z_{ip-1}+\varepsilon_{ip}$,
$\boldsymbol{\varepsilon}_{p}=(\varepsilon_{1p},\varepsilon_{2p})^{\prime}\sim
NID(0,\mathbf{I})$ and

$\displaystyle\int^{1}_{0}W(\lambda)dB(\lambda)\approx$
$\displaystyle\sum^{n}_{p=1}z_{2p}\varepsilon_{1p}$ (50)
$\displaystyle\int^{1}_{0}B(\lambda)dW(\lambda)\approx$
$\displaystyle\sum^{n}_{p=1}z_{1p}\varepsilon_{2p}.$ (51)

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 (47) - (51) and the substitution their
values into Equation (46) 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 (46) using Equations (47) - (51) with $\sigma_{x}=\sigma_{u}=1$ and
$n=10000$ can be found in Table 3 and the respective histogram in Figure 3.
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.

Figure 3: Simulated Sample of the Distribution $\alpha$ | Lower CV | Upper CV
---|---|---
0.01 % | -4.129 | 3.913
0.05 % | -3.021 | 2.912
0.10 % | -2.530 | 2.498
Table 3: Simulated Critical Values: $w_{ij}=\Delta x_{ij}$

####  5.2.2 Case 2: $w_{ij}=|\Delta x_{ij}|$

Similar derivation as above yields

$\begin{split}S(\mathbf{w})=&n^{-2}\left(\sum^{n}_{p=1}\sum^{p-1}_{q=1}\Delta
x_{pq}^{2}\right)\left(\sum^{n}_{i=2}\sum^{i-1}_{j=1}|\Delta
x_{ij}|\right)^{-1}\left(\sum^{n}_{i=2}\sum^{i-1}_{j=1}\text{sgn}\left(\Delta
x_{ij}\right)\Delta u_{ij}\right)\\\
&\qquad+n^{-2}\left(\sum^{n}_{p=2}\sum^{p-1}_{q=1}\Delta x_{pq}\Delta
u_{pq}\right)\,.\end{split}$ (52)

Given the absolute value of the weight, the asymptotic distribution of the
test statistics is more challenging to obtain. However, the jackknife
procedure similar to those proposed in Section 3.4.3 can easily be used to
obtain the critical value of the test statistics in this case as well.

Appendix B 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 B, the mean test statistics deviates from zeros as the
correlation between the regressor and the error term grows. Further
investigations would suggest that the covariance test has the right size with
relatively good power, even in small sample.

###  5.3 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:

  1. a)

Estimate the model, e.g., $y_{i}=\beta x_{i}+u_{i}$ with EwPO;

  2. b)

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_{i}y_{i}$, $w_{i}=g_{i}x_{i}$ and $v_{i}=g_{i}u_{i}$ with $g_{i}$ an
IV candidate that satisfies the same assumptions as $x_{i}$;

  3. c)

Select the model for inference with the lowest statistics value. This will
minimise 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.

##  6 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 and consistent estimation of a regression
model with 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 shows how to test directly for endogeneity. It also derives a simple new
consistent estimation procedure for the linear regression model with
endogeneity, which can easily be implemented in practice.

##  Appendix A Monte Carlo Simulations Setups and Some Simulation Results for
the EwPO Estimation

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 Chan et al. (2023)).

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

$y_{i}=\beta_{0}+\beta_{1}x_{i}+u_{i}$

and the MC experiments consider two possible distributions for $u_{i}$, namely

  1. 1.

$u_{i}\sim N(0,1)$,

  2. 2.

$u_{i}\sim$ skewed normal distribution,

where the skewed normal distribution is generated as

$u_{i}=\xi+\lambda|v_{i}|+z_{i},$

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 (Chan et al. (2023)).

Sorted MC | Estimates and MC standard errors
---|---
n=50 | Parameter | Estimate/S.e. | OLS | pairwise
$\hat{\beta_{0}}$ | Estimate | 0.9866 | 0.9866
S.e. | 0.1609 | 0.1609
$\hat{\beta_{1}}$ | Estimate | 0.4999 | 0.4999
| S.e. | 0.2928 | 0.2928
n = 500 | $\hat{\beta_{0}}$ | Estimate | 1.0001 | 1.0001
S.e. | 0.0504 | 0.0504
$\hat{\beta_{1}}$ | Estimate | 0.4971 | 0.4971
S.e. | 0.0948 | 0.0948
n = 5000 | $\hat{\beta_{0}}$ | Estimate | 0.9992 | 0.9992
S.e. | 0.0169 | 0.0169
$\hat{\beta_{1}}$ | Estimate | 0.4997 | 0.4997
S.e. | 0.0281 | 0.0281
Table 4: Sorted – full-pairwise MC, $x_{i}$ $\sim U(-10,10)$, $u_{i}\sim$ skewed normal distribution, $\Delta x$ weighted estimator Sorted MC | Estimates and MC standard errors
---|---
n=50 | Parameter | Estimate/S.e. | OLS | pairwise
$\hat{\beta_{0}}$ | Estimate | 0.9866 | 0.9866
S.e. | 0.1609 | 0.1609
$\hat{\beta_{1}}$ | Estimate | 0.4999 | 0.4999
| S.e. | 0.2928 | 0.2928
n = 500 | $\hat{\beta_{0}}$ | Estimate | 1.0001 | 1.0001
S.e. | 0.0504 | 0.0504
$\hat{\beta_{1}}$ | Estimate | 0.4971 | 0.4971
S.e. | 0.0948 | 0.0948
n = 5000 | $\hat{\beta_{0}}$ | Estimate | 0.9992 | 0.9992
S.e. | 0.0169 | 0.0169
$\hat{\beta_{1}}$ | Estimate | 0.4997 | 0.4997
S.e. | 0.0281 | 0.0281
Table 5: Non-sorted full-pairwise MC, $x_{i}$ $\sim U(-10,10)$, $u_{i}\sim$
skewed normal distribution, $\Delta x$ weighted estimator
Note: No mistake, the sorted and non-sorted full-pairwise estimation results given in Tables 4 and 5 are exactly the same.  Non-sorted MC | Estimates and MC standard errors
---|---
n=50 | Parameter | Estimate/S.e. | OLS | pairwise
$\hat{\beta_{0}}$ | Estimate | 1.0009 | 0.9967
S.e. | 0.1408 | 0.1997
$\hat{\beta_{1}}$ | Estimate | 0.5007 | 0.5018
| S.e. | 0.0251 | 0.0291
n = 500 | $\hat{\beta_{0}}$ | Estimate | 0.9967 | 0.9966
S.e. | 0.0446 | 0.0640
$\hat{\beta_{1}}$ | Estimate | 0.5001 | 0.4998
S.e. | 0.0081 | 0.0095
n = 5000 | $\hat{\beta_{0}}$ | Estimate | 1.0001 | 1.0004
S.e. | 0.0144 | 0.0192
$\hat{\beta_{1}}$ | Estimate | 0.4999 | 0.4999
S.e. | 0.0025 | 0.0031
Table 6: Non-sorted adjacent MC, $x_{i}\sim U(-10,10)$, $u_{i}\sim N(0,1)$, $|\Delta x|$ weighted estimator Full-pairwise MC | Estimates and MC standard errors
---|---
n=50 | Parameter | Estimate/S.e. | OLS | pairwise
$\hat{\beta_{0}}$ | Estimate | 1.0009 | 1.0009
S.e. | 0.1408 | 0.1408
$\hat{\beta_{1}}$ | Estimate | 0.5007 | 0.5007
| S.e. | 0.0251 | 0.0251
n = 500 | $\hat{\beta_{0}}$ | Estimate | 0.9967 | 0.9967
S.e. | 0.0446 | 0.0446
$\hat{\beta_{1}}$ | Estimate | 0.5001 | 0.5001
S.e. | 0.0081 | 0.0081
n = 5000 | $\hat{\beta_{0}}$ | Estimate | 1.0001 | 1.0001
S.e. | 0.0144 | 0.0144
$\hat{\beta_{1}}$ | Estimate | 0.4999 | 0.4999
S.e. | 0.0025 | 0.0025
Table 7: Sorted full-pairwise MC, $x_{i}\sim U(-10,10)$, $u_{i}\sim N(0,1)$,
$\Delta x$ weighted estimator

##  Appendix B Monte Carlo Simulations for the Covariance Test

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

  1. Step 1.

Generated the same model as above with one explanatory variable:

$y_{i}=\beta_{0}+\beta_{1}x_{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).

  2. Step 2.

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
(43).

For further simulation results please refer to the Online Supplement (Chan et
al. (2023)).

Full-pairwise MC | Average test statistics
---|---
n=50 | Parameter | Pairwise |  Standard deviation | Skewness | Kurtosis
Exogen | -0.0003 |  0.7842 | -0.1804 | 3.1856
$\rho=0.2$ | -0.5478 |  0.8647 | 0.0273 | 3.1118
$\rho=0.5$ | -0.6175 |  0.5483 | 0.0593 | 3.0018
| $\rho=0.8$ | -1.2625 |  0.4008 | 0.0122 | 3.1264
n = 500 | Exogen | -0.0058 |  0.2175 | -0.0497 | 3.0011
$\rho=0.2$ | -0.4241 |  0.2075 | -0.1293 | 3.0698
$\rho=0.5$ | -1.1112 |  0.2017 | 0.0756 | 3.0814
$\rho=0.8$ | -1.5581 |  0.1311 | -0.0891 | 2.8930
n = 5000 | Exogen | -0.0031 |  0.0666 | 0.0101 | 2.8834
$\rho=0.2$ | -0.3995 |  0.0684 | 0.0123 | 3.1566
$\rho=0.5$ | -0.9737 |  0.0604 | -0.0610 | 2.9180
$\rho=0.8$ | -1.5681 |  0.0421 | 0.0637 | 3.0883
Table 8: Average test statistics, full-pairwise MC, $x_{i}\sim N(0,5)$,
$\Delta x$ weighted estimator Figure 4: Distribution of the test-statistics,
$x_{i}\sim$ N(0,5), $\Delta x$ weighted full-pairwise estimator, $n=50$ Figure
5: Distribution of the test-statistics, $x_{i}\sim$ N(0,5), $\Delta x$
weighted full-pairwise estimator, $n=500$ Figure 6: Distribution of the test-
statistics, $x_{i}\sim$ N(0,5), $\Delta x$ weighted full-pairwise estimator,
$n=5000$

Full-pairwise MC | Average test statistics
---|---
n=50 | Parameter | Pairwise |  Standard deviation | Skewness | Kurtosis
Exogen | 0.0060 |  0.8890 | 0.0382 | 3.2387
$\rho=0.2$ | -0.9247 |  0.7920 | -0.0244 | 2.9620
$\rho=0.5$ | -2.6616 |  0.7437 | 0.0121 | 3.2542
| $\rho=0.8$ | -5.9118 |  0.5729 | -0.0305 | 3.0115
n = 500 | Exogen | -0.0234 |  0.2750 | -0.1429 | 3.1852
$\rho=0.2$ | -1.0846 |  0.2471 | 0.0120 | 2.7804
$\rho=0.5$ | -2.9665 |  0.2309 | 0.1625 | 3.2859
$\rho=0.8$ | -4.2625 |  0.1494 | 0.0345 | 3.1253
n = 5000 | Exogen | -0.0036 |  0.0785 | 0.0088 | 2.8812
$\rho=0.2$ | -1.1328 |  0.0786 | -0.0625 | 3.1951
$\rho=0.5$ | -2.8321 |  0.0710 | -0.0394 | 2.9287
$\rho=0.8$ | -4.5415 |  0.0494 | 0.1022 | 3.0687
Table 9: Average test statistics, full-pairwise MC, $x_{i}\sim U(-5,5)$,
$\Delta x$ weighted estimator

Full-pairwise MC | Average test statistics
---|---
n=50 | Parameter | Pairwise |  Standard deviation | Skewness | Kurtosis
Exogen | 0.0191 |  0.8476 | -0.1264 | 3.1170
$\rho=0.2$ | -1.2884 |  0.8520 | -0.1054 | 2.8218
$\rho=0.5$ | -2.8584 |  0.7488 | -0.0452 | 2.9483
| $\rho=0.8$ | -6.9515 |  0.6224 | -0.0069 | 3.2675
n = 500 | Exogen | 0.0010 |  0.2608 | -0.0092 | 2.9270
$\rho=0.2$ | -1.2375 |  0.2588 | 0.0632 | 2.7668
$\rho=0.5$ | -2.7598 |  0.2218 | -0.0266 | 2.9173
$\rho=0.8$ | -4.3437 |  0.1479 | 0.0420 | 2.9936
n = 5000 | Exogen | -0.0003 |  0.0784 | 0.1196 | 2.8769
$\rho=0.2$ | -1.1162 |  0.0790 | -0.1253 | 2.8993
$\rho=0.5$ | -2.8175 |  0.0703 | -0.1483 | 3.0588
$\rho=0.8$ | -4.5073 |  0.0490 | -0.0011 | 2.9966
Table 10: Average test statistics, full-pairwise MC, $x_{i}\sim U(-5,5)$,
$|\Delta x|$ weighted estimator

##  Appendix C Monte Carlo Results of the Residual Test and the Consistent
Estimation

The DGP in this case is

$y_{i}=0.5x_{i}+u_{i}$

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.

Sample size | Correlation | Mean | Variance | Skewness | Kurtosis
---|---|---|---|---|---
n=50 | Exogen | 0.0044 | 0.1604 | -0.0176 | 0.1476
$\rho=0.2$ | -0.4930 | 0.1546 | -0.0546 | 0.0601
$\rho=0.5$ | -1.2506 | 0.1180 | -0.0113 | 0.1834
$\rho=0.8$ | -2.0033 | 0.0606 | 0.0223 | 0.1587
n=500 | Exogen | 0.0032 | 0.0149 | 0.0839 | 0.0023
$\rho=0.2$ | -0.5005 | 0.0142 | 0.0264 | -0.0169
$\rho=0.5$ | -1.2517 | 0.0118 | 0.023 | 0.1509
$\rho=0.8$ | -2.0015 | 0.0051 | -0.0321 | -0.2207
n=1000 | Exogen | 0.0009 | 0.0075 | -0.0792 | 0.0818
$\rho=0.2$ | -0.5027 | 0.0073 | -0.0534 | 0.1680
$\rho=0.5$ | -1.2474 | 0.0056 | 0.0687 | 0.0894
$\rho=0.8$ | -2.0003 | 0.0028 | 0.0432 | -0.0320
n=5000 | Exogen | 0.0003 | 0.0015 | -0.0341 | 0.1900
$\rho=0.2$ | -0.4993 | 0.0014 | 0.0541 | 0.1994
$\rho=0.5$ | -1.2507 | 0.0012 | 0.0573 | 0.1016
$\rho=0.8$ | -2.0003 | 0.0005 | 0.0722 | -0.0429
Table 11: MC Results on the Distributions of $\hat{u}_{i}$. Sample size | Correlation | Mean | Variance | Skewness | Kurtosis
---|---|---|---|---|---
$n=50$ | Exogen | 0.4993 | 0.0056 | 0.0584 | 0.2268
$\rho=0.2$ | 0.5984 | 0.0054 | 0.0641 | 0.067
$\rho=0.5$ | 0.7505 | 0.004 | 0.013 | 0.0595
$\rho=0.8$ | 0.9001 | 0.0021 | -0.0344 | 0.3254
$n=500$ | Exogen | 0.4994 | 0.0005 | -0.0126 | 0.0119
$\rho=0.2$ | 0.5998 | 0.0005 | 0.0287 | -0.1062
$\rho=0.5$ | 0.7505 | 0.0004 | -0.0401 | 0.0729
$\rho=0.8$ | 0.9003 | 0.0002 | 0.0101 | -0.1041
$n=1000$ | Exogen | 0.4998 | 0.0003 | 0.0331 | 0.0474
$\rho=0.2$ | 0.6006 | 0.0003 | 0.0736 | 0.2013
$\rho=0.5$ | 0.7498 | 0.0002 | -0.0369 | -0.0212
$\rho=0.8$ | 0.9001 | 0.0001 | -0.0594 | -0.0057
$n=5000$ | Exogen | 0.5000 | 0.0001 | 0.0240 | 0.0152
$\rho=0.2$ | 0.5999 | 4.877e-5 | -0.0424 | 0.1421
$\rho=0.5$ | 0.7502 | 4.055e-5 | -0.0478 | 0.2321
$\rho=0.8$ | 0.9000 | 1.850e-5 | -0.0246 | -0.0227
Table 12: Monte Carlo Results on $\hat{\beta}$ Sample size | Correlation | Mean | Variance | Skewness | Kurtosis
---|---|---|---|---|---
$n=50$ | Exogen | 0.5002 | 0.0008 | 0.0503 | -0.0639
$\rho=0.2$ | 0.4998 | 0.0008 | 0.0861 | 0.2677
$\rho=0.5$ | 0.5004 | 0.0006 | 0.0015 | -0.0670
$\rho=0.8$ | 0.4995 | 0.0003 | 0.0377 | 0.0951
$n=500$ | Exogen | 0.5000 | 8.031e-5 | -0.0401 | 0.0146
$\rho=0.2$ | 0.4997 | 7.297e-5 | -0.1208 | 0.0535
$\rho=0.5$ | 0.5002 | 6.249e-5 | 0.0534 | 0.0321
$\rho=0.8$ | 0.5000 | 2.946e-5 | 0.0825 | 0.1877
$n=1000$ | Exogen | 0.5000 | 3.975e-5 | -0.0374 | -0.1315
$\rho=0.2$ | 0.5000 | 3.871e-5 | 0.0496 | -0.0301
$\rho=0.5$ | 0.5003 | 2.94e-5 | -0.1190 | 0.0958
$\rho=0.8$ | 0.5001 | 1.435e-5 | 0.0450 | 0.1941
$n=5000$ | Exogen | 0.5000 | 8.225e-6 | 0.0703 | -0.0280
$\rho=0.2$ | 0.5000 | 7.724e-6 | 0.0096 | 0.0661
$\rho=0.5$ | 0.5000 | 5.731e-6 | -0.0170 | -0.0944
$\rho=0.8$ | 0.5000 | 2.933e-6 | 0.0195 | 0.0369
Table 13: Bias-corrected $\hat{\beta}$

##  Appendix D Mathematical Appendix

Proof of Proposition 2.

Define

$\displaystyle B_{n}(\lambda)=$ $\displaystyle
n^{-1}\sum_{i=1}^{\lfloor\lambda n\rfloor}\frac{x_{i}-\mu_{x}}{\sigma_{x}}$
(53) $\displaystyle W_{n}(\lambda)=$ $\displaystyle
n^{-1}\sum_{i=1}^{\lfloor\lambda n\rfloor}\frac{u_{i}}{\sigma_{u}}$ (54)

where $1/n\leq\lambda\leq 1$ and $\lfloor x\rfloor$ denotes the largest
integer contains in $x$. Under Assumptions (1) to (4)
$B_{n}(\lambda)\overset{d}{\rightarrow}B(\lambda)$ and
$W_{n}(\lambda)\overset{d}{\rightarrow}W(\lambda)$ following the arguments in
McLeish (1975).

The proof first identifies the asymptotic distribution of each term in
Equation (43) and the result follow by an application of the Continuous
Mapping Theorem as shown in Billingsley (1999). From Equation (16), it is
straightforward to show that

$n^{-\frac{3}{2}}\sigma_{x}^{-1}\sum^{n}_{i=2}\sum^{i-1}_{j=1}\Delta
x_{ij}\overset{d}{\rightarrow}B(1)-2\int^{1}_{0}B(\lambda)d\lambda$ (56)

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

$n^{-\frac{3}{2}}\sigma_{u}^{-1}\sum^{n}_{i=2}\sum^{i-1}_{j=1}\Delta
u_{ij}\overset{d}{\rightarrow}W(1)-2\int^{1}_{0}W(\lambda)d\lambda.$ (57)

Now, note that $\Delta
x_{pq}^{2}=\left[(x_{p}-\mu_{x})-(x_{q}-\mu_{x})\right]^{2}$ and

$\displaystyle\sum^{n}_{i=2}\sum^{i-1}_{j=1}\Delta
x^{2}_{pq}=\sigma^{2}_{x}\sum^{n}_{i=2}\sum^{i-1}_{j=1}\left(\frac{x_{i}-\mu_{x}}{\sigma_{x}}-\frac{x_{j}-\mu_{x}}{\sigma_{x}}\right)^{2}=$
$\displaystyle=$
$\displaystyle\sigma^{2}_{x}\left[\sum^{n}_{i=2}\sum^{i-1}_{j=1}\left(\frac{x_{i}-\mu_{x}}{\sigma_{x}}\right)^{2}+\sum^{n}_{i=2}\sum^{i-1}_{j=1}\left(\frac{x_{j}-\mu_{x}}{\sigma_{x}}\right)^{2}-\sum^{n}_{i=2}\sum^{i-1}_{j=1}\frac{(x_{i}-\mu_{x})(x_{j}-\mu_{x})}{\sigma^{2}_{x}}\right]$
$\displaystyle=$
$\displaystyle\sigma^{2}_{x}n^{2}\left(B^{2}_{n}(1)-2\int^{1}_{0}B^{2}_{n}(\lambda)d\lambda\right)$

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

$n^{-2}\sum^{n}_{p=2}\sum^{p-1}_{q=1}\Delta
x^{2}_{pq}\overset{d}{\rightarrow}\sigma^{2}_{x}\left(B^{2}(1)-2\int^{1}_{0}B^{2}(\lambda)d\lambda\right).$
(58)

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

$\displaystyle\sum^{n}_{p=2}\sum^{p-1}_{q=1}\Delta x_{pq}\Delta u_{pq}=$
$\displaystyle\sum^{n}_{p=2}\sum^{p-1}_{q=1}x_{p}u_{p}+x_{q}u_{q}-x_{p}u_{q}-x_{q}u_{p}$
$\displaystyle=$ $\displaystyle
n^{2}\sigma_{x}\sigma_{u}\left(B_{n}(1)W_{n}(1)-\int^{1}_{0}W_{n}(\lambda)dB_{n}(\lambda)-\int^{1}_{0}B_{n}(\lambda)dW_{n}(\lambda)\right).$

Therefore

$n^{-2}\sum^{n}_{p=2}\sum^{p-1}_{q=1}\Delta
x_{pq}u_{pq}\overset{d}{\rightarrow}\sigma_{x}\sigma_{u}\left(B(1)W(1)-\int^{1}_{0}W(\lambda)dB(\lambda)-\int^{1}_{0}B(\lambda)dW(\lambda)\right).$
(59)

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

## REFERENCES

  * Billingsley (1999) Billingsley, Patrick (1999): _Convergence of Probability Measures, Second Edition_ , Wiley Inter-Science.
  * Chan et al. (2023) Chan, Felix, László Mátyás, and Kristóf Reizinger (2023): “Estimation with Pairwise Observations – Online Supplement,” _Github Repository_ , Click here to navigate to the online supplement.
  * Hausman (1978) Hausman, Jerry A. (1978): “Specification Tests in Econometrics,” _Econometrica_ , 1251–1271.
  * Johansen (1995) Johansen, S. (1995): _”Likelihood-based Inference in Cointegrated Vector Auto-regressive Models_ , Oxford University Press.
  * McLeish (1975) McLeish, D L (1975): “Invariance Principles for Dependent Variables,” _Zeitschrift Fur Wahrscheinlichkeitstheorie Und Verwandte Gebiete_ , 32 (3), 165–178.
  * Wooldridge (2002) Wooldridge, Jeffrey M. (2002): _Econometric Analysis of Cross Section and Panel Data_ , MIT Press.
  * Wooldridge and White (1988) Wooldridge, Jeffrey M and Halbert White (1988): “Some invariance principles and central limit theorems for dependent heterogeneous processes,” _Econometric Theory_ , 210–230, publisher: JSTOR.