function [b,V] = conley(y,X0,loc,year,locCutoff)

% Function calculates non-parametric spatial corrlation structure
%   using a panel data set (this one estimates the variance-covarinace matrix for each year) 
% See: TG Conley "GMM Estimation with Cross Sectional Dependence" 
%      Journal of Econometrics, Vol. 92 Issue 1 (September 1999) 1-45
%      http://www.elsevier.com/homepage/sae/econworld/econbase/econom/frame.htm
%
% Arguments: y: dependent variable (Nx1) vector 
%            X: independnet variables (Nxk) vector (INCLUDE constant as column)
%            loc: (Nxr) vector of locations (r=2 means plane, r=3 means 3D) 
%            year: (Nx1) vector of year indicators
%            locCutoff: (1xr) vector of cutoffs for nonparametric windows (one for each spatial dimension)
%
% By: Ruben Lebowski and Wolfram Schlenker

% rescale X 
Xs = X0;
% kill zeros
Xs(Xs==0) = 1;
Xmin = log10(abs(min(Xs)));
Xmax = log10(abs(max(Xs)));
Xscale = 10.^round(-(Xmin + Xmax)/2);
clear Xs Xmin Xmax;
X = X0.*repmat(Xscale,size(X0,1),1);

% first-stage OLS
b = X\y;

% get information on number of variables
K = size(X,2);
spatialDim = size(loc,2);    
% set variance-covariance matrix equal to zeros
XeeX = zeros(K);

% ------------------------------------------------------------
% loop over different years
yearUnique = sort(unique(year));
Nyears = max(size(yearUnique));
for ti = 1:Nyears
    rows = (year == yearUnique(ti));
    % get subset of variables
    X1 = X(rows,:);
    loc1 = loc(rows,:);
    e1 = y(rows) - X1*b;
    % get information on number of variables
    N = size(X1,1);
   
    % loop over all observations
    for i = 1:N
        % ----------------------------------------------------------------
        % step a: get non-parametric weight
        % Bartlett window of Newey-West
        window = ones(N,1);
        % loop over all spatial dimensions
        for j = 1:spatialDim
            window = window.*max( (1 - (abs(loc1(:,j)-loc1(i,j)))/locCutoff(j)) , 0);
        end
        % ----------------------------------------------------------------
        % step b: construct X'e'eX for given observation
        XeeXh = repmat(X1(i,:)',1,N).*repmat(e1(i)*e1'.*window',K,1) * X1;
        XeeX = XeeX + (XeeXh + XeeXh')/2;
    end
end
XeeX = XeeX/size(X,1);

% now get corrected variance-covariance matrix
invXX = inv(X'*X) * size(X,1);
V = invXX * XeeX * invXX / size(X,1);
% rescale
V = V.*(Xscale'*Xscale);
b = b.*Xscale';