new; cls;            
/* Bootstrapped Kernel Density Monte Carlo from Small Data        */
/*                                                                */
/* Notice that the data consist of 'mean'of each bootstrap loop,  */
/* so density, its mean and variance of random draws here are all */
/* slightly different from what we call density and its mean and  */
/* variance.   (Wait for a moment....to see 2nd random draw plot) */
let P[30,1]=
17.44 17.48 17.72 17.67 17.40 17.37 17.72 17.72 17.52 17.88
18.32 18.73 18.69 18.65 18.10 18.39 18.39 18.24 17.95 18.09
18.39 18.52 18.54 18.78 18.59 18.46 18.30 18.24 18.46 18.27
;
L=ln(P);
data=L[2:rows(L)]-L[1:rows(L)-1];        
n=200;                /* Should increase more, but it will take some hours to run. */
times=2000;           /* number of bootstrap */
data=bootdata(data,times);
call rndkern(data,n);


proc bootdata(data,times);
   local nobs,bdata,i;
   nobs=rows(data);
   bdata=zeros(times,1);
   i=1;
   do while i<=times;
      bdata[i]=meanc(data[ceil(nobs*rndu(nobs,1))]);
      i=i+1;
   endo;
   retp(bdata);
endp;

/*
**  rndkern.txt - kernel density-random number generator from real data.
** (C) Copyright 2003 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Computes kernel density-random numbers under estimated kernel density.
**
**  Format:     y = rndkern(data,n);   
**                     or
**              call rndkern(data,n);
**
**  Input:      data     vector, a column vector of real data
**
**              n        scalar, number of rows of resulting vector
**
**  Output:     y        n x 1 vector, random numbers. 
*/

proc rndkern(data,n);
   local m,nn,low,up,maxf,X,i,Y,U;
   m=plotdens(data);
   nn=rows(m);
   low=m[1,1]; up=m[nn,1];
   maxf=maxc(m[.,2]);
   library pgraph;
   graphset;
   pqgwin auto;
   title("kernel density f(x)");
   xy(m[.,1],m[.,2]);
   print "maxf=" maxf;
   print/rz "[" low "," up "]";
   print;
   print "Proceeding...(wait for a long time)";
/* acceptance rejection method */
   X=zeros(n,1);
   i=1;
   do while i<=n;
      Y=(up-low)*rndu(1,1)+low;
      U=rndu(1,1);
      if U<=f(Y)/maxf;
         X[i]=Y;
         call sleep(0.1);
         print/lz i "out of " n;
         i=i+1;
      endif;
   endo;
   title("histogram of random draw");
   call hist(X,19);
   retp(X);
endp;

proc f(x);
   local m,fx,index;
   m=plotdens(data);
   fx=m[.,2];
   index=minindc(abs(m[.,1]-x));
   retp(fx[index]);
endp;

proc plotdens(x);
   local npoints,n,h,points,kern,i,j;
   npoints=51;                      /* Cut in less number of segments */
   n=rows(x);        
   x=sortc(x,1);
   h=1.059*stdc(x)*n^(-0.2)*1;
   points=seqa(x[1],(x[n]-x[1])/(npoints-1),npoints);
   i=1; kern=zeros(npoints,1);
   do while i<=npoints;
      j=1;
      do while j<=rows(h);
         kern[i,j]=gkernel(points[i],x,h[j]);
         j=j+1;
      endo;
      i=i+1;
   endo;
   retp(points~kern);
endp;

proc gkernel(x0,x,h);
   local u,kv;
   u=(x0-x)/h;
   kv=(1/sqrt(2*pi))*exp(-0.5*u^2);
   retp( sumc(kv/(rows(x)*h)) );
endp;