new; cls;            
/* Technical note and correction!                                                              */
/* Emprical Monte Carlo using a liner interpolation is completely WRONG in a theoretical sense */
/* because we cannot convert the ordinal statistic to cardinal one. That's why this method is  */
/* so simple but not implemented. Especiallly, it is obvious if we consider the updated "cdf". */
/* There are several papers and a book in Japan and a software in US dealing with some variant */
/* of this, but all should be dropped. Especially, they should be dropped academically. That's */
/* partially my responsibility though I do not argue which one is first to come. Academically  */
/* worng is nothing but wrong. All we can do is to blur the realized ordinal points(Option 0). */

/* Random Number Generation by Inverse of Empirical CDF */
/*       Quasi-Random Sequence (van der Corput) version */
/* Data from Sheldon M. Ross(1999) Table 10.5 in        */
/* "An Introduction to Mathematical Finance"            */
let P[752,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
18.32 18.42 18.59 18.91 18.91 18.86 18.63 18.43 18.69 18.66
18.49 18.32 18.35 18.63 18.59 18.63 18.33 18.02 17.91 18.19
17.94 18.11 18.16 18.26 18.56 18.43 18.96 18.92 18.78 19.07
19.05 19.22 19.15 19.17 19.03 19.18 19.56 19.77 19.67 19.59
19.88 19.55 19.15 19.15 19.73 20.05 20.41 20.52 20.41 20.12
20.29 20.15 20.43 20.38 20.50 20.09 19.89 20.29 20.33 20.29
19.61 19.75 19.41 19.52 19.90 20.08 19.96 20.00 20.06 19.81

19.77 19.41 19.26 18.69 18.69 18.78 18.89 18.90 19.14 19.25
19.06 19.18 18.91 18.80 18.86 18.91 19.05 18.94 18.84 18.22
18.01 17.46 17.50 17.49 17.64 17.77 17.97 17.56 17.40 17.40
17.40 17.18 17.37 17.14 17.34 17.32 17.49 17.25 17.32 17.20
17.35 17.33 17.01 16.79 16.88 16.93 17.50 17.49 17.43 17.56
17.70 17.78 17.72 17.71 17.65 17.79 17.78 17.89 17.86 17.48
17.47 17.55 17.66 17.87 18.25 18.54 18.00 17.86 17.86 17.82
17.82 17.79 17.84 18.04 18.04 18.58 18.36 18.27 18.44 18.47
18.64 18.54 18.85 18.92 18.93 18.95 18.69 17.56 17.25 17.47
17.33 17.57 17.76 17.54 17.64 17.56 17.30 16.87 17.03 17.31

17.42 17.29 17.12 17.41 17.59 17.68 17.61 17.32 17.37 17.21
17.32 17.32 17.58 17.54 17.62 17.64 17.74 17.98 17.94 17.71
17.65 17.82 17.84 17.83 17.80 17.82 17.93 18.19 18.57 18.06
17.97 17.96 17.96 17.96 18.38 18.33 18.26 18.18 18.43 18.63
18.67 18.77 18.73 18.97 18.66 18.73 19.00 19.11 19.51 19.67
19.12 18.97 18.96 19.14 19.14 19.27 19.50 19.36 19.55 19.55
19.81 19.89 19.91 20.26 20.26 19.95 19.67 18.79 18.25 18.38
18.05 18.52 19.18 18.94 18.62 18.06 18.28 17.67 17.73 17.45
17.56 17.74 17.71 17.80 17.54 17.69 17.74 17.76 17.78 17.97
18.91 18.96 19.04 19.16 19.16 21.05 19.71 19.85 19.06 19.39

19.70 19.29 19.54 19.44 19.20 19.54 20.19 19.81 19.61 19.91
20.46 20.58 21.16 21.99 23.27 24.34 23.06 21.05 21.95 22.40
22.19 21.79 21.41 21.47 22.26 22.70 22.27 22.75 22.75 23.03
23.06 24.21 25.34 24.29 25.06 24.47 24.67 23.82 23.95 24.07
22.70 22.40 22.20 22.32 22.43 21.20 20.81 20.86 21.18 21.04
21.11 21.00 20.68 21.01 21.36 21.42 21.48 20.78 20.64 22.48
22.65 21.40 21.23 21.32 21.32 21.11 20.76 19.94 19.76 19.85
20.44 19.72 20.05 20.28 20.25 20.10 20.09 20.01 20.34 22.14
21.46 20.76 20.65 19.92 19.98 19.96 20.65 21.02 20.92 21.53
21.13 21.21 21.21 21.21 21.27 21.41 21.55 21.95 21.90 22.48

22.38 21.80 21.68 21.00 21.40 21.01 20.68 20.74 20.11 20.28
20.33 20.42 21.04 21.34 21.23 21.13 21.42 21.55 21.57 22.22
22.37 22.12 21.90 22.66 23.26 22.86 21.72 22.30 21.96 21.62
21.56 21.71 22.15 22.25 22.25 23.40 23.24 23.44 23.85 23.73
24.12 24.75 25.00 24.51 23.19 23.31 23.89 23.54 23.63 23.37
24.07 24.46 24.16 24.60 24.38 24.14 24.05 24.81 24.73 25.24
25.54 25.07 24.26 24.66 25.62 25.42 25.17 25.42 25.75 25.92
25.75 24.86 24.51 24.86 24.85 24.34 24.28 23.35 23.03 22.79
22.64 22.69 22.74 23.59 23.37 23.35 24.12 24.41 24.17 23.88
24.49 23.76 23.84 23.75 23.49 23.62 23.75 23.75 23.75 24.80

24.93 24.80 25.58 25.62 25.30 24.42 23.38 23.72 24.47 25.74
25.71 26.16 26.57 25.08 24.79 25.10 25.10 24.92 25.22 25.37
25.92 25.92 25.69 25.59 26.37 26.23 26.62 26.37 26.09 25.19
25.11 25.95 25.52 25.41 25.23 24.80 24.24 24.18 24.05 23.94
23.90 24.47 24.87 24.15 24.15 24.02 23.91 23.10 22.23 22.46
22.42 21.86 22.02 22.41 22.41 22.52 22.79 21.98 21.39 20.71
21.00 21.11 20.89 20.30 20.25 20.66 20.49 20.94 21.28 20.49
20.11 20.62 20.70 21.29 20.92 22.06 22.04 22.32 21.51 21.06
20.99 20.64 20.70 20.70 20.41 20.28 19.47 19.47 19.12 19.23
19.35 19.27 19.57 19.53 19.90 19.83 19.35 19.42 19.91 20.38

19.60 19.73 20.03 19.99 19.91 20.44 20.21 19.91 19.60 19.63
19.66 19.62 20.34 20.43 21.38 21.37 21.39 21.30 22.12 21.59
21.19 21.86 21.86 21.63 21.63 20.79 20.79 20.97 20.88 21.12
20.33 20.12 19.66 18.79 18.68 18.67 18.53 18.69 18.83 19.01
19.23 18.79 18.67 18.55 19.14 19.03 19.52 19.09 19.46 19.80
20.12 20.34 19.56 19.56 19.52 19.73 19.46 19.22 19.33 18.99
19.67 19.65 19.99 19.27 19.18 19.08 19.63 19.77 19.89 19.81
19.85 20.30 20.14 20.28 20.75 20.81 20.46 20.09 19.54 19.69
19.99 20.19 20.08 20.07 19.91 20.12 20.06 19.66 19.70 19.26
19.28 19.73 19.58 19.61 19.61 19.65 19.61 19.40 19.63 19.45

19.42 19.42 19.37 19.32 19.27 19.61 19.42 19.38 19.35 19.60
19.79 19.94 20.39 20.87 21.26 21.18 21.05 21.77 22.76 21.93
21.96 22.18 22.12 22.10 21.32 20.70 20.57 20.97 20.59 20.70
20.67 21.42 21.09 20.97 21.07 20.46 20.71 21.22 21.08 20.96
20.70 20.31 20.39 20.77 20.40 20.51 20.49 20.70 21.00 20.26
20.04 19.80
;
L=ln(P);
data=L[2:rows(L)]-L[1:rows(L)-1];

nn=10000;
library pgraph;
graphset;
call hist(rndempinv3(data,nn,0),49);


/*
**  rndempinv3.txt - Random number generation by inverse of empirical cdf from real data.
**                   (Quasi-Random Sequence (van der Corput) version)
** (C) Copyright 2005 Yosuke Amijima. All Rights Reserved.
**
**  Purpose:    Generates random numbers by inverse of empirical cdf.
**
**  Format:     y = rndempinv3(data,nn,linetype);   
**
**  Input:      data     n x 1 vector,    real data
**
**              nn       scalar,          number of resulting random numbers
**
**              linetype 0: original step   1: linear interpolation(approximate)
**
**  Output:     y        nn x 1 vector,   random numbers
**
*/
proc rndempinv3(data,nn,linetype);
   local bb,x,i,n0,ff,n1,r,n,p,y,j,index;
/* Quasi-Random Sequence (van der Corput) */
   bb=2;
   x=zeros(nn,1);
   i=1;
   do while i<=nn;
      n0=i;
      ff=1/bb;
      do while n0>0;
         n1=floor(n0/bb);
         r=n0-n1*bb;
         x[i]=x[i]+ff*r;
         ff=ff/bb;
         n0=n1;
      endo;
      i=i+1;
   endo;
/* main */
   data=sortc(data,1);
   n=rows(data);
   p=seqa(0,1,n)/(n-1);
   y=zeros(nn,1);
   j=1;
   do while j<=nn;
      if x[j]==1;
         y[j]=maxc(data);
      else;
         index=sumc(x[j].>=p);
         if linetype==0;
            y[j]=data[index]+(data[index+1]-data[index])*rndu(1,1); 
         else;
            y[j]=data[index]+(data[index+1]-data[index])*(x[j]-p[index])/(p[index+1]-p[index]);
         endif;
      endif;
      j=j+1;
   endo;
   retp(y);
endp;