Test for differences between parameters of models estimated from partially overlaping samples Suppose I have $n$ observations $(\pmb{x}_i,y_i)$ and two subsets of indexes of $\{1:n\}$, $S_1$ and $S_2$, with $S_1\neq S_2$, $\#\{S_1\}\neq\#\{S_2\}$, and $\{S_1\cap S_2\}\neq\emptyset$.
Suppose I call $\hat{\beta}_1$ the vector of OLS-estimated coefficients of the 
regression $y_i\tilde{}\pmb{x}_i|i\in S_1$ and likewise for $\hat{\beta}_2$. 
My question is  what is the right/best test for 
$$H_0:\beta_1-\beta_2=0$$
Some thoughts:


*

*My first inclination was for a Wald test, but how do we combine the two covariance matrices from the two regressions? Perhaps by pooling them?

*Another idea I thought about was to use the Bhattacharyya distance with $\beta_1$ against $\beta_2$ (together with their respective covariance matrices).

*A third possibility would be to merge the two samples and add an interaction effect on every coefficient with a dummy taking value 1 for all observations that are, say, members of $S_2$ only, but because of the overlap between, $S_1$ and $S_2$ i'm not sure this is really measuring what I want. 


Any thoughts?
EDIT1:
Here is a simple R code to implement the solution proposed 
by Bill below:
library("MASS")
n<-200
p<-5

b<-rnorm(p)
x<-mvrnorm(n,rep(0,p-1),diag(p-1))
y<-cbind(1,x)%*%b+rnorm(n)
L<-sample(1:200,100,replace=TRUE)
L<-list(mod1=L[1:60],mod2=L[61:100])

M1<-lm(y[L[[1]]]~x[L[[1]],])
M2<-lm(y[L[[2]]]~x[L[[2]],])

fx01<-function(M1,M2){
    X1<-cbind(1,M1$model[,-1])
     X2<-cbind(1,M2$model[,-1])
    nt<-nrow(X1)+nrow(X2)
    D1<-diag(nt)[1:nrow(X1),]   
    D2<-diag(nt)[1:nrow(X2),]   
    st<-crossprod(c(M1$resid,M2$resid))/(nt-2*ncol(X1))
    P1<-solve(crossprod(X1))%*%crossprod(X1,D1)-solve(crossprod(X2))%*%crossprod(X2,D2)
    P1<-tcrossprod(P1)*drop(st)
    mahalanobis(M1$coef,M2$coef,P1)
}

 A: The asymptotic Wald test is what you want.  You want to calculate:
$W=\left(\hat{\beta}_1-\hat{\beta}_2 \right)'\left( V(\hat{\beta}_1-\hat{\beta}_2 )\right)^{-1} \left( \hat{\beta}_1-\hat{\beta}_2 \right)$
How to calculate the variance in the middle, though?  Write the difference in the beta hats like this:
$\begin{align}
\hat{\beta}_1-\hat{\beta}_2 &= \left( X_1'X_1\right)^{-1}X_1'Y_1 -\left( X_2'X_2\right)^{-1}X_2'Y_2\\
&= \beta_1-\beta_2 + \left( X_1'X_1\right)^{-1}X_1'\epsilon_1 -\left( X_2'X_2\right)^{-1}X_2'\epsilon_2
\end{align}$
In the above, $\epsilon_1$ are the error terms from model 1, and $\epsilon_2$ are the errors from model 2.  The right hand side variables from model 1 are in the $(\#(S_1) \times K_1)$ matrix $X_1$.  Below, I'm going to use $\epsilon$ to be all the unique errors (i.e. an n by 1 vector).
Next, define a matrix $\Delta_1$ which "chooses" the observations in $S_1$ from the full set of observations.  The $(\#(S_1) \times n)$ matrix $\Delta_1$ looks like an identity matrix with the rows for the observations which are not in $S_1$ removed.  Substituting in . . .
$\begin{align}
\hat{\beta}_1-\hat{\beta}_2 &= \beta_1-\beta_2 + \left( X_1'X_1\right)^{-1}X_1'\Delta_1\epsilon -\left( X_2'X_2\right)^{-1}X_2'\Delta_2\epsilon\\
&= \beta_1-\beta_2 + \left( \left( X_1'X_1\right)^{-1}X_1'\Delta_1 -\left( X_2'X_2\right)^{-1}X_2'\Delta_2\right)\epsilon
\end{align}$
The variance of that is simple to calculate:
$V(\hat{\beta}_1-\hat{\beta}_2 )=\sigma^2\left( \left( X_1'X_1\right)^{-1}X_1'\Delta_1 -\left( X_2'X_2\right)^{-1}X_2'\Delta_2\right)\left( \left( X_1'X_1\right)^{-1}X_1'\Delta_1 -\left( X_2'X_2\right)^{-1}X_2'\Delta_2\right)'$
Sigma squared is the variance of the error term.  You have to estimate it.  Estimate it as:
$\frac{1}{\#(S_1)+\#(S_2)-K_1-K_2}\sum e_i^2$
where $K_1$ is the number of columns in $X_1$  The sum runs over all the residuals from both models (so there will be "double counting").
