function [est_vc,est_qr] = vcoef (outdata, toler, beta, q, h1,h2)

% this program used to estimate varying-coefficient function for both of
% mean regression and quantile regression; 
% model is Y=a(U)X+epsilon 
% for example
% p=size(outdata,2)-1 
% [est_vc,est_qr] = vcoef (outdata, 10^(-6), zeros(2*p,1), 0.5, 0.2,0.3)

%  input 
%     -- outdata=[Y,X,U], X=[X1,X2,...,Xp]
%     -- h1 is the bandwidth for mean regression and h2 for quantile regression.

%     -- q is the quantile desired (number in (0,1)).
%     -- toler is a small number giving the minimum change in the value
%        of the surrogate function before convergence is declared, used for quantile regression.
%     -- beta is the starting value of beta for quantile regression.
%  output 
%     -- est_vc return estimator for mean regression with vc in section2
%     -- est_qr return estimator for quantile regression with vc in section5
  

Y=outdata(:,1); 
U=outdata(:,end);
X=outdata(:,2:end-1); 
[n,p]=size(X);
XX=[ones(n,1),X];
p=p+1;


f=@(x) 3*(1-x.^2).*(abs(x)<=1)/4;  % Epanechnikov kernel function 
[n,m]=size(outdata);
func=@(XX,beta)(XX*beta);    
tn=toler/n; e0=-tn/log(tn); epsilon=(e0-tn)/(1+log(e0));
   %  epsilon is the smoothing parameter.  As suggested in the Hunter and Lange's paper, it
   %  is chosen to roughly satisfy -epsilon(log epsilon)=toler/n using a Newton-Raphson step above.

u=sort(U);
est_vc=[];  % save estimator for mean regression 
est_qr=[];  % save estimator for quantile regression 

for i=1:length(u)
    
    
    i
    u(i)
    if i>=2 & u(i)==u(i-1)
        ;
    else
        
    repu=repmat(U-u(i),1,p);
    ZX=[XX,XX.*repu];
    
    % ++++++++++++++++++++++++++++++++++++++++++++++++++
    % local kernal smooth estimator for mean regression 
    % (Fan & Gijbels,1996) in section 2 
    % ++++++++++++++++++++++++++++++++++++++++++++++++++
    Ker=diag(f((U-u(i))/h1));
    est0=pinv(ZX'*Ker*ZX)*(ZX'*Ker*Y);         
    est_vc=[est_vc;[u(i),est0(1:p)']];                 

    % ++++++++++++++++++++++++++++++++++++++++++++++++++
    % to estimate quanitle regreesion with vc in section 5 
    % by MM algorithm (Hunter and Lange, 2000)
    % ++++++++++++++++++++++++++++++++++++++++++++++++++
    iteration=0;
    change=realmax;     %parameters for controlling convergence
    Kh=f((U-u(i))/h2);
   
       %updating  res weights v and lastsurr
    res = Y -feval(func,ZX,beta);    % Compute residuals for current beta.
    weights = 1./(epsilon+abs(res)); 
    v = 1-2*q - res.*weights;
    lastsurr = -sum(Kh.*res.*v) - (1-2*q)*sum(Kh.*res); % Last value of surrogate (up to a constant), multiple 1/4 

    while (change > toler)
	   iteration = iteration + 1; 
       J=ZX;   %first difference f(x,beta) 
       WW=diag(Kh.*weights);
            %step = -inv(J'*(weights(:,ones(p,1)).*J)) * (J' * v);

       J2=J'*WW*J;
       J2=double(J2);

       J1=J' *(Kh.* v);
       J1=double(J1);
 
       step=-pinv(J2)*J1;       
       beta=beta+step;
        
       res=Y-feval(func,ZX,beta);
       surr = sum(Kh.*res.^2.*weights) + (4*q-2) * sum(Kh.*res);   % the value of surrogate(beta+step|beta)
    
     % Now do the step-halving procedure to ensure a decrease in the value of the surrogate function.    
       while (surr > lastsurr)
       		step=step/2;
       		beta=beta-step;
       		res=Y-feval(func,ZX,beta);
       		surr = sum(Kh.*res.^2.*weights) + (4*q-2) * sum(Kh.*res);
       end
       
     change=lastsurr-surr;
	 weights = 1./(epsilon+abs(res));
	 v = 1-2*q - res.*weights;
     lastsurr = -sum(Kh.*res.*v) - (1-2*q)*sum(Kh.*res);
    end    % end for while 
    
    est_qr=[est_qr;[u(i),beta(1:p)']];
    
    end  % end for if-else 
    
end