
function f = frontier(m, x, xdata, z,zdata,q,qdata,ydata,nobs)
%m: expected max out of $m$ draws (scalar)
%x: value of input to evaluate frontier (scalar)
%xdata: data on input x (vector)
%z: value of input to evaluate frontier (scalar)
%zdata: data on input z (vector)
%q: value of input to evaluate frontier (scalar)
%qdata: data on input q (vector)
%ydata: data on output (vector)
%nobs: exclude $nobs$ largest values of y 
	ydata=ydata(xdata<=x&zdata<=z&qdata<=q);    
    ydata=sortrows(ydata);
    
  if(length(ydata)<=nobs)
        f=0;
  else
    ydata=ydata(1:max(1,length(ydata)-nobs));
    nx=length(ydata);
    ydata_t=ones(size(ydata));
         ft=(cumsum(ydata_t)./nx).^m;
         f=ydata(1)*ft(1);  
         if(length(ydata)>=2)
         f=f+sum(diff(ft).*ydata(2:end));         
         end
         
   end
end
	