
/* This a collection of procedures to generate a Matrix n_ind (rows) and r (cols) of a Halton sequence with 
prime number X and after discarding "d" numbers.

We have:

halton3 that uses prime number 3
halton5 that uses prime number 5


For more information, read Kenneth Train "Discrete Choice Methods with Simulation" pp. 224-234

p must be a discrete number, eg, 2,3,5,7,11...

recomended "d" is 20 (see Cameron and Trivedi Microeconometrics Book

Notice, the result of:
g=halton3(2,5,10);
g;

coincides with the example of Kenneth Train page 227


Notice, to generate normal(0,1) from these Halton draws, you will use the inverse of the normal cdf, that is:

g=halton3(2,5,10);
normal_draws=cdfni(g);

*/

@----------------------Halton3-----------------------@

proc halton3(n_ind,r,d);

@Primero se definen las variables locales a la funcion de verosimilitud@
local t, halton_acu, a1, a2, a3;

t=0;
halton_acu=0;

@Notice, this loops produces the iterative formula in page 225 of the book@
do until rows(halton_acu)>n_ind*r+d;
  t=t+1;
  a1=halton_acu;
  a2=halton_acu+1/(3^t);
  a3=halton_acu+2/(3^t);
  halton_acu=a1|a2|a3;
endo;

@This is to discard the first d draws as the book recommends@
halton_acu=trimr(halton_acu,d,0);

@The following one deletes the excess of draws@
halton_acu=trimr(halton_acu,0,rows(halton_acu)-n_ind*r);

halton_acu=reshape(halton_acu,n_ind,r);

retp(halton_acu);
endp;


@----------------------Halton5-----------------------@

proc halton5(n_ind,r,d);

@Primero se definen las variables locales a la funcion de verosimilitud@
local t, halton_acu, a1, a2, a3, a4, a5;

t=0;
halton_acu=0;

@Notice, this loops produces the iterative formula in page 227 of the book@

do until rows(halton_acu)>n_ind*r+d;
  t=t+1;
  a1=halton_acu;
  a2=halton_acu+1/(5^t);
  a3=halton_acu+2/(5^t);
  a4=halton_acu+3/(5^t);
  a5=halton_acu+4/(5^t);
  halton_acu=a1|a2|a3|a4|a5;
endo;

@This is to discard the first d draws as the book recommends@
halton_acu=trimr(halton_acu,d,0);

@The following one deletes the excess of draws@
halton_acu=trimr(halton_acu,0,rows(halton_acu)-n_ind*r);

halton_acu=reshape(halton_acu,n_ind,r);

retp(halton_acu);
endp;



