function [beta, margEff, deriv_g] = Ichimura(X,Y,h)

% X: n-by-k matrix
% Y: n-by-1 vector
% h: a positive scalar

% [beta, margEff, deriv_g] = Ichimura(X,Y,h) implements the Ichimura estimator for a semiparametric
% regression model Y = g(X*b)+e, using kernel regression with the Quartic kernel.
% The Ichimura estimator can only identify beta up to location and scale, so the first coefficient is normalized to 1.

% Outputs allowed: coefficients, marginal effects, and the derivative of g function evaulated at the index mean, g'(Xbar*b). 

% beta = Ichimura(X,Y) implements the Ichimura estimator, using the optimal bandwidth computed by Silverman's rule.
% beta = Ichimura(X,Y,h) implements the Ichimura estimator, using the user provided bandwidth, h.


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

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

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

[n,k]= size(X);

%---------Obtain starting values of coeff.s in Y = g(X'b)+ e ----------------
% Let's use OLS estimates as starting values
X1 = [ones(n,1),X];
opt_beta0 = (X1'*X1)\(X1'*Y); 
beta0 = opt_beta0(2:end)/opt_beta0(2,1); % Need to remove the constant term and do scale normalization, here I assume the first coeff. of b is 1 

%------------Implement the Klein-Spady estimator----------------------
options=optimset('MaxFunEvals', 1000, 'MaxIter', 1000,'TolX',1e-10,'TolFun',1e-10, 'Display', 'iter' );
con = [1,zeros(1,k-1)]; % Normalize the first coeff. to 1
[beta] = fmincon(@(beta)Objfun(X,Y,beta,h), beta0, [], [],con, 1, -100, 100,[], options);
Xb = X*beta;
hY = 2.7799 * std(Xb) * n^(-1/5); % Optimal bandwidth for the Quartic kernel by Silverman's rule
margEff = mfx(X,Y,beta,hY);
deriv_g = margEff(1);


%%
    function value = Objfun(X,Y,beta,h)
        
        Xb = X*beta;
        n = length(Y);

        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);

        for i=1:n
            dis = Xb - Xb(i);
            K = Kernel(dis/h);
            F(i) = K'* Y/sum(K);
        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.*(Y-F))/sum(Flag);
        value(value~=value) = 1e+308;
    end

%%
    function Margeff = mfx(X,Y,beta,h)
        Xb = X*beta;
        Xbar = mean(Xb);
        dis = Xb - Xbar;
        u = dis/h;
        K = 15/16*(1-u.^2).^2.*(abs(u)<=1);
        dkdu = 15/4*u.*(u.^2-1).*(abs(u)<=1);
        f = -1/h*(dkdu'*Y/sum(K)-K'*Y*sum(dkdu)/sum(K)^2);
        if any(abs(u) == 1)
            f = nan;
        end

        Margeff = f*beta;
    end

end