function [beta, margEff, density] = KleinSpady(X,D,opt,h)

% X: n-by-k matrix
% D: n-by-1 vector
% h: a positive scalar
% opt: 'kr' or [] if (local constant) kernel regression, 'krll' if local linear kernel regression

% [beta, margEff, density] = KleinSpady(X,D,opt,h) implements the Klein-Spady estimator for a semiparametric
% binary choice model D=1(X*beta + epsilon >=0). The distribution function of epsilon is unknown and
% is estimated nonparametrically, i.e., by kernel or local linear regression using the Quartic kernel.

% The Klein and Spady estimator can only identify beta up to location and scale, so the last coefficient is normalized to 1.
% The bandwidth h, when not provided, is computed by Silverman's rule.

% Outputs allowed: coefficients, marginal effects, and the density evaulated at the index mean (Xbar*beta). 

% beta = KleinSpady(X,D) implements kernel regression to estimate the distribution funcion of epsilon.
% beta = KleinSpady(X,D,opt) implements the specified (in opt) nonparametric regression to estimate the distribution funcion of epsilon. 
% beta = KleinSpady(X,D,opt,h)) implements the specified (in opt) nonparametric regression to estimate
% the distribution funcion of epsilon and uses the provided bandwidth, h.

% Copyright: Yingying Dong at Boston College, July, 2008.

%%
% check inputs and outputs
error(nargchk(2,4,nargin));
error(nargchk(1,3,nargout));

if nargin < 4
    h =  [];
    display('The optimal bandwidth is computed using Silverman''s rule.')
elseif ~isscalar(h)
    error('The provided bandwidth must be a scalar.')
end

if nargin < 3
    opt = [];
    display('The error distribution function is estimated by kernel regression.')
    display('The optimal bandwidth is computed using Silverman''s rule.')
elseif ~(strcmp(opt, 'kr')||strcmp(opt,'krll')||isempty(opt))
    error('opt must be ''kr'', ''krll'' or [].')
end

%check and clean invalid data
inv = (sum(X, 2)~=sum(X, 2))|(D~=D);
X(inv) = [];
D(inv) = [];

[n,k]= size(X);

%---------Obtain starting values of coeff.s from probit----------------
options=optimset('MaxFunEvals', 1000, 'MaxIter', 1000,'TolX',1e-12,'TolFun',1e-12, 'Display', 'off' );
[opt_beta0] = fmincon(@(beta)lf_probit([ones(n,1),X],D,beta), zeros(k+1,1), [],[],[],[], -100, 100,[], options);
beta0 = opt_beta0(2:end)/opt_beta0(k+1,1); %creat starting values for the Klein & Spady estimator

%------------Implement the Klein-Spady estimator----------------------
options=optimset('MaxFunEvals', 1000, 'MaxIter', 1000,'TolX',1e-10,'TolFun',1e-10, 'Display', 'off' );
con = zeros(1,k); con(k)=1; % normalize the kth coeff. to 1
[beta] = fmincon(@(beta)lf(X,D,beta,opt,h), beta0, [], [],con, 1, -100, 100,[], options);
Xb = X*beta;
hD = 2.7799 * std(Xb) * n^(-1/5); % Optimal bandwidth for the Quartic kernel by Silverman's rule
margEff = mfx(X,D,beta,opt,hD);
density = margEff(k,1);


%%
    function value2 = lf_probit(X,D,beta)
        Xb = X*beta;
        F = normcdf(Xb);

        Lm=1e-323;
        F(F<=Lm) = Lm;
        nF = 1-F; nF(nF<=Lm) = Lm;

        value2 = -(D'*log(F) + (1-D)'*log(nF))/length(D);
        value2(isnan(value2)|isinf(value2)) = 1e+308;
    end
%%
    function value = lf(X,D,beta,opt,h)
        
        Xb = X*beta;
        n = length(D);
        
        if isempty(h)
            sigma = std(Xb);
            h = 2.7799 * sigma * n^(-1/5);  % Optimal bandwidth for Quartic Kernel
            %h = 2.3466 * sigma * n^(-1/5); % for Epanechnikov kernel
            %h = 1.06 * sigma * n^(-1/5);   % for Gaussian kernel
        end

        Kernel = @(u) 15/16 * (1-u.^2).^2.*(abs(u)<=1);
        F = zeros(n,1);
        
        if strcmp(opt, 'krll')
            for i=1:n
                dis = Xb - Xb(i);
                K = Kernel(dis/h);
                s1=dis'*K;
                s2=(dis.^2)'*K;
                F(i)=sum((s2-s1*dis).*K.*D)/(s2*sum(K)-s1^2);
            end
        else
            for i=1:n
                dis = Xb - Xb(i);
                K = Kernel(dis/h);
                F(i) = K'* D/sum(K);
            end
        end

        Flag = ~(F~=F);
        F(~Flag) = 1e+308;

        Lm=1e-323;
        F(F<=Lm) = Lm;
        nF = 1-F; nF(nF<=Lm) = Lm;

        value = -sum(Flag.*(D.*log(F) + (1-D).*log(nF)))/sum(Flag);
        value(value~=value) = 1e+308;
    end

%%
    function Margeff = mfx(X,D,beta,opt,h)
        Xb = X*beta;
        Xbar = mean(Xb);
        dis = Xb - Xbar;
        u = dis/h;
        K = 15/16*(1-u.^2).^2.*(abs(u)<=1);

        if strcmp(opt, 'krll')
            sw = sum(K);
            t1 = dis.*K;
            t2 = sum(t1);
            f =(sw*t1'*D - t2*K'*D)/(sw*t1'*dis - t2*t2);
        else
            dkdu = 15/4*u.*(u.^2-1).*(abs(u)<=1);
            f = -1/h*(dkdu'*D/sum(K)-K'*D*sum(dkdu)/sum(K)^2);
            if any(abs(u) == 1)
                f = nan;
            end
        end
        Margeff = f*beta;
    end

end