capture mata mata drop altrdvce()
version 15.0
mata
real matrix altrdvce(real matrix X, real matrix y, real matrix z, real scalar p, real scalar h, real scalar matches, string vce, string kernel)
{
n = length(X)
p1 = p+1
if (vce=="resid") {
mu0_phat_y = mu0_phat_z = J(n,1, .)
sigma    = J(n,1, .)
	for (k=1; k<=n; k++) {
		cutoff = X[k]
		W = kweight(X,cutoff,h,kernel)
		W_p = select(W, W:> 0)
		X_p = select(X, W:> 0)
		y_p = select(y, W:> 0)
		z_p = select(z, W:> 0)
		XX  = J(length(W_p),p1,.)
		for (j=1; j<=p1; j++) {
			XX[.,j] = (X_p:-cutoff):^(j-1)
		}
	m_p_y = invsym(cross(XX,W_p,XX))*cross(XX,W_p,y_p)
	m_p_z = invsym(cross(XX,W_p,XX))*cross(XX,W_p,z_p)
	mu0_phat_y[k] = m_p_y[1]
	mu0_phat_z[k] = m_p_z[1]
	sigma[k] = (y[k] - mu0_phat_y[k])*(z[k] - mu0_phat_z[k])
	}
}
else  {

X_range = abs(X[n]-X[1])

X_shift_mat = J(n,(2*matches),.)
y_shift_mat = J(n,(2*matches),.)
z_shift_mat = J(n,(2*matches),.)

X_shift_mat[,1] = .\X[1::(n-1)]
X_shift_mat[,(matches+1)] = X[2::n]\.

y_shift_mat[,1] = .\y[1::(n-1)]
y_shift_mat[,(matches+1)] = y[2::n]\.

z_shift_mat[,1] = .\z[1::(n-1)]
z_shift_mat[,(matches+1)] = z[2::n]\.

for (k=2; k<=matches; k++) {

	X_shift_mat[,k] = .\X_shift_mat[1::(n-1),(k-1)]
	y_shift_mat[,k] = .\y_shift_mat[1::(n-1),(k-1)]
	z_shift_mat[,k] = .\z_shift_mat[1::(n-1),(k-1)]

							}

for (k=(matches+2); k<=(2*matches); k++) {

	X_shift_mat[,k] = X_shift_mat[2::n,(k-1)]\.
	y_shift_mat[,k] = y_shift_mat[2::n,(k-1)]\.
	z_shift_mat[,k] = z_shift_mat[2::n,(k-1)]\.

							}

ones = J(1,(2*matches),1)

X_mat = X * ones

X_diff = abs(X_mat-X_shift_mat)

num_matches = J(n,1,0)

y_match = J(n,matches,.)
z_match = J(n,matches,.)

for (k=1; k<=matches; k++)	{

	min_ind = (X_diff:==(rowmin(X_diff)*ones))

	num_matches_lt3 = num_matches:<3

	min_ind = (num_matches_lt3*ones):*min_ind

	y_match[,k] = rowsum(min_ind:*y_shift_mat)
	z_match[,k] = rowsum(min_ind:*z_shift_mat)

	X_diff = min_ind*(X_range+1) + (1:-min_ind):*X_diff

	num_matches = num_matches + rowsum(min_ind)

							}

y_match_ave = rowsum(y_match):/num_matches
z_match_ave = rowsum(z_match):/num_matches


sigma = (num_matches:/(num_matches:+1)):*(y - y_match_ave):*(z - z_match_ave)

}
return(sigma)
}
mata mosave altrdvce(), replace
end