program warpregf;
{First written on November 11, 1994; last revised November 11, 1994.
Authors: Isaias Hazarmabeth Salgado-Ugarte, Makoto Shimizu
and Toru Taniuchi.
This version do not display results to screen, instead write them to
a text file in order to be used by Stata's command "infile" as a part
of a Stata's ado file (warpreg.ado).
This Pascal program performs Nayarada-Watson Nonparametric regression
employing Haerdle's WARPing approximation using modified algorithms
taken from Haerdle (1991) "Smoothing Techniques with Implementations
in S. Springer-Verlag Series in Statistics New York}

const
  maxsize = 1000;
  maxsize2 = 1000; {it is a really big limit!}

type
  vector = array[0..maxsize] of real;
  vectorint = array[0..maxsize2] of integer;

var
  indexstart,indexend,nl,n,i,M,numbin,kerneltype,j   : integer;
  start, delta,h,origin,maxval,minval                : real;
  bin,frequ,index                                    : vectorint;
  x,y,ysum,kwe,fm,rm,mm,midval                       : vector;
  inpval,data,resul                                  : text;


procedure InputValues;
begin
   assign(data,'_data2.raw');
   reset(data);
   i:=1;
   while (not eof(data)) do
    begin
       read(data,y[i],x[i]);
       i:=i+1;
    end;
   close(data);
   n:=i-2;
   assign(inpval,'_inpval.raw');
   reset(inpval);
   read(inpval,h,M,kerneltype);
   close(inpval);
end;


procedure InitialCalc;
var maxnoemp,ia,ib: integer;

begin
   maxval:=x[n];
   minval:=x[1];
   if kerneltype=6 then
      h:=4*h;
   delta:= h/M;
   start:=minval-h-delta*0.1;
   if (start<0) then
      origin:=(round((start/delta)-0.5)-0.5)*delta
   else
      origin:=(trunc(start/delta)-0.5)*delta;
   numbin:=trunc(((maxval+h+delta*0.1-origin)/delta)+0.5)+1;
   if (n+1<numbin) then
      maxnoemp:=n+1
   else
      maxnoemp:=numbin;
   indexstart:=trunc((minval-origin)/delta);
   indexend:=trunc((maxval-origin)/delta);
   for ia:=indexstart to indexend+1 do
         midval[ia]:=origin+((0.5+ia)*delta);
   nl:=0;
   j:=0;
   for ib:= 0 to numbin+1 do
      begin
        bin[ib]:=0;
        frequ[ib]:=0;
        index[ib]:=0;
      end;
end;

procedure BinningData;

begin
   for i:=1 to n do
    begin
      j:= trunc((x[i]-origin)/delta);
      if (bin[j]=0) then
         begin
            nl:= nl+1;
            index[nl]:=j;
            frequ[nl]:=1;
            ysum[nl]:=y[i];
            bin[j]:=nl;
         end
      else
         begin
           frequ[bin[j]]:=frequ[bin[j]]+1;
           ysum[bin[j]]:=ysum[bin[j]]+y[i];
         end
    end;
end;

procedure CreateWeight; {Note: weights are not normalized}
var
l: integer;

function Power(x:real; n:integer): real;
  begin
  if n=1 then
    Power:=x
  else if odd(n) then
    Power:=x*sqr(Power(x, n div 2))
  else
    Power:=sqr(Power(x, n div 2))
  end;

begin
  case kerneltype of
    1: begin      {uniform kernel}
         for l:=0 to M-1 do
           kwe[l]:= 1;
       end;
    2: begin {Triangle kernel (Original Average-Shifted-Histogram kernel)}
         for l:=0 to M-1 do
           kwe[l]:= 1 -l/M;
        end;
    3: begin     {Epanechnikov kernel}
         for l:=0 to M-1 do
           kwe[l]:= 1-Power((l/M),2);
       end;
    4: begin   {Quartic Kernel}
         for l:=0 to M-1 do
           kwe[l]:=Power((1-Power((l/M),2)),2);
       end;
    5: begin    {Triweight kernel}
         for l:=0 to M-1 do
           kwe[l]:=Power((1-Power((l/M),2)),3);
       end;
     6: begin   {Gaussian kernel}
          for l:=0 to M-1 do
           kwe[l]:=exp(-8.*(Power((l/M),2)));
        end;
    end;{case}
end;

procedure WeightBins;
var
kw,lw,iw,indexi,lwabs  : integer;

begin
  for iw:=0 to numbin do
     begin
        fm[iw]:=0;
        rm[iw]:=0;
     end;
  for kw:=1 to nl do
      for lw:=1-M to M-1 do
        begin
          indexi:=lw+index[kw];
          lwabs:=abs(lw);
          fm[indexi]:= fm[indexi]+kwe[lwabs]*frequ[kw];
          rm[indexi]:= rm[indexi]+kwe[lwabs]*ysum[kw];
        end;
end;

procedure CombiNumDen;
begin
   for j:=0 to numbin-1 do
      if fm[j]=0 then
         mm[j]:=0
      else
         mm[j]:=rm[j]/fm[j];
end;

procedure ResulFile;
var ir: integer;
begin
  assign(resul,'resfile');
  rewrite(resul);
  for ir:=indexstart to indexend do
    writeln(resul,mm[ir]:12:8,midval[ir]:12:6);
  close(resul);
end;

begin {main}
    InputValues;
    InitialCalc;
    BinningData;
    CreateWeight;
    WeightBins;
    CombiNumDen;
    ResulFile;
end.