program warpings(input,output);
{First written on August, 1994; last correction September 12, 1994.
Authors: Isaias Hazarmabeth Salgado-Ugarte, Makoto Shimizu and
Toru Taniuchi. Reference: Stata Technical Bulletin XX, snp6.1

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 three Stata's ado files (warpstep.ado, warpoly.ado and warping.ado).
This Pascal program performs Haerdle's WARPing density estimation
using modified algorithms from Haerdle (1991) "Smoothing Techniques
with Implementations in S. Springer-Verlag Series in Statistics
New York, and Scott (1992) "Multivariate Density Estimation: Theory,
Practice, and Visualization. John Wiley and Sons, New York.}

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

type
  vector = array[1..maxsize] of real;
  striname = string[30];
  vectorint = array[1..maxsize2] of integer;

var
  binumber,n,i,M,numbin,kerneltype,index   : integer;
  start, delta,h,origin,maxval,minval,cm   : real;
  binmesh,counts,J                         : vectorint;
  x,wM,fh,midval                           : vector;
  ch,answer1,answer3                       : char;
  name                                     : striname;
  data,resul                               : text;

procedure DataInput;
begin
   writeln;
   writeln('TYPE THE PATH, NAME AND EXTENSION OF TEXT DATA FILE');
   readln(name);
   assign(data,name);
   reset(data);
   i:=1;
   while (not eof(data)) do
     begin
        read(data,x[i]);
        i:=i+1;
     end;
   close(data);
   n:=i-2;
   writeln;
   writeln('THE NUMBER OF VALUES READ IS: ',n:10);
   writeln;
   writeln('!!!!!WARNING!!!!!');
   writeln('IF THIS IS NOT CORRECT PLEASE INTERRUPT AND');
   writeln('USE STATA COMMAND outfile TO GENERATE AN ASCII');
   writeln('FILE WITH THE DESIRED DATA VECTOR');
   writeln;
   writeln('GIVE THE VALUE OF THE BANDWIDTH ''h'' ');
   readln(h);
   writeln;
   writeln('GIVE THE NUMBER OF HISTOGRAMS TO SHIFT AND AVERAGE');
   readln(M);
   writeln;
end;

procedure SelectKernel;
begin
   kerneltype:=0;
   writeln;
   writeln('SPECIFY THE WEIGHT FUNCTION: ');
   writeln;
   writeln('1 = Uniform; 2 = Triangle (ASH); 3 = Epanechnikov');
   writeln('4 = Quartic; 5 = Triweight; 6 = Gaussian');
   readln(kerneltype);
   if (kerneltype>6) then
     begin
       writeln;
       writeln('Error in the choice of the kernel! Use a valid option');
       writeln;
       writeln('TO TRY AGAIN PUSH ENTER');
       readln;
       SelectKernel;
     end;
end;

procedure SortData;
var
   loc   : integer;
   temp  : real;

begin
   for loc:=1 to n-1 do
      for i:=loc+1 to n do
         if x[i]<x[loc] then
            begin
              temp:=x[loc];
              x[loc]:=x[i];
              x[i]:=temp;
            end;
end;

procedure InitialCalc;
var maxnoemp,ib: integer;

begin
   maxval:=x[n];
   minval:=x[1];
   if (kerneltype=6) then
      h:=4*h;
   delta:= h/M;
   numbin:= trunc((maxval-minval)/delta)+2*(M+1+round((M/10)+0.5));
   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;
   if (n+1<numbin) then
      maxnoemp:=n+1
   else
      maxnoemp:=numbin;
   binumber:=0;
   index:=0;
   for ib:= 0 to numbin+1 do
      begin
        binmesh[ib]:=0;
        counts[ib]:=0;
        J[ib]:=0;
      end;
end;

procedure BinmeshCalc;
var
   ia,bindx,bi,totcount : integer;
   sum1, sum2,interval : real;

begin
   for i:=1 to n do
    begin
      index:= trunc((x[i]-origin)/delta);
      if (binmesh[index]=0) then
         begin
            binumber:= binumber+1;
            J[binumber]:=index;
            counts[binumber]:=1;
            binmesh[index]:=binumber;
          end
      else
        counts[binmesh[index]]:=counts[binmesh[index]]+1;
    end;
    totcount:=0;
    for ia:=0 to numbin-1 do
         midval[ia]:=((0.5+ia)*delta)+origin;
end;

procedure CreateWeight;
var
l: integer;
part1,part2 : real;

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

begin
  case kerneltype of
    1: begin      {uniform kernel}
         cm:=M/((2*M-1)*(n*h));
         for l:=0 to M-1 do
           wM[l]:= cm;
       end;
    2: begin {Triangle kernel (Original Average-Shifted-Histogram kernel)}
         cm:=1/(n*h);
         for l:=0 to M-1 do
           wM[l]:= cm*(1 -l/M);
        end;
    3: begin     {Epanechnikov kernel}
         cm:=3*sqr(M)/((4*sqr(M)-1)*(n*h));
         for l:=0 to M-1 do
           wM[l]:= cm*(1-sqr(l/M));
       end;
    4: begin   {Quartic Kernel}
         cm:=0.9375/((1-0.0625/(sqr(M)*sqr(M)))*(n*h));
         for l:=0 to M-1 do
           wM[l]:=cm*sqr(1-sqr(l/M));
       end;
    5: begin    {Triweight kernel}
         part1:=1+0.14583333/CleverPower(M,4);
         part2:=0.052083333/CleverPower(M,6);
         cm:= 1.09375/((part1-part2)*(n*h));
         for l:=0 to M-1 do
           wM[l]:=cm*CleverPower(1-sqr(l/M),3);
       end;
     6: begin   {Gaussian kernel}
          cm:= 0.3989*4/(n*h);
          for l:=0 to M-1 do
           wM[l]:=cm*exp(-8.0*(sqr(l/M)));
        end;
    end;{case}
end;

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

begin
  for iw:=0 to numbin-1 do
     fh[iw]:=0;
  for kw:=1 to binumber do
      for lw:=1-M to M-1 do
        fh[J[kw]+lw]:= fh[J[kw]+lw]+wM[abs(lw)]*counts[kw];
end;

procedure ResulFile;
var rcount: integer;
begin
  assign(resul,'resfile');
  rewrite(resul);
  for rcount:=0 to numbin-1 do
    writeln(resul,fh[rcount]:12:8,midval[rcount]:12:6);
  close(resul);
end;

begin {main}
    DataInput;
    SelectKernel;
    SortData;
    InitialCalc;
    BinmeshCalc;
    CreateWeight;
    WeightBins;
    ResulFile;
end.